信号与系统_钱玲 课后MATLAB学习代码

t=-5:0.01:5;
y1=heaviside(t);
subplot(1,2,1);plot(t,y1);axis([-5,5,-0.5,1.5]);
xlabel('t');ylabel('u(t)');title('单位阶跃信号');

t1=-5*pi:0.01:5*pi;y2=sinc(t1);
subplot(1,2,2);plot(t1,y2);axis([-5*pi,5*pi,-0.5,1.5]);
xlabel('t');ylabel('sinc(t)');title('抽样信号');
line([-5*pi,5*pi],[0,0]);line([0,0],[-0.5,1.5]);

syms t;
%f=sym('heaviside(t+2)-heaviside(t-3)');
%f=sym('cos(2*pi*t+pi/3)');
%f=sym('(2*exp(-t)-exp(-2*t))*heaviside(t-1)');
%f=sym('(0.2*t-2)*heaviside(t)');
%f=sym('sinc(t)');
%f1=subs(f,t,2-2*t);
fplot(f1);axis([-4 4 -2 2]);

clear,clc
syms t;
f=sym('heaviside(t+1)-t*heaviside(t)+(t-1)*heaviside(t-1)');
f1=2*subs(f,t,0.5*t);
f2=subs(f,t,-2*t+3);
f3=diff(f);
f4=int(f);
subplot(2,3,1);fplot(f,[-3,3]);set(gca, 'YLim', [-1,2]);%axis(-2,2,-3,3);
subplot(2,3,2);fplot(f1,[-3,3]);set(gca, 'YLim', [-1,2]);%axis(-2,2,-3,3);
subplot(2,3,3);fplot(f2,[0,3]);set(gca, 'YLim', [-1,2]);%axis(-2,2,0,3);
subplot(2,3,4);fplot(f3,[-2,2]);
line([-1,-1],[0,1]);set(gca, 'YLim', [-2,2]);%axis(-2,2,-1.5,1.5);
subplot(2,3,5);fplot(f4,[-2,3]);set(gca, 'YLim', [-2,2]);

%%
a=[1,2,2];
b=[1 3];
sys=tf(b,a);
t=0:0.01:6;
f=exp(-t);
lsim(sys,f,t);
%gtext('系统激励');gtext('系统响应');
%%
a=[1 2 2 1];b=[3 2];sys=tf(b,a);
subplot(1,2,1);impulse(b,a,8);
subplot(1,2,2);step(b,a,10);

k1=0:0.01:5;k2=-1:0.01:3;p=0.01;%p为采样率
f1=heaviside(k1-1)-heaviside(k1-4);
f2=0.5*k2.*(heaviside(k2)-heaviside(k2-2));
f=conv(f1,f2);f=f*p;%conv离散序列卷积和，是对函数整体序列进行卷积，记得乘上采样率p
k0=k1(1)+k2(1);
k3=length(k1)+length(k2)-2;
k=k0:p:k0+k3*p;
subplot(1,3,1);plot(k1,f1);axis([0 5 -1 2]);
subplot(1,3,2);plot(k2,f2);axis([-1 3 -1 2]);
subplot(1,3,3);plot(k,f);axis([0 7 -1 2]);

%%
syms t;
a=[1 1 3];b=[2 1];
sys=tf(b,a);
t=0:0.01:5;
u=exp(-0.5*t).*heaviside(t);
lsim(sys,u,t);
%%
a=[1 3 -1];b=[1 -3];
sys=tf([1 -3],[1 3 -1]);
subplot(1,2,1);impulse(b,a,8);
subplot(1,2,2);step(b,a,10);
%%
k1=-1:0.01:6;k2=-1:0.01:6;p=0.01;
f1=heaviside(k1)-2*heaviside(k1-3);f2=k2.*(heaviside(k2-1)-heaviside(k2-4));
f=conv(f1,f2);f=f*p;
k=k1(1)+k2(1):p:k1(1)+k2(1)+(length(f)-1)*p;
subplot(1,3,1);plot(k1,f1);axis([-1 6 -2 2]);
subplot(1,3,2);plot(k2,f2);axis([-1 6 -1 5]);
subplot(1,3,3);plot(k,f);%axis([-2 10 -1 2]);

