fir.m

Matlab code for FIR filters

% This script illustrates basic processing and plotting functions
% for examining filter properties.  This script uses the FIR filters
% FIR1 (windowing method) and FIRPM (the Parks-McClellan) filter design
% program to create high-pass filters, FREQZ to compute the filter
% transfer functions, FILTER to compute the impulse response, and ROOTS
% to compute zero locations of the filters.

clear
% General filter specification
n = 44;  %  Filter order
fs = 1;  %  Sampling frequency
impdur = n+10;  %  Impulse response duration in samples

% Window method FIR parameters
wwin = kaiser(n+1,4);  %  Kiaser window with beta parameter 4
fc = .2;  %  Cutoff frequency

% Window method FIR parameters
fa = [0 .18 .22 .5]/(fs/2);  %  Select passband and stopband
am = [0 0 1 1]; %  Select corresponding amplitudes at these frequencies


%  Compute filter coefficient
b1 = fir1(n,fc/(fs/2),'high',wwin);  %  windowing method
bpm = firpm(n,fa,am);  % Parks-McClellan


%  Compute impulse response
imps = [1, zeros(1,round(impdur)-1)];  %  Impulse signal input
y1 = filter(b1,1,imps);   %  Impulse response (windowing method)
ypm = filter(bpm,1,imps);  %  Impulse response (Parks-McClellan)

%  Plot impuse response
figure(1)
stem([0:length(y1)-1],y1/max(abs(y1)), 'gx');  %  Normalize amplitude for plot
hold on
stem([0:length(ypm)-1],ypm/max(abs(ypm)), 'ro'); %  Normalize amplitude for plot
hold off
legend({'windowing', 'Parks-McClellan'})
xlabel('Time sample index')
ylabel('Amplitude')
title('Comparison of impulse response')

%  Plot Transfer Functions
figure(2)
[h1,f1] = freqz(b1,1,1024,fs);
[hpm,fpm] = freqz(bpm,1,1024,fs);
plot(f1,abs(h1),'k--', fpm, abs(hpm), 'b:')
legend({'windowing', 'Parks-McClellan'})
ylabel('Linear Amplitude')
xlabel('Hz')
title('Comparison of TF magnitudes (linear scale)')
figure(3)
plot(f1,20*log10(abs(h1)),'k--', fpm, 20*log10(abs(hpm)), 'b:')
legend({'windowing', 'Parks-McClellan'})
ylabel('Linear Amplitude')
xlabel('Hz')
title('Comparison of TF magnitudes (dB scale)')


% Plot pole zero diagram
zp1 = roots(b1);
zppm = roots(bpm);

figure(4)
%  Generate unit circle for reference
w = [-1:.001:1]*pi;
cpts = exp(j*w);
%  Plot unit circle
plot(real(cpts),imag(cpts),'k-')
hold on
%  Plot zeros as individual markers
%  for window method
plot(real(zp1), imag(zp1), 'og', 'MarkerSize', 9)
hold off
xlabel('Real')
ylabel('Imaginary')
title('Window method filter, Zeros in Z-plane')

figure(5)
%  Plot unit circle
plot(real(cpts),imag(cpts),'k-')
hold on
% Plot zeros as individual markers
%  for Parks-McClellan method
plot(real(zppm), imag(zppm), 'or', 'MarkerSize', 9)
hold off
xlabel('Real')
ylabel('Imaginary')
title('Parks-McClellan method filter, Zeros in Z-plane')