F = [0 Wp Ws 1]; A = [1 1 0 0]; h = firpm(N,F,A);
% Show the Numerator Coefficients
disp('Numerator Coefficients are ' );disp(h); % Compute and plot the gain response [g, w] = gain(h,[1]); figure(1);
plot(w/pi,g);grid;
xlabel('\\omega /\\pi' ); ylabel('Gain in dB' ); title('Gain Response' );
% Compute the frequency response w2 = 0:pi/511:pi; Hz = freqz(h,[1],w2);
% TEST: did we meet the spec? MagH = abs(Hz);
T1 = 1.005*ones(1,length(w2)); T2 = 0.995*ones(1,length(w2)); T3 = 0.005*ones(1,length(w2)); figure(4);
plot(w2/pi,MagH,w2/pi,T1,w2/pi,T2,w2/pi,T3);grid; % Find and plot the phase figure(2);
Phase = angle(Hz);
plot(w2/pi,Phase);grid;
xlabel('\\omega /\\pi' ); ylabel('Phase (rad)' ); title('Phase Response' ); figure(3);
UPhase = unwrap(Phase); plot(w2/pi,UPhase);grid;
xlabel('\\omega /\\pi' ); ylabel('Unwrapped Phase (rad)' ); title('Unwrapped Phase Response' ); 低通滤波器系数：
0.0028 -0.0022 -0.0046 -0.0006 0.0053 0.0019 -0.0073 -0.0058 0.0079 0.0106 -0.0069 -0.0170 0.0032 0.0243 0.0045 -0.0319 -0.0182 0.0390 0.0422 -0.0448 -0.0924 0.0486 0.3136 0.4501 0.3136 0.0486 -0.0924 -0.0448 0.0422 0.0390 -0.0182 -0.0319 0.0045 0.0243 0.0032 -0.0170 -0.0069 0.0106 0.0079 -0.0058 -0.0073 0.0019 0.0053 -0.0006 -0.0046 -0.0022 0.0028 增益和相位响
从图中可以看出此时的滤波器不满足指标。欲满足指标，应调节N=47. Q7.23 用凯泽窗设计一个有限冲激响应低通滤波器。 程序：
% Program Q7_23
% Use Kaiser window to design a linear phase Lowpass
% FIR Digital Filter meeting the design specification given % in Q7.23.
%
% It is clear from the statement of the question that Mitra % wants us to use (7.36) and (7.37) for this problem. That
% isn't the greatest thing to try because kaiserord already does % exactly what we need.... but that's Q7_24! So here goes! % - Print out the numerator coefficients % for the transfer function.
% - Compute and plot the gain function.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% clear;
% Design spec as given in Q7.23. Wp = 0.31; Ws = 0.41;
Wn = Wp + (Ws-Wp)/2; As = 50;
Ds = 10^(-As/20);
Dp = Ds; %Kaiser window design has equal ripple in % passband and stopband.
% estimate order using (7.37) if As > 21
N = ceil((As-7.95)*2/(14.36*(abs(Wp-Ws)))+1) else
N = ceil(0.9222*2/abs(Wp-Ws)+1) end
% Use (7.36) to get Beta if As > 50
BTA = 0.1102*(As-8.7); elseif As >= 21
BTA = 0.5842*(As-21)^0.4+0.07886*(As-21); else BTA = 0; end
Win = kaiser(N+1,BTA); h = fir1(N,Wn,Win);
% Show the Numerator Coefficients
disp('Numerator Coefficients are ' );disp(h); % Compute and plot the gain response [g, w] = gain(h,[1]); figure(1);
plot(w/pi,g);grid; axis([0 1 -80 5]);
xlabel('\\omega /\\pi' ); ylabel('Gain in dB' ); title('Gain Response' );
% Compute the frequency response
w2 = 0:pi/511:pi; Hz = freqz(h,[1],w2); % Find and plot the phase figure(2);
Phase = angle(Hz);
plot(w2/pi,Phase);grid;
xlabel('\\omega /\\pi' ); ylabel('Phase (rad)' ); title('Phase Response' ); figure(3);
UPhase = unwrap(Phase); plot(w2/pi,UPhase);grid;
xlabel('\\omega /\\pi' ); ylabel('Unwrapped Phase (rad)' ); title('Unwrapped Phase Response' ); 低通滤波器系数：
0.0003 0.0008 0.0003 -0.0011 -0.0017 0.0000 0.0026 0.0027 -0.0010 -0.0049 -0.0035 0.0033 0.0080 0.0034 -0.0074 -0.0119 -0.0018 0.0140 0.0161 -0.0027 -0.0241 -0.0201 0.0127 0.0406 0.0236 -0.0354 -0.0754 -0.0258 0.1214 0.2871 0.3597 0.2871 0.1214 -0.0258 -0.0754 -0.0354 0.0236 0.0406 0.0127 -0.0201 -0.0241 -0.0027 0.0161 0.0140 -0.0018 -0.0119 -0.0074 0.0034 0.0080 0.0033 -0.0035 -0.0049 -0.0010 0.0027 0.0026 0.0000 -0.0017 -0.0011 0.0003 0.0008 0.0003
增益和相位响应如下：
从图中可以看出设计的滤波器满足要求。N=60.
Q7.24 用函数kaiserord和firl重做习题Q7.23
程序：
% Use Kaiser window to design a linear phase Lowpass
% FIR Digital Filter meeting the design specification given % in Q7.23. Use kaiserord and fir1.
% - Print out the numerator coefficients % for the transfer function.
% - Compute and plot the gain function.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% clear;
% Design spec as given in Q7.23. Wp = 0.31; Ws = 0.41; As = 50;
Ds = 10^(-As/20); % Design the Filter F = [Wp Ws]; A = [1 0]; DEV = [Ds Ds];
[N,Wn,BTA,Ftype] = kaiserord(F,A,DEV); Win = kaiser(N+1,BTA); h = fir1(N,Wn,Ftype,Win);
% Show the Numerator Coefficients
disp('Numerator Coefficients are ' );disp(h); % Compute and plot the gain response [g, w] = gain(h,[1]); figure(1);
plot(w/pi,g);grid; axis([0 1 -80 5]);
xlabel('\\omega /\\pi' ); ylabel('Gain in dB' ); title('Gain Response' );
% Compute the frequency response w2 = 0:pi/511:pi; Hz = freqz(h,[1],w2); % Find and plot the phase figure(2);
Phase = angle(Hz);
plot(w2/pi,Phase);grid;
xlabel('\\omega /\\pi' ); ylabel('Phase (rad)' ); title('Phase Response' ); figure(3);
UPhase = unwrap(Phase); plot(w2/pi,UPhase);grid;
xlabel('\\omega /\\pi' ); ylabel('Unwrapped Phase (rad)' );