MATLAB离散系统的频响曲线与零极点图绘制

%画出因果系统y(n)=0.9y(n-1)+x(n)的零极点图和频率响应曲线,并求单位抽样响应.
N=50;
a=[1 0];
b=[1 -0.9];
figure(1)
subplot(2,2,1);
zplane(a,b);                 %根据系统函数返回零极点图
xlabel('Real part');
ylabel('Imaginary part');
grid on; 
              
x=[1 zeros(1,N-1)];
n=0:N-1;
subplot(2,2,2);
stem(n,filter(a,b,x));       %根据系统函数返回单位冲激响应(离散的)
xlabel('n');
ylabel('h(n)');
grid on;
                     
[H,w]=freqz(a,b,N);          %根据系统函数返回0-pi间等间隔的N个频率w相应的频率响应H
subplot(2,2,3);
plot(w,abs(H));
xlabel('w');
ylabel('|H(exp(jw))|');
grid on;
subplot(2,2,4);
plot(w,angle(H));
xlabel('w');
ylabel('arg(H(exp(jw)))');
grid on;

% 一个三阶低通滤波器由下面差分方程描述:
% y(n)=0.0181x(n)+0.0543x(n-1)+0.0543x(n-2)+0.0181x(n-3)+1.76y(n-1)-1.1829y(n-2)+0.2781y(n-3)
% 画出这个滤波器的幅度和相位响应.

%不用MATLAB提供的函数
w=-2*pi:pi/100:2*pi;
H1=(0.0181+0.0543*exp(-1j*w)+0.0543*exp(-1j*w*2)+0.0181*exp(-1j*w*3))./(1-1.76*exp(-1j*w)+1.1829*exp(-1j*w*2)-0.2781*exp(-1j*w*3));
figure(1)
subplot(2,2,1);
plot(w,abs(H1),'.');
xlabel('w');
ylabel('|H(exp(jw))|');
grid on;

subplot(2,2,2);
plot(w,angle(H1),'.');
xlabel('w');
ylabel('ArgH(exp(jw))');
grid on;


% 用MATLAB提供函数
% N=500;
a=[0.0181 0.0543 0.0543 0.0181];
b=[1 -1.76 1.1829 -0.2781];
[H2,w]=freqz(a,b,-2*pi:pi/100:2*pi);      %自定义区间   
% [H2,w]=freqz(a,b,N,'whole');              %0-2*pi区间
subplot(2,2,3);
plot(w,abs(H2));
xlabel('w');
ylabel('|H(exp(jw))|');
grid on;
subplot(2,2,4);
plot(w,angle(H2));
xlabel('w');
ylabel('arg(H(exp(jw)))');
grid on;

% 已知离散系统函数H(z)有一个零点在z=-2,两个极点在z=0.5*exp(j*2*pi/3)及其共轭位置.若其直流增益为1.
% 求:
% (1)H(z)的系数,单位冲激响应,并画出系统零极点图和频率响应曲线.
% (2)将零点移至镜像位置,重复(1)并与之比较,看看有哪些区别.

% (1)
B=poly([0.5*exp(1j*2*pi/3) 0.5*exp(-1j*2*pi/3)]); %poly函数由根求出多项式的系数.
k=1.75/3;
a=[0 k 2*k];
b=[B(1) B(2) B(3)];

figure(1)
subplot(3,1,1)
zplane(a,b);                 %根据系统函数返回零极点图
xlabel('Real part');
ylabel('Imaginary part');
grid on; 


subplot(3,1,2)
[H,w]=freqz(a,b,-2*pi:pi/100:2*pi);      %自定义区间   
plot(w,abs(H));
xlabel('w');
ylabel('|H(exp(jw))|');
grid on;

subplot(3,1,3);
plot(w,angle(H));
xlabel('w');
ylabel('arg(H(exp(jw)))');
grid on;

%(2)
B=poly([0.5*exp(1j*2*pi/3) 0.5*exp(-1j*2*pi/3)]); %poly函数由根求出多项式的系数.
k=1.75*2/3;
a=[0 k 1/2*k];
b=[B(1) B(2) B(3)];

figure(2)
subplot(3,1,1)
zplane(a,b);                 %根据系统函数返回零极点图
xlabel('Real part');
ylabel('Imaginary part');
grid on; 


subplot(3,1,2)
[H,w]=freqz(a,b,-2*pi:pi/100:2*pi);      %自定义区间   
plot(w,abs(H));
xlabel('w');
ylabel('|H(exp(jw))|');
grid on;

subplot(3,1,3);
plot(w,angle(H));
xlabel('w');
ylabel('arg(H(exp(jw)))');
grid on;