0001 function vb_plot_filter_response(ftype,fcut,fs,Norder,fmax,tmax)
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016
0017
0018
0019
0020
0021
0022
0023
0024 if nargin<4,
0025 help vb_plot_filter_response
0026 error('input arguments is wrong')
0027 end
0028
0029 Nfreq = 2000;
0030 f = (0:Nfreq)/Nfreq;
0031 dt = 1/fs;
0032
0033 switch ftype
0034 case 'low'
0035 Fmax = min(fcut(1)*3,fs/2);
0036 Tmax = 10/fcut(1);
0037 Tlim = 3/fcut(1);
0038 case {'high', 'high_exp'}
0039 Fmax = min(fcut(1)*5,fs/2);
0040 Tmax = 10/fcut(1);
0041 Tlim = 3/fcut(1);
0042 case {'band', 'stop'}
0043 Fmax = min(fcut(2)*3,fs/2);
0044 Tmax = 10;
0045 Tlim = 5/fcut(1);
0046 end
0047
0048 if exist('fmax','var'), Fmax = fmax; end;
0049 if exist('tmax','var'), Tlim = fmax; end;
0050 if Fmax > fs/2, Fmax = fs/2; end;
0051
0052 t0 = fix(0.5*Tmax/dt);
0053 t = 0:dt:Tmax;
0054 f = f*Fmax;
0055 X = zeros(length(t),1);
0056 X(t0) = 1;
0057
0058 if Norder <= 0
0059 switch ftype
0060 case 'low'
0061 [X ,B] = eegfilt(X', fs, 0, fcut(1) );
0062 case 'high'
0063 [X ,B] = eegfilt(X', fs, fcut(1), 0 );
0064 case 'band'
0065 [X ,B] = eegfilt(X', fs, fcut(1), fcut(2) );
0066 case 'stop'
0067 if Norder==0, Norder = fix(3*(fs/(fcut(2)-fcut(1)))); end;
0068 [B, A] = fir1( abs(Norder), fcut./(fs/2), 'stop');
0069 X = filtfilt(B, A, X);
0070 end
0071 A = 1;
0072 H = freqz(B,A,f,fs);
0073 elseif strcmp(ftype,'high_exp')==1
0074
0075 a = vb_lowpass_init(fs,fcut(1),1);
0076 A = [1 -a];
0077 B = 1-a;
0078 H = 1 - freqz(B,A,f,fs);
0079 X = vb_online_highpass_cut(X, fs, fcut(1));
0080 else
0081 switch ftype
0082 case 'low'
0083 [B,A] = butter(Norder, fcut(1)/(fs/2) ,'low');
0084 case 'high'
0085 [B,A] = butter(Norder, fcut(1)/(fs/2) ,'high');
0086 case 'band'
0087 [B,A] = butter(Norder, fcut(1:2)/(fs/2) );
0088 case 'stop'
0089 [B,A] = butter(Norder, fcut(1:2)/(fs/2) , 'stop');
0090 end
0091 H = freqz(B,A,f,fs);
0092 X = filtfilt(B,A, X );
0093 end
0094
0095 subplot 221
0096 plot(f,abs(H))
0097 xlabel('Freq[Hz]')
0098 title('Gain of frequency response')
0099 ylim([ 0 1.1])
0100 xlim([f(1) f(end)])
0101 hold on
0102
0103 subplot 222
0104 plot(f,angle(H))
0105 xlabel('Freq[Hz]')
0106 title('Angle of frequency response')
0107 hold on
0108 xlim([f(1) f(end)])
0109
0110 subplot 223
0111 plot(t,X)
0112 xlabel('Time[sec]')
0113 title('Impulse response')
0114 hold on
0115
0116 xlim([Tmax/2-Tlim, Tmax/2+Tlim])
0117
0118 subplot 224
0119 plot(B)
0120 hold on
0121 plot(A,'--')
0122 xlabel('Tap number')
0123 title('Filter coefficient')
0124 xlim([1 max(length(A),length(B))])
0125 fprintf('# of A and B = %d, %d\n',length(A),length(B))