%%
n=0:15;
x1=0.6.^n;x2=(-0.6).^n;x3=1.2.^n;x4=(-1.2).^n;
subplot(141);stem(n,x1,'filled');
subplot(142);stem(n,x2,'filled');
subplot(143);stem(n,x3,'filled');
subplot(144);stem(n,x4,'filled');
%%
n=-4:8;
x1=jyxl(n+2)-jyxl(n-4);x2=jyxl(n)-jyxl(n-5);
n1=0:20;
x3=sin(n1*pi/5).*jyxl(n1);x4=sin(5*n1);
subplot(141);stem(n,x1,'filled');
subplot(142);stem(n,x2,'filled');
subplot(143);stem(n1,x3,'filled');
subplot(144);stem(n1,x4,'filled');
%%
clear;clc
n1=0:15;x=n1;
n2=0:11;h=0.8.^n2;
y=conv(n1,n2);n=0:length(y)-1;
subplot(141);stem(n1,x,'filled');
subplot(142);stem(n2,h,'filled');
subplot(143);stem(n,y,'filled');
p=get(gca,'position');p(3)=2.4*p(3);set(gca,'position',p);

%%
clear;clc
n1=0:15;x=n1;
n2=0:11;h=(0.8).^n2;
y=conv(n1,h);n=0:length(y)-1;
subplot(141);stem(n1,x,'filled');
subplot(142);stem(n2,h,'filled');
subplot(143);stem(n,y,'filled');
p=get(gca,'position');p(3)=2.4*p(3);set(gca,'position',p);
%%
a=[2 -3 2];b=[1 1];n=0:50;
x=4*sin(n*pi/5);y=filter(b,a,x);
subplot(211);stem(n,x,'filled');
subplot(212);stem(n,y,'filled');

%%
n=-1:7;
x1=jyxl(n)-jyxl(n-5);x2=jyxl(n-2)-jyxl(n-6);
n2=0:20;
x3=sin(n2*pi/3);x4=sin(3*n2);
% subplot(141);stem(n,x1,'filled');
% subplot(142);stem(n,x2,'filled');
% subplot(143);stem(n2,x3,'filled');
% subplot(144);stem(n2,x4,'filled');
%%
clear;clc
zx=[1 -0.5];zy=[0 0];
px=[0.5 0.5];py=[0.7 -0.7];
plot(zx,zy,'o',px,py,'x');
%%
b=[1 -1.5 0.5];a=[1 -1 0.74];
n=0:50;
y=filter(b,a,jyxl(n));
subplot(211);stem(n,jyxl(n),'filled');
subplot(212);stem(n,y,'filled');
%%
n1=0:8;x=(0.9).^n1;
n2=0:11;h=(n2+2);
y=conv(x,h);n=0:length(y)-1;
subplot(141);stem(n1,x,'filled');
subplot(142);stem(n2,h,'filled');
subplot(143);stem(n,y,'filled');
p=get(gca,'position');p(3)=2.4*p(3);set(gca,'position',p);

%%
n=-10:10;w1=0.4*pi;
n1=-10:-1;ft1=sin(0.2*pi*n1)./(pi*n1);%直接使用sinc函数有问题，可能问题在于sinc使用'/'未使用'./'
n2=1:10;ft2=sin(0.2*pi*n2)./(pi*n2);
ft=[ft1 0.2 ft2];
fn=abs(ft);phase=angle(ft);
% subplot(141);stem(n,fn);%指数形式的幅频和相频特性
% subplot(142);stem(n,phase);
syms t;s1=0.2;s2=0.2;
for k1=1:5
    s1=s1+2*sin(k1*pi/5)*cos(w1*t*k1)/pi./k1;end
for k2=1:500
    s2=s2+2*sin(k2*pi/5)*cos(w1*t*k2)/pi./k2;end
% subplot(143);fplot(s1);
% subplot(144);fplot(s2);
%%
syms t;f1=sin(t)/t;f2=pi*sym('heaviside(t+1)-heaviside(t-1)');
F1=simplify(fourier(f1));F2=simplify(fourier(f2));
% subplot(141);fplot(f1,[-10 10]);
% subplot(142);fplot(F1,[-2 2]);
% subplot(143);fplot(f2,[-2 2]);
% subplot(144);fplot(F2,[-10 10]);
%%
clear;clc
syms t w s;
f1=exp(-t)*sin(2*t)*heaviside(t);F1=simplify(fourier(f1));
F2=2*sinc(2*w);f2=simplify(ifourier(F2));
f3=exp(-2*t)*cos(3*t);L3=laplace(f3);
L4=(s+2)/(s^3+9*s); f4=ilaplace(L4);
%%
clear clc
w=-10:0.03:10;Fw=2.*sin(w).*exp(-1i*w)./w;
subplot(121);plot(w,abs(Fw));
% syms t;f1=sym('heaviside(t)-heaviside(t-2)');%利用傅里叶变换函数画出图像
% F1=simplify(fourier(f1));                    %此方法画出的图像在w=0出有区别，存疑
% subplot(122);fplot(abs(F1),[-10 10]);
x=-0:0.07:5;y=-10:0.07:10;
[x,y]=meshgrid(x,y);z=x+1i*y;
z=abs((1-exp(-2*z))./z);
% subplot(122);mesh(x,y,z);surf(x,y,z);
axis([-0 5 -10 10 0 2]);

