南邮dsp实验报告(8)

Rp = 0.5; % max passband ripple, dB Rs = 40; % min stopband attenuation, dB

% Convert spec to normalized digital frequencies omega_p = 2*pi*Fp/FT;

Wp = 2*Fp/FT; % omega_p/pi omega_s = 2*pi*Fs/FT;

Ws = 2*Fs/FT; % omega_s/pi % Estimate the Filter Order

[N, Wn] = cheb1ord(Wp, Ws, Rp, Rs); % Design the Filter

[num,den] = cheby1(N,Rp,Wn); % Display the transfer function

disp('Numerator Coefficients are ');disp(num); disp('Denominator Coefficients are ');disp(den); % Compute the gain response [g, w] = gain(num,den); % Plot the gain response figure(1);

plot(w/pi,g);grid; axis([0 1 -60 5]);

xlabel('\\omega /\\pi'); ylabel('Gain in dB');

title('Gain Response of a Type 1 Chebyshev Lowpass Filter'); % Find and plot the phase figure(2);

w2 = 0:pi/511:pi;

Hz = freqz(num,den,w2); Phase = unwrap(angle(Hz)); plot(w2/pi,Phase);grid;

xlabel('\\omega /\\pi'); ylabel('Unwrapped Phase (rad)');

title('Unwrapped Phase Response of a Type 1 Chebyshev Lowpass Filter');

% Find and plot the group delay figure(3);

GR = grpdelay(num,den,w2); plot(w2/pi,GR);grid;

xlabel('\\omega /\\pi'); ylabel('Group Delay (sec)');

title('Group Delay of a Type 1 Chebyshev Lowpass Filter');

编写的MATLAB程序如下（计算未畸变的相位响应及群延迟响应）：

% Program Q7_6

Ws = [0.4 0.6]; Wp = [0.2 0.8]; Rp = 0.4; Rs = 50; % Estimate the Filter Order

[N1, Wn1] = buttord(Wp, Ws, Rp, Rs); % Design the Filter

[num,den] = butter(N1,Wn1,'stop');

% Find the frequency response; find and plot unwrapped phase wp = 0:pi/1023:pi; wg = 0:pi/511:pi;

Hz = freqz(num,den,wp); Phase = unwrap(angle(Hz)); figure(1);

plot(wp/pi,Phase); grid;

% axis([0 1 a b]);

xlabel('\\omega /\\pi'); ylabel('Unwrapped Phase (rad)'); title('Unwrapped Phase Response of a Butterworth Bandstop Filter'); % Find and plot the group delay GR = grpdelay(num,den,wg); figure(2);

plot(wg/pi,GR); grid;

%axis([0 1 a b]);

xlabel('\\omega /\\pi'); ylabel('Group Delay (sec)'); title('Group Delay of a Butterworth Bandstop Filter');

0-2-4Unwrapped Phase Response of a Butterworth Bandstop FilterUnwrapped Phase (rad)-6-8-10-12-14-16-18-2000.20.4? /?0.60.81

1098Group Delay of a Butterworth Bandstop FilterGroup Delay (sec)76543200.20.4? /?0.60.81

Q7.9 使用sinc编写一个MATLAB程序，以产生截止频率在?C=0.4?处，长度分别为81,61,41,和21的四个零相位低通滤波器的冲激响应系数，然后计算并画出它们的幅度响应。使用冒号“：”运算符从长度为81的滤波器的冲激响应系数中抽出较短长度滤波器的冲激响应系数。在每一个滤波器的截止频率两边研究频率响应的摆动行为。波纹的数量与滤波器的长度有什么关系？最大波纹的高度与滤波器的长度有什么关系？将怎样修改上述程序以产生一个偶数长度的零相位低通滤波器的冲激响应系数？

（1）长度为81的零相位低通滤波器的冲激响应：

1.41.2Magnitude Response for Length=8110.8|H(ej?)|0.60.40.2000.20.4? /?0.60.81

（2）长度为61的零相位低通滤波器的冲激响应：

1.41.2Magnitude Response for Length=6110.8|H(ej?)|0.60.40.2000.20.4? /?0.60.81

（3）长度为41的零相位低通滤波器的冲激响应：

1.41.2Magnitude Response for Length=4110.8|H(ej?)|0.60.40.2000.20.4? /?0.60.81

（4）长度为21的零相位低通滤波器的冲激响应：

1.41.2Magnitude Response for Length=2110.8|H(ej?)|0.60.40.2000.20.4? /?0.60.81

答：从这些图形可以证明不同情况下的幅值振动都是由于Gibb效应。波纹的数量的减少和滤波器的长度成正比的。最大波纹的高度与滤波器的长度没有关系。 修改程序产生一个偶数长度的零相位低通滤波器的冲激响应系数：

% Program Q7_9

n = -39.5:39.5; % this gives us a length of 80

hn_80 = 0.4 * sinc(0.4*n); % the length-80 impulse response omega = 0:pi/1023:pi; % radian frequency vector W = omega/pi; % Matlab normalized freq vector

Hz_80 = abs(freqz(hn_80,1,omega)); % 1024 samles of |H(e^jw)| figure(1);

plot(W,Hz_80); grid;

xlabel('\\omega /\\pi'); ylabel('|H(e^{j\\omega})|'); title('Magnitude Response for Length=80'); % Reduce length to 60 and repeat hn_60 = hn_80(11:70);

Hz_60 = abs(freqz(hn_60,1,omega)); figure(2);

plot(W,Hz_60); grid;

xlabel('\\omega /\\pi'); ylabel('|H(e^{j\\omega})|'); title('Magnitude Response for Length=60'); % Reduce length to 40 and repeat hn_40 = hn_60(11:50);

Hz_40 = abs(freqz(hn_40,1,omega)); figure(3);

plot(W,Hz_40); grid;

xlabel('\\omega /\\pi'); ylabel('|H(e^{j\\omega})|'); title('Magnitude Response for Length=40');