% FIR filter design examples
% demonstrating the use of the Kaiser
% window.
%
% R Kakarala
% U C Berkeley Extension

% Example 1 parameters
Fpe = 0.2;  Fse = 0.22;   % passband and stopband edges
eps_p = 0.01; As = 40;   % passband variation and stopband
                          % attenuation in dB
                          
eps_s = 10^(As/-20);      % compute equivalent stopband limit
if (eps_s > eps_p)        % do we need to increase attenuation
                          % to get passband criterion? 
   eps_s = eps_p;
   As = -20*log10(eps_s);
   
end;

DeltaF = Fse - Fpe;       % compute number of coefficients
M=(As-8)/(2.285*(2*pi*DeltaF));
M=ceil(M);                % increase to next largest integer

w = kaiserwindow(M,As);   % compute Kaiser window
Fc = (Fse+Fpe)/2;         % ideal filter cutoff
n = 0:(M-1);
hi = 2*Fc*sinc(2*Fc*(n - (M-1)/2)); % ideal filter

hw = w.*hi;               % windowed sinc;

% frequency response;
[Hw,F] = iirfreqresp(hw);  
figure,plot(F,abs(Hw)),title('example 1');

% check if criteria were met

Hstop = Hw(abs(F)>=Fse);  % response in stopband
fprintf(1,'Example 1 results...\n');
fprintf(1,'Number of coefficients used = %4d\n',M);
fprintf(1,'Min Stopband Attenuation (dB)=%8.5f\n',min(-20*log10(abs(Hstop)+eps)) );
fprintf(1,'Target was=%8.5f (dB)\n',As);
Hpass = Hw(abs(F)<=Fpe);  % response in stopband
fprintf(1,'Max Passband Variation =%8.5f\n',max(abs(abs(Hpass)-1)) );
fprintf(1,'Target was=%8.5f (dB)\n',eps_p);

% Example 2 parameters
Fpe=100/400; Fse=120/400; % passband and stopband edges
eps_p = 0.01; As = 60;   % passband variation and stopband
                          % attenuation in dB
                          
eps_s = 10^(As/-20);      % compute equivalent stopband limit
if (eps_s > eps_p)        % do we need to increase attenuation
                          % to get passband criterion? 
   eps_s = eps_p;
   As = -20*log10(eps_s);
   
end;

DeltaF = Fse - Fpe;       % compute number of coefficients
M=(As-8)/(2.285*(2*pi*DeltaF));

M=ceil(M);                % increase to next largest integer
M=M+1;                    % take next number (even) to get null at Nyquist

w = kaiserwindow(M,As);   % compute Kaiser window
Fc = (Fse+Fpe)/2;         % ideal filter cutoff
n = 0:(M-1);
hi = 2*Fc*sinc(2*Fc*(n - (M-1)/2)); % ideal filter

hw = w.*hi;               % windowed sinc;

% frequency response;
[Hw,F] = iirfreqresp(hw);  
figure,plot(F,abs(Hw)),title('example 2');

% check if criteria were met

Hstop = Hw(abs(F)>=Fse);  % response in stopband
fprintf(1,'Example 2 results...\n');
fprintf(1,'Number of coefficients used = %4d\n',M);
fprintf(1,'Min Stopband Attenuation (dB)=%8.5f\n',min(-20*log10(abs(Hstop)+eps)) );
fprintf(1,'Target was=%8.5f (dB)\n',As);
Hpass = Hw(abs(F)<=Fpe);  % response in stopband
fprintf(1,'Max Passband Variation =%8.5f\n',max(abs(abs(Hpass)-1)) );
fprintf(1,'Target was=%8.5f (dB)\n',eps_p);