%%
n=-10:10;w1=pi/4;
n1=-10:-1;ft1=sin(pi*n1/4)./(pi*n1);
n2=1:10;ft2=sin(pi*n2/4)./(pi*n2);
ft=[ft1 1/4 ft2];
fn=abs(ft);phase=angle(ft);
% subplot(141);stem(n,fn);
% subplot(142);stem(n,phase);
syms t;s1=0.25;s2=0.25;
for k1=1:5
    s1=s1+2*sin(k1*pi/4)*cos(w1*t*k1)/pi./k1;end
for k2=1:500
    s2=s2+2*sin(k2*pi/4)*cos(w1*t*k2)/pi./k2;end
% subplot(143);fplot(s1,[-10 10]);
% subplot(144);fplot(s2,[-10 10]);
%%
clear clc
syms t;
x1=exp(-t+1)*cos(0.5*t)*heaviside(t);F1=simplify(fourier(x1));
x2=exp(-t)*sin(2*t)/2/t;F2=simplify(fourier(x2));
syms w;
F3=sinc(2*w)+sinc(pi*w);f3=simplify(ifourier(F3));
F4=cos(w)*(heaviside(w+2)-heaviside(w-2));f4=simplify(ifourier(F4));
%%
clear clc
syms t;
x=exp(-3*t*1i)*(heaviside(t+2)-heaviside(t-2));
F=simplify(fourier(x));
subplot(121);fplot(abs(F));axis([-10 10 0 4]);
subplot(122);fplot(angle(F),[-10 10]);axis([-10 10 0 4]);

%%
syms w;H=sym('exp(-2*w*1i)*(heaviside(w+4)-heaviside(w-4))');
h=simplify(ifourier(H));
% subplot(121);fplot(h);grid;
% axis([-6 6 -0.3 1.4]);g=int(h);
% subplot(122);fplot(g);grid;
%%
n1=-8:4/3:8;f1=sinc(n1);t1=-8:0.1:8;f11=sinc(t1);
% subplot(221);hold on;stem(n1,f1);plot(t1,f11,':');
n2=-8:4/5:8;f2=sinc(n2);t2=-8:0.1:8;f12=sinc(t2);
% subplot(222);hold on;stem(n2,f2);plot(t2,f12,':');
x1=-4.5*pi:0.001:4.5*pi;d1=-4.5*pi:1.5*pi:4.5*pi;
% subplot(223);y1=pulstran(x1+0.75*pi,d1,'rectpuls',0.5*pi);plot(x1/pi,y1+1);
%pulstran产生脉冲串，x1+0.75*pi为时间轴，d1为采样间隔，rectpuls使其图像为矩形图像
%0.5*pi确定单个图像信号的宽度
x2=-5*pi:0.001:5*pi;d2=-5*pi:2.5*pi:5*pi;
% subplot(224);y2=pulstran(x2,d2,'rectpuls',2*pi);plot(x2/pi,y2);
%%
clear clc;
syms t;g1='1+0.6*cos(pi*t)';g2='1+1.4*cos(pi*t)';
subplot(221);fplot(g1);subplot(222);fplot(g2,[0 4]);
s1='cos(500*t)*(1+0.6*cos(pi*t))';s2='cos(500*t)*(1+1.4*cos(pi*t))';
subplot(223);fplot(s1,[0 4]);subplot(224);fplot(s2,[0 4]);

