matlab 8.4,《DSP using MATLAB》Problem 8.42

代码：

%% ------------------------------------------------------------------------

%% Output Info about this m-file

fprintf(‘\n***********************************************************\n‘);

fprintf(‘ Problem 8.42 \n\n‘);

banner();

%% ------------------------------------------------------------------------

% Digital Filter Specifications: Elliptic bandstop

wsbs = [0.40*pi 0.48*pi]; % digital stopband freq in rad

wpbs = [0.25*pi 0.75*pi]; % digital passband freq in rad

Rp = 1.0 % passband ripple in dB

As = 80 % stopband attenuation in dB

Ripple = 10 ^ (-Rp/20) % passband ripple in absolute

Attn = 10 ^ (-As/20) % stopband attenuation in absolute

% Calculation of Elliptic filter parameters:

[N, wn] = ellipord(wpbs/pi, wsbs/pi, Rp, As);

fprintf(‘\n ********* Elliptic Filter Order is = %3.0f \n‘, N)

% Digital Elliptic bandstop Filter Design:

[bbs, abs] = ellip(N, Rp, As, wn, ‘stop‘);

[C, B, A] = dir2cas(bbs, abs)

% Calculation of Frequency Response:

[dbbs, magbs, phabs, grdbs, wwbs] = freqz_m(bbs, abs);

% ---------------------------------------------------------------

% find Actual Passband Ripple and Min Stopband attenuation

% ---------------------------------------------------------------

delta_w = 2*pi/1000;

Rp_bs = -(min(dbbs(1:1:ceil(wpbs(1)/delta_w+1)))); % Actual Passband Ripple

fprintf(‘\nActual Passband Ripple is %.4f dB.\n‘, Rp_bs);

As_bs = -round(max(dbbs(ceil(wsbs(1)/delta_w)+1:1:ceil(wsbs(2)/delta_w)+1))); % Min Stopband attenuation

fprintf(‘\nMin Stopband attenuation is %.4f dB.\n\n‘, As_bs);

%% -----------------------------------------------------------------

%% Plot

%% -----------------------------------------------------------------

figure(‘NumberTitle‘, ‘off‘, ‘Name‘, ‘Problem 8.42 Elliptic bs by ellip function‘)

set(gcf,‘Color‘,‘white‘);

M = 1; % Omega max

subplot(2,2,1); plot(wwbs/pi, magbs); axis([0, M, 0, 1.2]); grid on;

xlabel(‘Digital frequency in \pi units‘); ylabel(‘|H|‘); title(‘Magnitude Response‘);

set(gca, ‘XTickMode‘, ‘manual‘, ‘XTick‘, [0, 0.25, 0.40, 0.48, 0.75, M]);

set(gca, ‘YTickMode‘, ‘manual‘, ‘YTick‘, [0, 0.01, 0.8913, 1]);

subplot(2,2,2); plot(wwbs/pi, dbbs); axis([0, M, -120, 2]); grid on;

xlabel(‘Digital frequency in \pi units‘); ylabel(‘Decibels‘); title(‘Magnitude in dB‘);

set(gca, ‘XTickMode‘, ‘manual‘, ‘XTick‘, [0, 0.25, 0.40, 0.48, 0.75, M]);

set(gca, ‘YTickMode‘, ‘manual‘, ‘YTick‘, [ -80, -40, 0]);

set(gca,‘YTickLabelMode‘,‘manual‘,‘YTickLabel‘,[‘80‘; ‘40‘;‘ 0‘]);

subplot(2,2,3); plot(wwbs/pi, phabs/pi); axis([0, M, -1.1, 1.1]); grid on;

xlabel(‘Digital frequency in \pi nuits‘); ylabel(‘radians in \pi units‘); title(‘Phase Response‘);

set(gca, ‘XTickMode‘, ‘manual‘, ‘XTick‘, [0, 0.25, 0.40, 0.48, 0.75, M]);

set(gca, ‘YTickMode‘, ‘manual‘, ‘YTick‘, [-1:0.5:1]);

subplot(2,2,4); plot(wwbs/pi, grdbs); axis([0, M, 0, 50]); grid on;

xlabel(‘Digital frequency in \pi units‘); ylabel(‘Samples‘); title(‘Group Delay‘);

set(gca, ‘XTickMode‘, ‘manual‘, ‘XTick‘, [0, 0.25, 0.40, 0.48, 0.75, M]);

set(gca, ‘YTickMode‘, ‘manual‘, ‘YTick‘, [0:20:50]);

% ------------------------------------------------------------

% PART 2

% ------------------------------------------------------------

% Discrete time signal

Ts = 1; % sample intevals

n1_start = 0; n1_end = 200;

n1 = [n1_start:n1_end]; % [0:200]

xn1 = sin(0.44*pi*n1); % digital signal

% ----------------------------

% DTFT of xn1

% ----------------------------

M = 500;

[X1, w] = dtft1(xn1, n1, M);

%magX1 = abs(X1);

angX1 = angle(X1); realX1 = real(X1); imagX1 = imag(X1);

magX1 = sqrt(realX1.^2 + imagX1.^2);

