Home > functions > device > filter_func > vb_plot_filter_response.m

vb_plot_filter_response

PURPOSE ^

plot filter response in frequency and Impulse response

SYNOPSIS ^

function vb_plot_filter_response(ftype,fcut,fs,Norder,fmax,tmax)

DESCRIPTION ^

 plot filter response in frequency and Impulse response
  vb_plot_filter_response(ftype,fcut,fs,Norder,fmax,tmax)
 --- Input (required)
  ftype ; filter type
    'low'
    'high'
    'band'
    'high_exp' : Exponential highpass filter

  fs   : Sampling frequency [Hz]
  fcut : Cutoff frequency [Hz]
       = [f1 f2] Bandpass
       = [f1]    Lowpass/Highpass
  Norder = 0: FIR filter (eegfilt)
  Norder > 0: filter order of Butterworth filter
 --- Input (optional)
 fmax : Max frequency in frequency plot [Hz]
 tmax : Max time in Impulse plot [sec]

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

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

0001 function    vb_plot_filter_response(ftype,fcut,fs,Norder,fmax,tmax)
0002 % plot filter response in frequency and Impulse response
0003 %  vb_plot_filter_response(ftype,fcut,fs,Norder,fmax,tmax)
0004 % --- Input (required)
0005 %  ftype ; filter type
0006 %    'low'
0007 %    'high'
0008 %    'band'
0009 %    'high_exp' : Exponential highpass filter
0010 %
0011 %  fs   : Sampling frequency [Hz]
0012 %  fcut : Cutoff frequency [Hz]
0013 %       = [f1 f2] Bandpass
0014 %       = [f1]    Lowpass/Highpass
0015 %  Norder = 0: FIR filter (eegfilt)
0016 %  Norder > 0: filter order of Butterworth filter
0017 % --- Input (optional)
0018 % fmax : Max frequency in frequency plot [Hz]
0019 % tmax : Max time in Impulse plot [sec]
0020 %
0021 % Copyright (C) 2011, ATR All Rights Reserved.
0022 % License : New BSD License(see VBMEG_LICENSE.txt)
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     % Exponential online highpass filter
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))

Generated on Tue 27-Aug-2013 11:46:04 by m2html © 2005