%%
syms w;H=sym('exp(-w*1i)*(heaviside(w+3)-heaviside(w-3))');
h=simplify(ifourier(H));
% subplot(121);fplot(h);
g=int(h);
% subplot(122);fplot(g);
%%
syms t;g1=sym('1+0.3*sin(2*pi*t)');g2=sym('1+1.3*sin(2*pi*t)');
s1=sym('cos(400*pi*t)*(1+0.3*sin(2*pi*t))');
s2=sym('cos(400*pi*t)*(1+1.3*sin(2*pi*t))');
figure,subplot(121);fplot(g1);axis([0 4 0.7 1.3]);
subplot(122);fplot(g2,[0 4]);
figure,subplot(121);fplot(s1,[0 4]);
subplot(122);fplot(s2,[0 4]);

%%
num=[2 1];den=[1 3 5];
sys=tf(num,den);w=0:0.5:20;
subplot(141);pzmap(sys);%给出系统函数，利用pzmap画出零极点图
subplot(142);impulse(sys);grid;%求冲激响应
h=freqs(num,den,w);%freqs对num及den组成的系统函数求频率响应，w为频率值范围
h1=abs(h);h2=angle(h);
subplot(143);plot(w,h1);grid;
subplot(144);plot(w,h2*180/pi);grid;
%%
w=0:0.5:15;
num1=2;num2=[1 0];num3=[2 0];num4=[1 0 5];
den1=[1 2];den2=[1 2];den3=[1 2 5];den4=[1 2 5];
h1=freqs(num1,den1,w);subplot(141);plot(w,abs(h1));grid;
h2=freqs(num2,den2,w);subplot(142);plot(w,abs(h2));grid;
h3=freqs(num3,den3,w);subplot(143);plot(w,abs(h3));grid;
h4=freqs(num4,den4,w);subplot(144);plot(w,abs(h4));grid;
%%
clear clc;
w=0:0.05:4*pi;
num1=[1 1 -2];num2=[1 3 2];num3=[1 -1 -2];num4=[1 -3 2];
den=[1 4 13];
h1=freqs(num1,den,w);subplot(141);plot(w/pi,angle(h1)*180/pi);grid;title('|H1(j\Omega)|');xlabel('\Omega/pi');
h2=freqs(num2,den,w);subplot(142);plot(w/pi,angle(h2)*180/pi);grid;title('|H2(j\Omega)|');xlabel('\Omega/pi');
h3=freqs(num3,den,w);subplot(143);plot(w/pi,angle(h3)*180/pi);grid;title('|H3(j\Omega)|');xlabel('\Omega/pi');
h4=freqs(num4,den,w);subplot(144);plot(w/pi,angle(h4)*180/pi);grid;title('|H4(j\Omega)|');xlabel('\Omega/pi');

%%
w=0:0.5:20;
%num1=1;den1=[1 0];num2=[2 3];den2=[1 2 6];
num=[2 3];den=[1 2 6 0];
%sys1=tf(num1,den1);sys2=tf(num2,den2);
%sys=sys1*sys2;
sys=tf(num,den);
subplot(141);pzmap(sys);
subplot(142);impulse(sys);grid;
h=freqs(num,den,w);
subplot(143);plot(w,abs(h));grid;
subplot(144);plot(w,angle(h)*180/pi);grid;
%%
w=0:0.5:15;
num1=[2 0];num2=[1 -1];num3=[2 1];num4=[1 0 4];
den1=[1 3];den2=[1 3];den3=[1 3 7];den4=[1 3 7];
h1=freqs(num1,den1,w);subplot(141);plot(w,abs(h1));grid;
h2=freqs(num2,den2,w);subplot(142);plot(w,abs(h2));grid;
h3=freqs(num3,den3,w);subplot(143);plot(w,abs(h3));grid;
h4=freqs(num4,den4,w);subplot(144);plot(w,abs(h4));grid;
%%
w=0:0.5:15;den=[1 6 10];
num1=[1 2.5 -1.5];num2=[1 3.5 1.5];num3=[1 -2.5 -1.5];num4=[1 -3.5 1.5];
sys1=tf(num1,den);sys2=tf(num2,den);sys3=tf(num3,den);sys4=tf(num4,den);
subplot(141);pzmap(sys1);subplot(142);pzmap(sys2);
subplot(143);pzmap(sys3);subplot(144);pzmap(sys4);
% h1=freqs(num1,den,w);subplot(141);plot(w/pi,angle(h1)*180/pi);grid;title('|H1(j\Omega)|');xlabel('\Omega/pi');
% h2=freqs(num2,den,w);subplot(142);plot(w/pi,angle(h2)*180/pi);grid;title('|H2(j\Omega)|');xlabel('\Omega/pi');
% h3=freqs(num3,den,w);subplot(143);plot(w/pi,angle(h3)*180/pi);grid;title('|H3(j\Omega)|');xlabel('\Omega/pi');
% h4=freqs(num4,den,w);subplot(144);plot(w/pi,angle(h4)*180/pi);grid;title('|H4(j\Omega)|');xlabel('\Omega/pi');

