fb_para

function fb_para
%
% set parameter for filters in filterbank.m as gobal variable
%     to avoid repeating the computing of filter parameters 
%
% Author: Yunbin Deng, work performed at Johns Hopkins Univeristy
%         [e-mail] or [e-mail]
% 
% Note: The parameters need to be customized for specific application
%
% All right reserved, free for academic education and research only.
%


global bhd ahd;     % pre amphesis high pass filter
global bbd abd;     % bandpass filter 
global bld ald ;    % lowpass filter
global bbd_last abd_last; % bandpass
global numChannel;

%Fs    = 20000;       % sampling frequency for tidigit
%Fs     = 8000;        % sampling frequency for yoho
Fs     = 5000;        % sampling frequency for EMG 

% step 1: high pass filter
fh     = 50;
wh     = 2*pi*fh;
[bhd, ahd] = butter(2, wh/Fs, 'high'); % digital high pass filter parameter
                                       % the naming convertion for each letter 
                                           % b:numerator a: denominator
                                           % h:highpass l:lowpass
                                           % d: digital filter parameter

% step 2: N channel Mel/BARK-scale filter bank
fstart  = 50;         % start corner frequency
flinear = 1000;       % linear sacle rang
fstop   = 2500;       % stop corner frequency
numChannel = 24;      % with 16 channel, performance decrease only 0.6%
Q       = 4; % 7 for speaker recognition, 4 for speech recognition

fo   = [linspace(fstart, flinear, 18)   logspace(log10(flinear+100), log10(fstop), 6) ];  
% fo  = [linspace(fstart, flinear, 10) logspace(log10(flinear+150), log10(fstop), 6) ];  
                 % center frequencies approximate mel-scale spaced
%fo = [70 100 150 250 350 450 570   700 840 1000 1170 1370 1600  1850 2150 2500 2900 3400 4000   4800 5800 7000 8500 10500 ];
wo   = 2 * pi * fo;


% step 4: low pass filter cut-off,  100 frequency better than 1000?
fc   = 1000*ones(1,numChannel);   
wc   = 2 * pi * fc;

% nonlinear+lowpass: mimiced by band pass filter at last stage
Q_last = 0.7;
f_last = [linspace(6,20,numChannel) ];
w_last = 2*pi*f_last;

   % to plot the filterbank response for paper
%figure;
%   f = 10:10:100000;
%   w = 2*pi*f;

for j = 1 : numChannel,
   ab  = [1 wo(j)/Q wo(j)^2];    % second order biquad band pass
   bb  = [wo(j)/Q 0];

   % to plot the filterbank response for paper
%   [mag, phase] = bode(bb,ab,w);
%   loglog(f,mag);
%   hold on;

   [bbdj, abdj] = bilinear (bb, ab, Fs);
   bbd = [bbd; bbdj];
   abd = [abd; abdj];

   % step 3: absolute value

   % step 4: lowpass filter envelp extractor
   al   = [1 wc(j)]; % analog low pass biquad filtr parameter
   bl   = wc(j);
   [bldj, aldj] = bilinear(bl,al,Fs); % digital low pass filter parameter
   bld = [bld; bldj]; 
   ald = [ald; aldj];

   % step 5: non-linear compression simulated by bandpass filter
   % y = log(y+0.00001); % prevent log of zero; 99.09% on clean
   ab_last = [1 w_last(j)/Q_last w_last(j)^2];
   bb_last = [w_last(j)/Q_last 0];
   [bbd_lastj, abd_lastj] = bilinear (bb_last, ab_last, Fs);
   bbd_last = [ bbd_last; bbd_lastj];
   abd_last = [ abd_last; abd_lastj];

   % step 6:  decimation
   
end

% step 7: dct

   % to plot the filterbank response for paper
%SetFont('Palatino', 18, 16, 18, 15);
%xlabel('Frequency');
%ylabel('Amplitude');
%hold off;