Examples of FIR filter design using Parks-McClellan algorithm
Use the Matlab function 'firpm'
Contents
Start
clear
Design FIR filter using Parks-McClellan algorithm
Low-pass filter design
fp = 0.1; fs = 0.15; [h, del] = firpm(30, [0 fp fs .5]*2, [1 1 0 0]); length(h) del figure(1) stem(0:30, h, 'filled') title('Impulse response')
ans = 31 del = 0.0241
Compute amplitude response A
L = 512; [A, om] = firamp(h, 1, L); f = om/(2*pi); figure(1) clf plot(f, A) xlim([0 0.5]) title('Amplitude response') xlabel('Frequency (normalized)')
Plot with delta
figure(1) clf plot(f, A, [0 fp], (1-del)*[1 1],'r', [0 fp], (1+del)*[1 1], 'r', [fs 0.5], -del*[1 1], 'r', [fs 0.5], del*[1 1], 'r') xlim([0 0.5]) title('Amplitude response') xlabel('Frequency (normalized)')
Use firpm with weighting
The weight function allows one to put more weight in one band than in the other band.
fp = 0.1; fs = 0.15; W = [1 10]; [h, del] = firpm(30, [0 fp fs .5]*2, [1 1 0 0], W); del
del = 0.0845
Compute amplitude response A
L = 512; [A, om] = firamp(h, 1, L); f = om/(2*pi); % Oops. Pass-band and stop-band ripples are different in size figure(1) clf plot(f, A, [0 fp], (1-del)*[1 1],'r', [0 fp], (1+del)*[1 1], 'r', [fs 0.5], -del*[1 1], 'r', [fs 0.5], del*[1 1], 'r') xlim([0 0.5]) title('Amplitude response') xlabel('Frequency (normalized)')
Plot with delta
Correct plot . . . Note: The stop-band ripple is one-tength the pass-band ripple
figure(1) clf plot(f, A, [0 fp], (1-del)*[1 1],'r', [0 fp], (1+del)*[1 1], 'r', [fs 0.5], -del/W(2)*[1 1], 'r', [fs 0.5], del/W(2)*[1 1], 'r') xlim([0 0.5]) title('Amplitude response') xlabel('Frequency (normalized)')
High-pass filter
fs = 0.1; fp = 0.15; [h, del] = firpm(30, [0 fs fp .5]*2, [0 0 1 1]); L = 512; [A, om] = firamp(h, 1, L); figure(1) clf plot(f, A) xlim([0 0.5]) title('Amplitude response') xlabel('Frequency (normalized)')
High-pass filter (Even-length)
Why does this produce an error?
fs = 0.1; fp = 0.15; [h, del] = firpm(31, [0 fs fp .5]*2, [0 0 1 1]); % Because 31 -> length(h) = 32 -> h is Type II -> Hf(pi) = 0 % which is inconsistent with the specification!
Warning: Odd order symmetric FIR filters must have a gain of zero at the Nyquist frequency. The order is being increased by one. Alternatively, you can pass a trailing 'h' argument, as in firpm(N,F,A,W,'h'), to design a type 4 linear phase filter.
Band-pass filter
fs1 = 0.05; fp1 = 0.1; fp2 = 0.2; fs2 = 0.25; [h, del] = firpm(50, [0 fs1 fp1 fp2 fs2 .5]*2, [0 0 1 1 0 0]); L = 512; [A, om] = firamp(h, 1, L); figure(1) clf plot(f, A) xlim([0 0.5]) title('Amplitude response') xlabel('Frequency (normalized)')
Band-stop filter
fp1 = 0.05; fs1 = 0.1; fs2 = 0.2; fp2 = 0.25; [h, del] = firpm(50, [0 fp1 fs1 fs2 fp2 .5]*2, [1 1 0 0 1 1]); L = 512; [A, om] = firamp(h, 1, L); figure(1) clf plot(f, A) xlim([0 0.5]) title('Amplitude response') xlabel('Frequency (normalized)')
Custom filter
f1 = 0.15; f2 = 0.2; [h, del] = firpm(42, [0 f1 f2 .5]*2, [0 1 1 1]); L = 512; A = firamp(h, 1, L); figure(1) clf plot(f, A) xlim([0 0.5]) title('Amplitude response') xlabel('Frequency (normalized)')
Type III and Type IV FIR filters
To obtain an anti-symmetric impulse response, use 'hilbert' in firpm. In this case, we must have Hf(0) = 0.
fs = 0.1; fp = 0.15; [h, del] = firpm(31, [0 fs fp .5]*2, [0 0 1 1], 'hilbert'); L = 512; A = firamp(h, 3, L); figure(1) clf subplot(2, 1, 1) stem(h, 'filled') subplot(2, 1, 2) plot(f, A) xlim([0 0.5]) title('Amplitude response') xlabel('Frequency (normalized)')