%%
z=0.2;p=[0.8*exp(pi*1i/4);0.8*exp(-pi*1i/4)];
k=1;subplot(121);zplane(z,p);%给出系统零极点位置，并使用zplane函数将其画出零极点图
[b,a]=zp2tf(z,p,k);%zp2tf将极点和零点转换为传递函数形式。
%subplot(122);zplane(b,a);%给出传递函数也可使用zplane画出
%subplot(122);impz(b,a,20);%impz画出系统的冲激响应，20表示画出n=0~19点的值
%%
syms n z;x1=sym('(1/3)^n');x2=cos(n*pi/2);
z1=ztrans(x1);%ztrans z变换
z2=ztrans(x2);
z3=z/((z+1)*(z-2));z4=(z.^2)/((z-0.5)*(z-0.25));
x3=iztrans(z3);%iztrans z逆变换
x4=iztrans(z4);
%%
n1=0:10;h1=(0.5).^n1;h2=(0.8).^n1;x=ones(1,11);%ones,表示返回一个1x11的数组，x=u[n]-u[n-11]
h=conv(h1,h2);y=conv(h,x);
subplot(221);stem(0:length(y)-1,y,'filled');
b=[1 0 0];a=[1 -1.3 0.4];subplot(222);zplane(b,a);
[H, w]=freqz(b,a);%注意！求离散系统的频响特性使用的是freqz，连续系统使用的是freqs
subplot(223);plot(w/pi,abs(H));
subplot(224);plot(w/pi,angle(H)*180/pi);

%%
clear clc;
syms n z;x1=2.^n;x2=sym('sin(n*pi/3)');
z1=ztrans(x1);
z2=ztrans(x2);
z3=1/((z+1)*(z+2));z4=z/((z-0.4)^2*(z-0.5));
x3=iztrans(z3);
x4=iztrans(z4);
%%
n=0:10;h1=(0.6).^n;h2=(0.9).^n;h=h1+h2;
x=ones(1,10);y=conv(h,x);subplot(221);stem(0:length(y)-1,y);
b=[2 -1.5 0];a=[1 -1.5 0.54];
subplot(222);zplane(b,a);
[H,w]=freqz(b,a);subplot(223);plot(w/pi,abs(H));
subplot(224);plot(w/pi,angle(H)*180/pi);

%%
b1=[1 2];a1=[1 3 2];[A1,B1,C1,D1]=tf2ss(b1,a1);
b2=[1 -3];a2=[1 2 3 6];[A2,B2,C2,D2]=tf2ss(b2,a2);
%%
clear;
clc;
A=[1 0;1 -3];B=[1 0;0 1];C=[1 -1;0 -1];D=[1 1;1 0];
[b1,a1]=ss2tf(A,B,C,D,1);
[b2,a2]=ss2tf(A,B,C,D,2);
r0=[1 -1];dt=0.01;t=0:dt:2;
x(:,1)=ones(length(t),1);x(:,2)=exp(-3*t);
sys=ss(A,B,C,D);y=lsim(sys,x,t,r0);
plot(t,y(:,1),'r');text(1,6,'y1(t)');
hold on;plot(t,y(:,2));text(1,1,'y2(t)');
%%
A=[-1 3;-2 4];B=[2;1];C=[-1 2;1 -1];D=0;
r0=[1;-1];N=10;x=ones(1,N);sys=ss(A,B,C,D,[]);y=lsim(sys,x,[],r0);
subplot(211);y1=y(:,1);stem((0:N-1),y1);
subplot(212);y2=y(:,2);stem((0:N-1),y2);