%% --------------------------------------------------------------------

%% START X(w)‘s mag ang real imag

%% --------------------------------------------------------------------

figure(‘NumberTitle‘, ‘off‘, ‘Name‘, ‘Problem 8.42 X1 DTFT‘);

set(gcf,‘Color‘,‘white‘);

subplot(2,1,1); plot(w/pi,magX1); grid on; %axis([-1,1,0,1.05]);

title(‘Magnitude Response‘);

xlabel(‘digital frequency in \pi units‘); ylabel(‘Magnitude |H|‘);

set(gca, ‘XTickMode‘, ‘manual‘, ‘XTick‘, [0, 0.2, 0.44, 0.6, 0.8, 1.0, 1.2, 1.4, 1.56, 1.8, 2]);

subplot(2,1,2); plot(w/pi, angX1/pi); grid on; %axis([-1,1,-1.05,1.05]);

title(‘Phase Response‘);

xlabel(‘digital frequency in \pi units‘); ylabel(‘Radians/\pi‘);

figure(‘NumberTitle‘, ‘off‘, ‘Name‘, ‘Problem 8.42 X1 DTFT‘);

set(gcf,‘Color‘,‘white‘);

subplot(2,1,1); plot(w/pi, realX1); grid on;

title(‘Real Part‘);

xlabel(‘digital frequency in \pi units‘); ylabel(‘Real‘);

set(gca, ‘XTickMode‘, ‘manual‘, ‘XTick‘, [0, 0.2, 0.44, 0.6, 0.8, 1.0, 1.2, 1.4, 1.56, 1.8, 2]);

subplot(2,1,2); plot(w/pi, imagX1); grid on;

title(‘Imaginary Part‘);

xlabel(‘digital frequency in \pi units‘); ylabel(‘Imaginary‘);

%% -------------------------------------------------------------------

%% END X‘s mag ang real imag

%% -------------------------------------------------------------------

% ------------------------------------------------------------

% PART 3

% ------------------------------------------------------------

yn1 = filter(bbs, abs, xn1);

n2 = n1;

% ----------------------------

% DTFT of yn1

% ----------------------------

M = 500;

[Y1, w] = dtft1(yn1, n2, M);

%magY1 = abs(Y1);

angY1 = angle(Y1); realY1 = real(Y1); imagY1 = imag(Y1);

magY1 = sqrt(realY1.^2 + imagY1.^2);

%% --------------------------------------------------------------------

%% START Y1(w)‘s mag ang real imag

%% --------------------------------------------------------------------

figure(‘NumberTitle‘, ‘off‘, ‘Name‘, ‘Problem 8.42 Y1 DTFT‘);

set(gcf,‘Color‘,‘white‘);

subplot(2,1,1); plot(w/pi,magY1); grid on; %axis([-1,1,0,1.05]);

title(‘Magnitude Response‘);

xlabel(‘digital frequency in \pi units‘); ylabel(‘Magnitude |H|‘);

subplot(2,1,2); plot(w/pi, angY1/pi); grid on; %axis([-1,1,-1.05,1.05]);

title(‘Phase Response‘);

xlabel(‘digital frequency in \pi units‘); ylabel(‘Radians/\pi‘);

figure(‘NumberTitle‘, ‘off‘, ‘Name‘, ‘Problem 8.42 Y1 DTFT‘);

set(gcf,‘Color‘,‘white‘);

subplot(2,1,1); plot(w/pi, realY1); grid on;

title(‘Real Part‘);

xlabel(‘digital frequency in \pi units‘); ylabel(‘Real‘);

subplot(2,1,2); plot(w/pi, imagY1); grid on;

title(‘Imaginary Part‘);

xlabel(‘digital frequency in \pi units‘); ylabel(‘Imaginary‘);

%% -------------------------------------------------------------------

%% END Y1‘s mag ang real imag

%% -------------------------------------------------------------------

figure(‘NumberTitle‘, ‘off‘, ‘Name‘, ‘Problem 8.42 x(n) and y(n)‘)

set(gcf,‘Color‘,‘white‘);

subplot(2,1,1); stem(n1, xn1);

xlabel(‘n‘); ylabel(‘x(n)‘);

title(‘xn sequence‘); grid on;

subplot(2,1,2); stem(n1, yn1);

xlabel(‘n‘); ylabel(‘y(n)‘);

title(‘yn sequence‘); grid on;

运行结果：

我自己假设通带1dB，阻带衰减80dB。

在此基础上设计指标，绝对单位，

ellip函数(MATLAB工具箱函数)得到Elliptic带阻滤波器，阶数为5，系统函数串联形式系数如下图。

要想得到题目中的10阶的话，阻带衰减估计需要达到160dB左右，觉得没必要那么大。

Elliptic带阻滤波器，幅度谱、相位谱和群延迟响应

输入离散时间信号x(n)的谱如下，可看出，频率分量0.44π

通过带阻滤波器后，得到的输出y(n)的谱，好像变乱了，o(╥﹏╥)o

输入和输出的离散时间序列如下图

原文：https://www.cnblogs.com/ky027wh-sx/p/11808896.html