《DSP using MATLAB》Problem 4.21

Posted 2020-10-27 沧海一粟

        快到龙抬头，居然下雪了，天空飘起了雪花，温度下降了近20°。

代码：

%% ------------------------------------------------------------------------
%%            Output Info about this m-file
fprintf(\'\\n***********************************************************\\n\');
fprintf(\'        <DSP using MATLAB> Problem 4.21 \\n\\n\');

banner();
%% ------------------------------------------------------------------------


% ----------------------------------------------------
%               1        H1(z)
% ----------------------------------------------------

b = [3/4, 5/4, 1, 1, 5/4, 3/4]; a = [1];                             %  

[R, p, C] = residuez(b,a)

Mp = (abs(p))\'
Ap = (angle(p))\'/pi

%% ------------------------------------------------------
%%   START a    determine H(z) and sketch    
%% ------------------------------------------------------
figure(\'NumberTitle\', \'off\', \'Name\', \'P4.21 H(z) its pole-zero plot\')
set(gcf,\'Color\',\'white\'); 
zplane(b,a);
title(\'pole-zero plot\'); grid on;

%% ----------------------------------------------
%%    END
%% ----------------------------------------------

% ------------------------------------
%                  h(n)  
% ------------------------------------

[delta, n] = impseq(0, 0, 19); 
h_check = filter(b, a, delta);                                             % check sequence


%% --------------------------------------------------------------
%%    START    b   |H|   <H
%%    3rd form of freqz
%% --------------------------------------------------------------
w = [-500:1:500]*2*pi/500;     H = freqz(b,a,w); 
%[H,w] = freqz(b,a,200,\'whole\');                 % 3rd form of freqz

magH  = abs(H);  angH  = angle(H);  realH  = real(H);  imagH  = imag(H);

%% ================================================
%%              START H\'s  mag ang real imag
%% ================================================
figure(\'NumberTitle\', \'off\', \'Name\', \'P4.21  DTFT and Real Imaginary Part \');
set(gcf,\'Color\',\'white\'); 
subplot(2,2,1); plot(w/pi,magH); grid on;  %axis([0,1,0,1.5]); 
title(\'Magnitude Response\');
xlabel(\'frequency in \\pi units\'); ylabel(\'Magnitude  |H|\'); 
subplot(2,2,3); plot(w/pi, angH/pi); grid on; % axis([-1,1,-1,1]);
title(\'Phase Response\');
xlabel(\'frequency in \\pi units\'); ylabel(\'Radians/\\pi\');

subplot(\'2,2,2\'); plot(w/pi, realH); grid on;
title(\'Real Part\');
xlabel(\'frequency in \\pi units\'); ylabel(\'Real\');
subplot(\'2,2,4\'); plot(w/pi, imagH); grid on;
title(\'Imaginary Part\');
xlabel(\'frequency in \\pi units\'); ylabel(\'Imaginary\');
%% ==================================================
%%             END H\'s  mag ang real imag
%% ==================================================



% --------------------------------------------------------------
%        x(n) through the filter, we get output y(n)
% --------------------------------------------------------------
N = 200;
nx = [0:1:N-1];
 x = sin(pi*nx/2) + 5 * cos(pi*nx);

[y, ny] = conv_m(h_check, n, x, nx);

figure(\'NumberTitle\', \'off\', \'Name\', \'P4.21  Input & h(n) Sequence\');
set(gcf,\'Color\',\'white\'); 
subplot(3,1,1); stem(nx, x); grid on;  %axis([0,1,0,1.5]); 
title(\'x(n)\');
xlabel(\'n\'); ylabel(\'x\'); 
subplot(3,1,2); stem(n, h_check); grid on;  %axis([0,1,0,1.5]); 
title(\'h(n)\');
xlabel(\'n\'); ylabel(\'h\'); 
subplot(3,1,3); stem(ny, y); grid on;  %axis([0,1,0,1.5]); 
title(\'y(n)\');
xlabel(\'n\'); ylabel(\'y\');

figure(\'NumberTitle\', \'off\', \'Name\', \'P4.21  Output Sequence\');
set(gcf,\'Color\',\'white\'); 
subplot(1,1,1); stem(ny, y); grid on;  %axis([0,1,0,1.5]); 
title(\'y(n)\');
xlabel(\'n\'); ylabel(\'y\'); 


% ----------------------------------------
%                yss Response
% ----------------------------------------
ax = conv([1,0,1], [1,2,1])
bx = conv([0,1], [1,2,1]) + conv([5,5], [1,0,1])

by = conv(bx, b) 
ay = ax

zeros = roots(by)

[R, p, C] = residuez(by, ay)

Mp_Y = (abs(p))\'
Ap_Y = (angle(p))\'/pi

%% ------------------------------------------------------
%%   START a    determine Y(z) and sketch    
%% ------------------------------------------------------
figure(\'NumberTitle\', \'off\', \'Name\', \'P4.21 Y(z) its pole-zero plot\')
set(gcf,\'Color\',\'white\'); 
zplane(by, ay);
title(\'pole-zero plot\'); grid on;


% ------------------------------------
%                  y(n)  
% ------------------------------------
LENGH = 100;
[delta, n] = impseq(0, 0, LENGH-1); 
y_check = filter(by, ay, delta);                                             % check sequence

y_answer0 = 4.75*delta;

[delta_1, n1] = sigshift(delta, n, 1);
y_answer1 = 2.25*delta_1;

[delta_2, n2] = sigshift(delta, n, 2);
y_answer2 = 2.75*delta_2;

[delta_3, n3] = sigshift(delta, n, 3);
y_answer3 = 3.75*delta_3; 

[delta_4, n4] = sigshift(delta, n, 4);
y_answer4 = 4.50*delta_4;

y_answer5 = (2*(-0.5)*cos(pi*n/2) + 2*0.5*sin(pi*n/2) ).*stepseq(0,0,LENGH-1);


[y01, n01] = sigadd(y_answer0, n,  y_answer1, n1);
[y02, n02] = sigadd(y_answer2, n2, y_answer3, n3);

[y03, n03] = sigadd(y01, n01, y02, n02);
[y04, n04] = sigadd(y03, n03, y_answer4, n4);

[y_answer, n_answer] = sigadd(y04, n04, y_answer5, n);


figure(\'NumberTitle\', \'off\', \'Name\', \'P4.21  Yss and Y \');
set(gcf,\'Color\',\'white\'); 
subplot(2,1,1); stem(n, y_answer5); grid on;  %axis([0,1,0,1.5]); 
title(\'Steady-State Response\');
xlabel(\'n\'); ylabel(\'Yss\'); 
subplot(2,1,2); stem(n, y_check); grid on; % axis([-1,1,-1,1]);
title(\'Total Response\');
xlabel(\'n\'); ylabel(\'Y\');

　　运行结果：

        系统函数H(z)的系数：

        系统的DTFT，注意当ω=π/2和π时的振幅谱、相位谱的值。

        当有输入时，输出的Y(z)进行部分分式展开，留数及对应的极点如下：

        单位圆上z=-1处，极点和零点相互抵消，稳态响应只和正负j有关。