%%
b1=[2 3];a1=[1 2 3 4];[A1,B1,C1,D1]=tf2ss(b1,a1);
b2=[1 -3 0];a2=[1 2 2 4];[A2,B2,C2,D2]=tf2ss(b2,a2);
%%
clear;
clc;
A=[0 1;-6 5];B=[0 2;1 1];C=[2 0;1 -1];D=[1 -1;0 1];dt=0.01;t=0:dt:2;
r0=[2;-1];x(:,1)=ones(length(t),1);x(:,2)=exp(-0.5*t);
sys=ss(A,B,C,D);y=lsim(sys,x,t,r0);
plot(t,y(:,1),'r');text(1,6,'y1');
hold on;plot(t,y(:,2));text(1,1,'y2');
%%
clear;
A=[-5 -1;3 -1];B=[-1;2];C=[2 1;-1 -3];D=0;
r0=[2;1];N=10;n=0:N;x=exp(-n);
sys=ss(A,B,C,D);y=lsim(sys,x,n,r0);
subplot(121);stem(n,y(:,1));
subplot(122);stem(n,y(:,2));

%第一章离散时间与系统
%%
%求实部和虚部
n0=-1;n2=10;
n=n0:n2;
x=exp((0.4+0.6*1j)*n);
% subplot(211);stem(n,real(x),'.');%求实部
% subplot(212);stem(n,imag(x),'filled');%求虚部
%%
%求将x延拓5个周期的序列
clc;clear all
x=[1,2,3,4];N=length(x);k=5;
nx=0:N-1;
ny=0:(k*N-1);
y=x(mod(ny,N)+1);
% subplot(211);stem(nx,x,'filled');
% subplot(212);stem(ny,y,'filled');
%%
%求两个序列的和以及乘积
clc;clear all
x1=[1,3,5,7,6,4,2,1];ns1=-3;
x2=[4,0,2,1,-1,3];ns2=1;
nf1=ns1+length(x1)-1;
nf2=ns2+length(x2)-1;
n1=ns1:nf1;
n2=ns2:nf2;
n=min(ns1,ns2):max(nf1:nf2);
y1=zeros(1,length(n));y2=y1;
y1(find((n>=ns1)&(n<=nf1)==1))=x1;
y2(find((n>=ns2)&(n<=nf2)==1))=x2;
ya=y1+y2;ym=y1.*y2;
% subplot(221);stem(n1,x1,'filled');
% subplot(222);stem(n2,x2,'filled');
% subplot(223);stem(n,ya,'filled');
% subplot(224);stem(n,ym,'filled');
%%
%两个序列卷积
clc;clear all
x=[1,2,3,-1,-2];nx=-1:3;
h=[2,2,1,-1,4,-2];nh=-3:2;
ny=(nx(1)+nh(1)):(nx(end)+nh(end));
y=conv(x,h);
stem(ny,y,'filled');
%%
%差分方程求解零状态响应
b=1;a=[1,-0.5];x=[1,0,0,0,0];
y=filter(b,a,x);
%若要求解全响应，给出初始状态ys=[y(-1),y(-2),....,y(-N)];xs=[x(-1),x(-2),....,x(-N)]
%使用xic=filtic(b,a,ys,xs);y=filter(b,a,x,xic);求解

