vb_plot_filter_response

 plot filter response in frequency and Impulse response
  vb_plot_filter_response(ftype,fcut,fs,Norder,fmax,tmax,online)
 --- Input (required)
  ftype ; filter type
    'low'  : Lowpass  filter
    'high' : Highpass filter
    'band' : Bandpass filter
    'stop' : Stopband filter

  fs   : Sampling frequency [Hz]

  fcut : Cutoff frequency [Hz]
       = [f1]    Lowpass/Highpass
       = [f1 f2] Bandpass/Stopband

  Norder = 0: FIR filter (eegfilt)
  Norder > 0: filter order of Butterworth filter

 --- Input (optional)
 fmax : Max frequency in frequency plot [Hz]
        if it is empty or not given, it is automatically determined
 tmax : Max time in Impulse plot [sec]
        if it is empty or not given, it is automatically determined

 online = 1: online filter with time delay
        = 0: 'filtfilt' is used for filter with no time delay 
        if this is not given, 'filtfilt' is used

  2012-9-19 Masa-aki Sato
 

 Copyright (C) 2011, ATR All Rights Reserved.
 License : New BSD License(see VBMEG_LICENSE.txt)

0001 function    vb_plot_filter_response(ftype,fcut,fs,Norder,fmax,tmax,online)
0002 % plot filter response in frequency and Impulse response
0003 %  vb_plot_filter_response(ftype,fcut,fs,Norder,fmax,tmax,online)
0004 % --- Input (required)
0005 %  ftype ; filter type
0006 %    'low'  : Lowpass  filter
0007 %    'high' : Highpass filter
0008 %    'band' : Bandpass filter
0009 %    'stop' : Stopband filter
0010 %
0011 %  fs   : Sampling frequency [Hz]
0012 %
0013 %  fcut : Cutoff frequency [Hz]
0014 %       = [f1]    Lowpass/Highpass
0015 %       = [f1 f2] Bandpass/Stopband
0016 %
0017 %  Norder = 0: FIR filter (eegfilt)
0018 %  Norder > 0: filter order of Butterworth filter
0019 %
0020 % --- Input (optional)
0021 % fmax : Max frequency in frequency plot [Hz]
0022 %        if it is empty or not given, it is automatically determined
0023 % tmax : Max time in Impulse plot [sec]
0024 %        if it is empty or not given, it is automatically determined
0025 %
0026 % online = 1: online filter with time delay
0027 %        = 0: 'filtfilt' is used for filter with no time delay
0028 %        if this is not given, 'filtfilt' is used
0029 %
0030 %  2012-9-19 Masa-aki Sato
0031 %
0032 %
0033 % Copyright (C) 2011, ATR All Rights Reserved.
0034 % License : New BSD License(see VBMEG_LICENSE.txt)
0035 
0036 if nargin<4, 
0037     help vb_plot_filter_response
0038     error('input arguments is wrong')
0039 end
0040 
0041 Nfreq = 2000;
0042 f  = (0:Nfreq)/Nfreq;
0043 dt = 1/fs;
0044 
0045 switch    ftype
0046 case    'low'
0047     Fmax = min(fcut(1)*3,fs/2);
0048     Tmax = 10/fcut(1);
0049 case    {'high', 'high_exp'}
0050     Fmax = min(fcut(1)*5,fs/2);
0051     Tmax = 10/fcut(1);
0052 case    {'band', 'stop'}
0053     Fmax = min(fcut(2)*3,fs/2);
0054     Tmax = 10;
0055 end
0056 
0057 if exist('fmax','var')&&~isempty(fmax), Fmax = fmax; end;
0058 if Fmax > fs/2, Fmax = fs/2; end;
0059 
0060 t0 = fix(0.5*Tmax/dt);
0061 t  = 0:dt:Tmax;
0062 f  = f*Fmax;
0063 X  = zeros(length(t),1);
0064 X(t0) = 1;
0065 t  = t - t(t0);
0066 
0067 if Norder <= 0
0068     switch    ftype
0069     case    'low'
0070         [X ,B] = eegfilt(X', fs, 0, fcut(1) );
0071     case    'high'
0072         [X ,B] = eegfilt(X', fs, fcut(1), 0 );
0073     case    'band'
0074         [X ,B] = eegfilt(X', fs, fcut(1), fcut(2) );
0075     case    'stop'
0076         if Norder==0, Norder = fix(3*(fs/(fcut(2)-fcut(1)))); end;
0077         [B, A] = fir1( abs(Norder), fcut./(fs/2), 'stop');
0078         X = filtfilt(B, A, X);
0079     end
0080     A = 1;
0081     H = freqz(B,A,f,fs);
0082 elseif strcmp(ftype,'high_exp')==1
0083     % Exponential online highpass filter
0084     a = vb_lowpass_init(fs,fcut(1),1);
0085     A = [1 -a];
0086     B = 1-a;
0087     H = 1 - freqz(B,A,f,fs);
0088     X = vb_online_highpass_cut(X, fs, fcut(1));
0089 else
0090     switch    ftype
0091     case    'low'
0092         [B,A] = butter(Norder, fcut(1)/(fs/2) ,'low');
0093     case    'high'
0094         [B,A] = butter(Norder, fcut(1)/(fs/2) ,'high');
0095     case    'band'
0096         [B,A] = butter(Norder, fcut(1:2)/(fs/2) );
0097     case    'stop'
0098         [B,A] = butter(Norder, fcut(1:2)/(fs/2) , 'stop');
0099     end
0100     H = freqz(B,A,f,fs);
0101     
0102     if exist('online','var') && online==1
0103         fprintf('online filter with time delay\n')
0104         X = filter(B,A, X );
0105     else
0106         X = filtfilt(B, A, X );
0107     end
0108 end    
0109 
0110 
0111 subplot 221
0112 plot(f,abs(H))
0113 xlabel('Freq[Hz]')
0114 title('Gain of frequency response')
0115 ylim([ 0  1.1])
0116 xlim([f(1) f(end)])
0117 hold on
0118 
0119 subplot 222
0120 plot(f,angle(H))
0121 xlabel('Freq[Hz]')
0122 title('Angle of frequency response')
0123 hold on
0124 xlim([f(1) f(end)])
0125 
0126 subplot 223
0127 plot(t,X)
0128 xlabel('Time[sec]')
0129 title('Impulse response')
0130 hold on
0131 
0132 Xmax = max(X);
0133 plot([0 0],[0 Xmax*1.5],'r--')
0134 
0135 if ~exist('tmax','var') || isempty(tmax), 
0136     ix = find( X > Xmax*0.001 );
0137     if length(ix) < 100,
0138         imax = t0+100;
0139         if imax > length(t)
0140             tmax = max(t);
0141         else
0142             tmax = t(imax);
0143         end
0144     else
0145         tmax = t(ix(end));
0146     end
0147 end
0148 
0149 xlim([-tmax, tmax])
0150 
0151 subplot 224
0152 plot(B)
0153 hold on
0154 plot(A,'--')
0155 xlabel('Tap number')
0156 title('Filter coefficient')
0157 xlim([1 max(length(A),length(B))])
0158 fprintf('# of A and B = %d, %d\n',length(A),length(B))