%第二章离散时间傅里叶变换
%%
%求x的DTFT，并画出幅频特性及相频特性
clc;clear all
n=0:4;x=[2 3 4 3 2];
w=0:pi/500:2*pi;X=x*(exp(-1j)).^(n'*w);
magX=abs(X);angX=angle(X);
% subplot(221),stem(n,x,'filled');
% subplot(223),plot(w/pi,magX);
% subplot(224),plot(w/pi,angX);
%若直接使用函数freqz，将其前三行改为以下语句
% clc;clear all
% n=0:4;a=1;b=[2 3 4 3 2];
% N=1000;
% [X,w]=freqz(b,a,N);%其中N表示在w(0~2pi)范围内取1000点
%%
%出现0/0情况时，令输入变量+1e-10
clc;clear all
N=5;M=2*N+1;n=-20:20;
x=[zeros(1,15),ones(1,M),zeros(1,15)];
omega=[-pi:0.01*pi:pi]+1e-10;
X=sin(0.5*M*omega)./sin(0.5*omega);
% subplot(211);stem(n,x,'.');axis([-20 20 -0.2 1.5]);
% subplot(212);plot(omega,X);axis([-pi pi min(X),max(X)]);grid;
%%
%求零极点图以及幅频相频图
clc;clear all
b=[1 -1.8 -1.44 0.64];a=[1 -1.64853 1.03882 -0.288];
rp=roots(a) %求极点
rz=roots(b) %求零点
[H,w]=freqz(b,a,1024,'whole');
figure(1)
zplane(b,a);
figure(2)
subplot(211);plot(w,abs(H));
subplot(212);plot(w,angle(H));
%%
%判断系统稳定性
clc;clear all
a=[1 -1.4 0.84 -0.288];
zp=roots(a);
zpm=max(abs(zp));
if zpm<1
    disp('系统稳定');
else
    disp('系统不稳定');
end

%第三章DFT
%%
%圆周移位
clc;clear all
x=[1 2 3 4 5];
y1=cirshift(x,3,5);
y2=cirshift(x,-3,6);
%%
%圆周翻褶
clc;clear all
x=[2 3 4 5 6];N=8;
x=[x,zeros(1,N-length(x))];nx=0:N-1;
y=x(mod(-nx,N)+1);
subplot(211);stem(0:N-1,x);
subplot(212);stem(0:N-1,y);
%%
%求N=4的DFT矩阵
w4=dftmtx(4);
%%
%N点IDFT可以表示为DFT矩阵取共轭再乘以1/N
w4I=conj(dftmtx(4))/4;
%%
%DFT的列向量X等于DFT矩阵左乘离散时域x列向量
% x=[4;3;2;1];
% X=dftmtx(4)*x
%求IDFT
X=[10;2-2i;2;2+2i];
x=conj(dftmtx(4))/4*X;
%%
%用圆周卷积计算啊线性卷积
clc;clear all
x1=[2 4 3 1];x2=[2 1 3];
y=circonvt(x1,x2,4);
%在频域计算圆周卷积
% X1=fft(x1,N);
% X2=fft(x2,N);
% Y=X1*X2;
% y=ifft(Y,N);

%第四章 FFT
%%
clc;clear all
x=[2 1 -1 2 3];nx=0:4;K=128;dw=2*pi/K;
k=floor((-K/2+0.5):(K/2-0.5));
X=x*exp(-1j*dw*nx'*k);
subplot(221);plot(k*dw,abs(X));hold on;
Xd=fft([2 1 -1 2 3]);stem([0:4]*2*pi/5,abs(Xd));
title('5点序列的DTFT和FFT');
Xd1=fftshift(Xd);
subplot(222);plot(k*dw,abs(X));hold on;
stem([-2:2]*2*pi/5,abs(Xd1));
ylabel('幅度响应');title('FFT移位后');
subplot(223);plot(k*dw,angle(X));
ylabel('相位响应');title('FFT移位后');
%%
%使用FFT算法求解两有限长序列的线性卷积
clc;clear all
x=[2 1 3 2 1 5 1];h=[1 2 -1 -3];
N=length(x)+length(h)-1;
n=0:N-1;
x=[x,zeros(1,N-length(x))];
h=[h,zeros(1,N-length(h))];
X=fft(x);H=fft(h);
Y=X.*H;y=ifft(Y);
% subplot(221);stem(n,x);title('x(n)');
% subplot(222);stem(n,h);title('h(n)');
% subplot(223);stem(n,y);title('y(n)');
%%
%求长度为9的x的频谱DTFT
clc;clear all
c=[9,16,32,512];T=0.4;
for i=1:4
    L=c(i);D=2*pi/(L*T);
    x=[ones(1,5),zeros(1,L-9),ones(1,4)];
    k=floor(-(L-1)/2:(L-1)/2);
    X=fftshift(fft(x,L));
%     subplot(2,2,i);plot(k*D,real(X));
    axis([min(k*D),max(k*D),-5,10]);
    Str=['N=',num2str(L)];title(Str);
end
%%
%求x的频谱
clear,clc
T0=[0.05,0.02,0.01,0.01];
L0=[10,10,10,20];
for i=1:4
    T=T0(i);N=L0(i)/T;
    D=2*pi/(N*T);
    n=0:N-1;x=exp(-0.02*n*T).*cos(6*pi*n*T)+2*cos(14*pi*n*T);
    k=floor(-(N-1)/2:(N-1)/2);
    X=T*fftshift(fft(x));
    subplot(2,2,i);plot(k*D,abs(X));
    axis([min(k*D),max(k*D),0,inf]);
    str=['T=',num2str(T),'N=',num2str(N)];
    title(str);
end
%%
% clc;clear;close all
% N=input('N=');T=0.05;n=1:N;
% D=2*pi/(N*T);
% xa=cos(10*n*T);
% Xa=T*fftshift(fft(xa,N));
% k=floor(-(N-1)/2:(N-1)/2);
% TITLE=sprintf('N=%i,L=%i',N,N*T);
% plot(k*D,abs(Xa));
% axis([-20 20 0 max(abs(Xa)+2)]);
% title(TITLE);
%%
clc;clear;close all
T=1;c=[4,8];
for i=1:2
    N=c(i);D=2*pi/N;
    kl=floor(-(N-1)/2:-1/2);
    kh=floor(0:(N-1)/2);
    w=[kh,kl]*D;
    X=3+exp(-1i*w)-2*exp(-1i*3*w)+exp(-1i*4*w);
    x=ifft(X);
    subplot(1,2,i);stem(T*[0:N-1],x);
    axis([-1 9 -3 5]);
    str=['N=',num2str(N)];title(str);
end
%%
clc;clear;close all
wc=5;T=0.1*pi/wc;
N=100*2*pi/wc/T;
D=2*pi/(N*T);w=[0:N-1]*D;
M=floor(wc/D);
H=[ones(1,M+1),zeros(1,fix(N-2*M-1)),ones(1,M)];
h=ifft(H/T);
plot([0:N-1]*T,h);

%第七章无限长数字滤波器设计
%%
%巴特沃斯滤波器
clc;clear all
OmegaP=2*pi*3000;OmegaS=2*pi*6000;Rp=2;As=30;
N=ceil(log10((10^(As/10)-1)/(10^(Rp/10)-1))/(2*log10(OmegaS/OmegaP)));
OmegaC=OmegaP/((10^(Rp/10)-1)^(1/(2*N)));
[z0,p0,k0]=buttap(N);
p=p0*OmegaC;a=real(poly(p));
k=k0*OmegaC^N;b0=real(poly(z0));
b=k*b0;
w0=[OmegaP,OmegaS];
[H,w]=freqs(b,a);
Hx=freqs(b,a,w0);
dbHx=-20*log10(abs(Hx)/max(abs(H)));
plot(w/(2*pi)/1000,20*log10(abs(H)));
xlabel('f(kHz)');ylabel('dB');axis([-1 12 -55 1]);
set(gca,'xtickmode','manual','xtick',[0 1 2 3 4 5 6 7 8 9 10]);
set(gca,'ytickmode','manual','ytick',[-50 -40 -30 -20 -10 0]);grid;
%%
%切比雪夫低通滤波器
clc;clear;close all;
OmegaP=2*pi*3000;OmegaS=2*pi*6000;Rp=2;As=30;
g=sqrt((10^(As/10)-1)/(10^(Rp/10)-1));
OmegaR=OmegaS/OmegaP;
N=ceil(log10(g+sqrt(g*g-1))/log10(OmegaR+sqrt(OmegaR*OmegaR-1)));
OmegaC=OmegaS;
[z0,p0,k0]=cheb2ap(N,As);
a0=real(poly(p0));
aNn=a0(N+1);p=p0*OmegaC;a=real(poly(p));
aNu=a(N+1);b0=real(poly(z0));M=length(b0);
bNn=b0(M);z=z0*OmegaC;b=real(poly(z));
bNu=b(M);k=k0*(aNu*bNn)/(aNn*bNu);
b=k*b;
w0=[OmegaP,OmegaS];
[H,w]=freqs(b,a);
Hx=freqs(b,a,w0);
dbHx=-20*log10(abs(Hx)/max(abs(H)));
plot(w/(2*pi)/1000,20*log10(abs(H)));
xlabel('f(kHz)');ylabel('dB');axis([-1 12 -55 1]);
set(gca,'xtickmode','manual','xtick',[0 1 2 3 4 5 6 7 8 9 10]);
set(gca,'ytickmode','manual','ytick',[-50 -40 -30 -20 -10 0]);grid;
%%

function y=circonvt(x1,x2,N)
%求x1与x2的圆周卷积
if length(x1)>N
    error('N长度必须大于x1长度');
end
if length(x2)>N
    error('N长度必须大于x2长度');
end
x1=[x1 zeros(1,N-length(x1))];
x2=[x2 zeros(1,N-length(x2))];
m=0:N-1;
x2=x2(mod(-m,N)+1);
H=zeros(N,N);
for n=1:N
    H(n,:)=cirshift(x2,n-1,N);
end

y=x1*H';

function y=cirshift(x,m,N)
%长度为N的x序列作m点圆周移位
%x为输入序列，长度应小于N
%m=移位点数
%N=圆周移位的长度
if length(x)>N
    error('N长度必须大于x长度')
end
x=[x zeros(1,N-length(x))];
n=0:N-1;
n=mod(n-m,N);
y=x(n+1);