DSP实验Matlab程序实例

1、2离散时间信号（序列）的产生利用MATLAB产生和绘制下列有限长序列：单位脉冲序列，单位阶跃序列，矩形序列x=-10:10;y=; for i=1:21 if x(i)=0 y(i)=1; else y(i)=0; endendsubplot(3,1,1);%在子图画单位冲激序列stem(x,y,'*');title('冲激'); xlabel('n'); hold on; for i=1:21 if x(i)<0 y(i)=0; else y(i)=1; endendsubplot(3,1,2); %在子图画单位阶跃序列stem(x,y,

2、'*'); title('阶跃'); xlabel('n');hold on; for i=1:21 if x(i)<-4 y(i)=0; else if x(i)>4 y(i)=0; else y(i)=1; end end end subplot(3,1,3);%在子图画矩形脉冲序列stem(x,y,'*'); title('矩形'); xlabel('n'); hold on在一幅图上绘出曲线，A=2, =0.5, f=2Hz，A1=1, A2=0.5, A3=0.2，f=2Hz。

3、t=0:0.5:720*pi/180;x1=5*sin(pi*t/5+pi/3);x2=5*cos(t-pi/3);x3=2*sin(3*t).*cos(2*t-pi/3);subplot(3,1,1); %在子图中画出一系列正弦余弦曲线plot(t,x1,'-red',t,x2,'-gr',t,x3,'-bl');title('正弦余弦曲线');xlabel('t');hold on;A=2;a=0.5;f=2;x4=A*exp(-a*t).*sin(2*pi*f*t);subplot(3,1,2);%在子图中画

4、出正弦衰减信号plot(t,x4)title('正弦衰减信号');xlabel('t');hold on;A1=1;A2=0.5;A3=0.2;f=2;x5=A1*sin(2*pi*f*t)+A2*sin(2*pi*2*f*t)+A3*sin(2*pi*3*f*t);subplot(3,1,3);%在子图中画出谐波信号plot(t,x5)title('谐波信号');xlabel('t');hold on;3序列的运算生成下列序列：(1) 利用MATLAB编程完成上述两序列的卷积，并绘制运算后序列的波形。x=1 2 3 4 5;h=

5、1 2 1 2;n1=0 1 2 3 4;n2=2 3 4 5;y=conv(x,h); %计算序列x与h的卷积和序列yn0=n1(1)+n2(1); %计算卷积和序列y的起点位置n3=length(x)+length(h)-2+n0; %计算卷积和序列y的终止位置n=n0:1:n3; %确定卷积和y非零样值的时间向量subplot(3,1,1);stem(n1,x,'.'); %在子图绘x(n)title('x(n)');xlabel('n');ylabel('x(n)');subplot(3,1,2);stem(n2,h,&#

6、39;.'); %在子图绘h(n)title('h(n)');xlabel('n');ylabel('h(n)');subplot(3,1,3);stem(n,y,'.'); %在子图绘卷积和y(n)title('y(n)=x(n)*h(n)');xlabel('n');ylabel('y(n)');（2）编写一个完成卷积的函数，输入是两个序列，输出是结果。function y,n=jsjuanji(x,n1,h,n2)%计算序列卷积和y(n)=x(n)*h(n)% y: 卷

7、积和y(n)对应的非零样值向量% n：卷积和y(n)对应的时间向量% x: x(n)对应的非零样值向量% n1: x(n)对应的时间向量% h: h(n)对应的非零样值向量% n2: h(n)对应的时间向量n3=n1(1)+n2(1); %卷积和序列y的起始位置n4=length(x)+length(h)-1;%卷积和序列y的长度m=0;for i=1:n4 %求卷积和y if i>length(x) x(i)=0; end for j=1:i if j>length(h) h(j)=0; end m=m+x(i+1-j)*h(j); end y(i)=m; m=0;endn1=n

8、1(1):1:n1(1)+n4-1; n1: % x(n)对应的新的时间向量n2=n2(1):1:n2(1)+n4-1; % h(n)对应的新的时间向量n=n3:1:n3+n4-1; %y(n)对应的时间向量subplot(3,1,1);stem(n1,x,'.'); %在子图绘x(n)title('x(n)');xlabel('n');ylabel('x(n)');subplot(3,1,2);stem(n2,h,'.'); %在子图绘h(n)title('h(n)');xlabel('n

9、');ylabel('h(n)');subplot(3,1,3);stem(n,y,'.'); %在子图绘卷积和y(n)title('y(n)=x(n)*h(n)');xlabel('n');ylabel('y(n)');4采样定理的研究分别令采样周期Ts为不同值，绘出不同采样周期下x(t) =sin(t)的频谱（直接用FFT函数），观察频谱混叠现象和防止混叠的采样周期。fs=2; %设定采样频率T=1/fs; %采样周期N=128; %采样点数n=0:N-1; t=n*T; %采样时间点x=sin(t);

10、 %生成正弦信号subplot(4,2,1);plot(t,x);%画正弦信号的时域波形xlabel('t');ylabel('x');title('正弦信号x=sin(t)时域波形'); y=fft(x,N); %进行fft变换mag=sqrt(y.*conj(y); %求幅值f=(0:N-1)*fs/N; %进行对应的频率转换subplot(4,2,2);plot(f,mag); %画频谱图axis(-10,50,0,80);xlabel('频率');ylabel('幅值');title('正弦信号x=

11、sin(t)幅频谱图 fs=2');hold on;fs=10; %设定采样频率T=1/fs; %采样周期N=1024; %采样点数n=0:N-1;t=n*T; %采样时间点x=sin(t); %生成正弦信号subplot(4,2,3);plot(t,x);%画正弦信号的时域波形xlabel('t');ylabel('x');title('正弦信号x=sin(t)时域波形');y=fft(x,N); %进行fft变换mag=sqrt(y.*conj(y); %求幅值f=(0:N-1)*fs/N; %进行对应的频率转换subplot(4,2

12、,4);plot(f,mag); %画频谱图axis(-10,50,0,400);xlabel('频率');ylabel('幅值');title('正弦信号x=sin(t)幅频谱图 fs=10');hold on;fs=20; %设定采样频率T=1/fs; %采样周期N=1024; %采样点数n=0:N-1;t=n*T; %采样时间点x=sin(t); %生成正弦信号subplot(4,2,5);plot(t,x);%画正弦信号的时域波形xlabel('t');ylabel('x');title('正弦信号

13、x=sin(t)时域波形');y=fft(x,N); %进行fft变换mag=sqrt(y.*conj(y); %求幅值f=(0:N-1)*fs/N; %进行对应的频率转换subplot(4,2,6);plot(f,mag); %画频谱图axis(-10,50,0,400);xlabel('频率');ylabel('幅值');title('正弦信号x=sin(t)幅频谱图 fs=20');hold on;fs=30; %设定采样频率T=1/fs; %采样周期N=1024; %采样点数n=0:N-1; t=n*T; %采样时间点x=sin(

14、t); %生成正弦信号subplot(4,2,7);plot(t,x);%画正弦信号的时域波形xlabel('t');ylabel('x');title('正弦信号x=sin(t)时域波形'); y=fft(x,N); %进行fft变换mag=sqrt(y.*conj(y); %求幅值f=(0:N-1)*fs/N; %进行对应的频率转换subplot(4,2,8);plot(f,mag); %画频谱图axis(-10,50,0,400);xlabel('频率');ylabel('幅值');title('正弦

15、信号x=sin(t)幅频谱图 fs=30');五、思考题1如何产生方波信号序列和锯齿波信号序列？答：方波信号由square函数可得，锯齿波信号由循环语句产生。%方波T=0:0.001:2*pi; %方波的时间向量y=square(6*T); %方波w=6,周期为T=2*pi/6subplot(2,1,1);plot(T,y); %画方波title('方波');axis(T(1)-1 T(end)+1 -2 2); %坐标轴区域hold on%锯齿波x=-10:10;y=;for i=1:21 %rem为取余函数，此循环使得y(i)交替取值1和-1if rem(i,2)=

16、0 y(i)=-1;else y(i)=1;endendsubplot(2,1,2);plot(x,y);title('锯齿波');1编写一个 DFT 计算程序，要求该程序具有正变换及反变换的双重功能。令 x(n)=sin(n)n(-1)n, n=0,1,.,31, 这是一个震荡信号，由它来验证所编程序是否正确（即：对 x(n)作DFT 得 X(k)，再对 X(k)作反变换，检验是否得到同样的 x(n)）。function X,x2=bianhuan(x,N)%DFT X=zeros(1,N);for m=1:N for n=1:N X(m)=X(m)+x(n)*exp(-j*

17、2*pi*n*m/N); endend %IDFT x2=zeros(1,N);for m=1:N for n=1:N x2(m)=x2(m)+X(n)*exp(j*2*pi*n*m/N); end x2(m)=x2(m)/N;endclearN=32;x=;for p=1:N x(p)=sin(p)*p*(-1)pendX,x2=bianhuan(x,N) n=0:1:N-1;subplot(3,1,1);stem(n,x,'*');title('original');subplot(3,1,2);stem(n,X,'*');title(

18、9;after DFT');subplot(3,1,3);stem(n,x2,'*');title('after IDFT');2. 序列 x(n)的 N 点 DTFT 的物理意义是对 X(ej)的在0，2上进行 N 点等间隔采样clear;N1=1024;w1=(0:2*pi/1024:2*pi);X1=(1-exp(-j*4*w1)./(1-exp(-j*w1); N2=8;for i=1:N2 X2(i)=X1(i*N1/N2); w2(i)=2*pi*i/N2;end N3=16;for i=1:N3 X3(i)=X1(i*N1/N3); w3(

19、i)=2*pi*i/N3;end subplot(3,2,1); plot(w1,abs(X1);titlesubplot(3,2,2), plot(w1,angle(X1);titlesubplot(3,2,3); stem(w2,abs(X2);titlesubplot(3,2,4), stem(w2,angle(X2);titlesubplot(3,2,5); stem(w3,abs(X3);titlesubplot(3,2,6), stem(w3,angle(X3);title 4. 混叠现象研究衰减正弦信号 x(t)=1+sin(7f0t)cos(2f0t), f0=120Hz, f

20、s=200Hz,采样点数 N=64，用 DFT 计算信号的幅频和相频，观察混叠现象clear;f0=120;N=64; fs1=200;t1=(0:1:N-1)/fs1;x1=(1+sin(7*pi*f0*t1).*cos(2*pi*f0*t1);y1=fft(x1); fs2=1500;t2=(0:1:N-1)/fs2;x2=(1+sin(7*pi*f0*t2).*cos(2*pi*f0*t2);y2=fft(x2); subplot(2,1,1);stem(abs(y1); title('fs=200');subplot(2,1,2);stem(abs(y2); title

21、('fs=1500');5. 泄漏现象 令余弦信号 x(t)= cos(2f0t),f0=60Hz,fs=200Hz,采样点数 N=64要求：先用 DFT 计算采样信号 x(n)的频谱 X(k)，显示结果图形。然后将采样信号 x(n)乘上哈明窗函数：x4(t)=0.54-0.46cos(2n/N),n=0,1,2,N-1用 DFT 计算加窗后信号的频谱 XW(k)，分别在屏幕上显示信号幅频 X(k)，窗函数幅频 X4(k)及加窗后信号的频谱 XW(k)，比较频谱 X(k)与 XW(k)的不同，并给出解释。clear;f0=60;fs=200;N=64;t=(0:1:N-1)/f

22、s;x=cos(2*pi*f0*t);X=fft(x); n=0:1:N-1;x1=0.54-0.46.*cos(2*pi*n/N);X1=fft(x1); x2=x.*x1;X2=fft(x2); subplot(3,1,1);stem(abs(X);title；subplot(3,1,2);stem(abs(X1);title;subplot(3,1,3);stem(abs(X2);title;6. 栅栏效应将内容1中的信号补零至L, L 分别为 16，32，64，计算 x(n)的频谱 X(k)，并与 X(k)比较，将频谱图 X(k)与 X(k)显示出来，观察补零的效果。clear;N=8

23、; x=;for p=1:N x(p)=sin(p)*p*(-1)pendX=zeros(1,N);for p=1:N for q=1:N X(p)=X(p)+x(q)*exp(-j*2*pi*p*q/N); endend L1=16; x1=;for p=1:L1 if p<=8 x1(p)=sin(p)*p*(-1)p else x1(p)=0 endend X1=zeros(1,L1); for p=1:L1 for q=1:L1 X1(p)=X1(p)+x1(q)*exp(-j*2*pi*p*q/L1); endend L2=32; x2=;for p=1:L2 if p<

24、=8 x2(p)=sin(p)*p*(-1)p else x2(p)=0 endend X2=zeros(1,L2); for p=1:L2 for q=1:L2 X2(p)=X2(p)+x2(q)*exp(-j*2*pi*p*q/L2); endend L3=64; x3=;for p=1:L3 if p<=8 x3(p)=sin(p)*p*(-1)p else x3(p)=0 endend X3=zeros(1,L3); for p=1:L3 for q=1:L3 X3(p)=X3(p)+x3(q)*exp(-j*2*pi*p*q/L3); endend n=0:1:N-1;n1=0

25、:1:L1-1;n2=0:1:L2-1;n3=0:1:L3-1;subplot(4,1,1);stem(n,X,'*');title;subplot(4,1,2);stem(n1,X1,'*');title('X1:L=16');subplot(4,1,3);stem(n2,X2,'*');title('X2:L=32');subplot(4,1,4);stem(n3,X3,'*');title('X3:L=64');1.基于 FFT 的卷积计算对上述序列做基于 FFT 的卷积计算，

26、做 FFT 时注意延拓序列长度，将结果显示出来，并和直接卷积结果比较.clear;N1=40;N2=30;N=41+31-1; x1=cos(0.025*pi*(0:N1).2);x=x1 zeros(1,N-N1); h1=sin(2.5*pi*(0:N2);h=h1 zeros(1,N-N2); y1=conv(x1,h1); X=fft(x); H=fft(h)Y=X.*H;y=ifft(Y); subplot(2,1,1);stem(y,'.');title('基于FFT卷积计算');subplot(2,1,2);stem(y1,'.')

27、;title(直接卷积计算);2.基于 FFT 的相关计算请计算 h(n)和 x(n)的基于 FFT 的相关函数，将结果显示出来，并和直接相关结果比较。clear;n1=0:128;n2=129:192;n3=193:256;n4=0:31;N=257+32-1; x1=cos(0.25*pi*n1.2)+sin(0.5*pi*n1);x2=cos(pi*n2)+sin(0.5*pi*n2)/6;x3=sin(0.5*pi*n3)/6+(-1).n3;x=x1 x2 x3 zeros(1,(N-257); h1=cos(pi*n4)+sin(0.5*pi*n4)/6;h=h1 zeros(1,

28、(N-32); r1=xcorr(x1 x2 x3,h1); X=fft(x); H=fft(h);R=X.*conj(H);r=ifft(R); subplot(2,1,1);stem(r,'.');title('FFT相关计算');subplot(2,1,2);stem(r1,'.');title('直接计算');3.功率谱计算信号设为：f1=42Hz, f2=45Hz, f3=47Hz,采样频率 fs=300Hz,采样点数 N=256。计算信号的功率谱，将功率谱图显示出来。功率谱定义为 P(k)= |X(k)|2clear;

29、f1=42;f2=45;f3=47;fs=300;N=64;t=(0:1:N-1)/fs;x=sin(2*pi*f1*t)+cos(2*pi*f2*t)+cos(2*pi*f3*t); X=abs(fft(x);P=X.2;stem(P,'.');title(功率谱)1给定差分方程y(n)-y(n-1)+0.9y(n-2)=x(n)（1）计算并画出单位脉冲响应h(n)，n=0,1,127；（2）计算并画出单位阶跃响应s(n)，n=0,1,127；（3）判断系统是否稳定。clear;num=1 0 0;den=1 -1 0.9;subplot(2,1,1);dimpulse(nu

30、m,den,128);title('单位冲激响应'); %绘subplot(2,1,2);dstep(num,den,128);title('单位阶跃响应');y1=dimpulse(num,den,128) %求y2=dstep(num,den,128)以下为定义法% x1=zeros(1,130); %单位冲激函数% x1(3)=1;% y1(1)=0;y1(2)=0;% for k=1:128% y1(k+2)=x1(k+2)+y1(k+1)-0.9*y1(k);% end% % x2=ones(1,130); %单位阶跃函数% y2(1)=0;y2(2)

31、=0;% for k=1:128% y2(k+2)=x2(k+2)+y2(k+1)-0.9*y2(k);% end% % figure% subplot(2,1,1);plot(y1);title('单位脉冲响应') % subplot(2,1,2);plot(y2);title('单位阶跃响应')% y1 % y22给定系统函数（1）画出H(z)零-极点示意图；（2）求出并绘出H(z)的幅频响应、相频响应；（3）求出并绘出该系统的单位抽样响应 h(n)；（4）求系统的单位阶跃响应s(n)。clear;num=-0.2 0;den=1 0 0.7;sys=tf(

32、num,den);figure;pzmap(sys);title('零极点分布图')figure;freqz(num,den);title('频率响应')figure;subplot(2,1,1);dimpulse(num,den);title('单位冲激响应');subplot(2,1,2);dstep(num,den);title('单位阶跃响应');y1=dimpulse(num,den,128) %求y2=dstep(num,den,128)% 以下为定义法% x1=zeros(1,130); %单位冲激函数% x1(3)

33、=1;% y1(1)=0;y1(2)=0;% for k=1:128% y1(k+2)=0.2*x1(k+1)-0.8*y1(k);% end% % x2=ones(1,130); %单位阶跃函数% y2(1)=0;y2(2)=0;% for k=1:128% y2(k+2)=0.2*x2(k+1)-0.8*y2(k);% end% % figure;% subplot(2,1,1);plot(y1);title('单位脉冲响应')% subplot(2,1,2);plot(y2);title('单位阶跃响应')（1）H(z)=H1(z)H2(z)（2）H(z)

34、=H1(z)+H2(z)重复 2 中的内容。clear;num1=1 1;den1=1 -0.83;num2=1 -0.2 -0.8;den2=1 -0.9 -0.83;%H(Z)=H1(Z)H2(Z)num=conv(num1,num2); den=conv(den1,den2);sys=tf(num,den);figure;pzmap(sys);title('H(Z)=H1(Z)H2(Z)零极点分布图')figure;freqz(num,den);title('H(Z)=H1(Z)H2(Z)频率响应')h=dimpulse(num,den,100)s=dst

35、ep(num,den,100)figure;subplot(2,1,1);dimpulse(num,den,100);title('H(Z)=H1(Z)H2(Z)单位冲激响应');subplot(2,1,2);dstep(num,den,100);title('H(Z)=H1(Z)H2(Z)单位阶跃响应');% %H(Z)=H1(Z)+H2(Z)num3=conv(den1,num2)+conv(den2,num1) den3=conv(den1,den2)sys3=tf(num3,den3);figure;pzmap(sys3);title('H(Z)

36、=H1(Z)+H2(Z)零极点分布图')figure;freqz(num3,den3);title('H(Z)=H1(Z)+H2(Z)频率响应')h3=dimpulse(num3,den3,100)s3=dstep(num3,den3,100)figure;subplot(2,1,1);dimpulse(num3,den3,100);title('H(Z)=H1(Z)+H2(Z)单位冲激响应');subplot(2,1,2);dstep(num3,den3,100);title('H(Z)=H1(Z)+H2(Z)单位阶跃响应');4.梳状

37、滤波器零极点和频率特性梳状滤波器系统函数有如下2种类型：分别对N=8，a=0.85，0.9，0.95计算H1(z) H 2 (z) 和的零点、极点及频率特性曲线。clear;%FIR型num1=-1 0 0 0 0 0 0 0 1;den1=1;sys1=tf(num1,den1);figure;freqz(num1,den1);title('FIR型频率响应')%IIR型a1=0.8;a2=0.9;a3=0.98;num2=-1 0 0 0 0 0 0 0 1;den21=-a18 0 0 0 0 0 0 0 1;den22=-a28 0 0 0 0 0 0 0 1;den2

38、3=-a38 0 0 0 0 0 0 0 1;%a=0.8sys21=tf(num2,den21);figure;freqz(num2,den21);title('IIR型频率响应a=0.8')% a=0.9sys22=tf(num2,den22);figure;freqz(num2,den22);title('IIR型频率响应a=0.9')% a=0.98sys23=tf(num2,den23);figure;freqz(num2,den23);title('IIR型频率响应a=0.98')%figure;subplot(2,2,1);pzma

39、p(sys1);title('FIR型零极点分布图')subplot(2,2,2);pzmap(sys21);title('IIR型零极点分布图a=0.8')subplot(2,2,3);pzmap(sys22);title('IIR型零极点分布图a=0.9')subplot(2,2,4);pzmap(sys23);title('IIR型零极点分布图a=0.98')1、分别用冲激响应不变法和双线性变换法将模拟滤波器转换为数字滤波器，fs=1，10，100。要求：a) 写出 H(z)表达式。b) 绘出数字滤波器的幅频响应及相频响应。

40、c) 比较两种映射法的特点clear;%冲击响应不变法%fs=1z2 - 0.22088 zH(Z)=- z2 - 0.18512 z + 0.0067379fs=10 0.1 z2 - 0.089664 zH(Z)= - z2 - 1.5595 z + 0.60653fs=100 0.01 z2 - 0.0098995 zH(Z)=- z2 - 1.9506 z + 0.95123c=0 1 1;d=1 5 6;fs=1;b,a=impinvar(c,d,fs);printsys(b,a,'z');h,f=freqz(b,a,200);figure;subplot(2,3,1

41、);plot(f/pi,abs(h);title('fs=1幅频特性');subplot(2,3,4);plot(f/pi,angle(h);title('fs=1相频特性');fs=10;b,a=impinvar(c,d,fs);printsys(b,a,'z');h,f=freqz(b,a,200);subplot(2,3,2);plot(f/pi,abs(h);title('fs=10幅频特性');subplot(2,3,5);plot(f/pi,angle(h);title('fs=10相频特性');fs=

42、100;b,a=impinvar(c,d,fs);printsys(b,a,'z');h,f=freqz(b,a,200);subplot(2,3,3);plot(f/pi,abs(h);title('fs=100幅频特性');subplot(2,3,6);plot(f/pi,angle(h);title('fs=100相频特性');clear;%双线性变换法%fs=1 0.15 z2 + 0.1 z - 0.05H(Z)= - z2 + 0.2 z - 2.2204e-017fs=10 0.041502 z2 + 0.0039526 z - 0

43、.037549H(Z)= -z2 - 1.5573 z + 0.60474fs=100 0.0049017 z2 + 4.8773e-005 z - 0.0048529H(Z)= -z2 - 1.9506 z + 0.95123c=0 1 1;d=1 5 6;fs=1;b,a=bilinear(c,d,fs);printsys(b,a,'z');h,f=freqz(b,a,200);figure;subplot(2,3,1);plot(f/pi,abs(h);title('fs=1幅频特性');subplot(2,3,4);plot(f/pi,angle(h);

44、title('fs=1相频特性');fs=10;b,a=bilinear(c,d,fs);printsys(b,a,'z');h,f=freqz(b,a,200);subplot(2,3,2);plot(f/pi,abs(h);title('fs=10幅频特性');subplot(2,3,5);plot(f/pi,angle(h);title('fs=10相频特性');fs=100;b,a=bilinear(c,d,fs);printsys(b,a,'z');h,f=freqz(b,a,200);subplot(2

45、,3,3);plot(f/pi,abs(h);title('fs=100幅频特性');subplot(2,3,6);plot(f/pi,angle(h);title('fs=100相频特性');2、设计切比雪夫型高通数字滤波器，令通带截止频率为250Hz，阻带截止频率为50Hz，Rp=3dB,As=35dB,Fs=1KHz。显示滤波器的幅频、相频曲线。要求：a) 确定DF的系统函数H(z)，显示N=256 点的幅频及相频曲线。b) 滤波器频带检验：x(t)=cos(2f1t)+sin(2f2t)，f1=10Hz，f2=300Hz将信号x(t)通过设计好的低通数字

46、滤波器系统，显示滤波器输出的幅频曲线，与信号的幅频比较，分析信号通过滤波器后的频谱变化a) 0.21301 z2 - 0.42603 z + 0.21301H(Z)= -z2 + 0.24825 z + 0.45181clear;%切比雪夫滤波器%wp=0.5;ws=0.1;Rp=3;As=35;N,wn=cheb1ord(wp,ws,Rp,As);b,a=cheby1(N,Rp,wn,'high');printsys(b,a,'z');h,f=freqz(b,a,256);figure;subplot(2,1,1);plot(f/pi,abs(h);title

47、('chebyshev高通滤波器幅频特性');subplot(2,1,2);plot(f/pi,angle(h);title('chebyshev高通滤波器相频特性');%验证%f1=10;f2=300;t=1:4:1024;x=cos(2*f1*t)+sin(2*f2*t);X=fft(x);Y=X'.*h;figure;subplot(2,2,1);plot(abs(X);title('原信号X幅频特性');subplot(2,2,3);plot(angle(X);title('原信号X相频特性');subplot(2

48、,2,2);plot(abs(Y);title('高通滤波后信号Y幅频特性');subplot(2,2,4);plot(angle(Y);title('高通滤波后信号Y相频特性');3、令f1=100Hz,f2=200Hz,f3=400Hz,f4=500Hz,fs=5KHz,设计Butterworth 型带通数字滤波器，Rp=3dB,As=40dB。要求：a) 显示滤波器的幅频、相频曲线。b) 分别令As=20dB, 30dB, 40dB,比较滤波器的阶次、幅频曲线、相频曲线。c) 滤波器频带检验，校验信号同2（b）。clear;%巴特沃兹滤波器（带通）%f1=

49、100;f2=200;f3=400;f4=500;fs=5000;ws1=f1/fs*2;wp1=f2/fs*2;wp2=f3/fs*2;ws2=f4/fs*2;ws=ws1 ws2;wp=wp1 wp2;Rp=3;As=40;N,wn=buttord(wp,ws,Rp,As);b,a=butter(N,wn);h,f=freqz(b,a,256);figure;subplot(2,1,1);plot(f/pi,abs(h);title('butterworth带通滤波器幅频特性');subplot(2,1,2);plot(f/pi,angle(h);title('bu

50、tterworth带通滤波器相频特性');%讨论As变化%As=20;N,wn=buttord(wp,ws,Rp,As);b,a=butter(N,wn);h,f=freqz(b,a,256);figure;subplot(2,3,1);plot(f/pi,abs(h);title('As=20带通滤波器幅频特性');subplot(2,3,4);plot(f/pi,angle(h);title('As=20带通滤波器相频特性');As=30;N,wn=buttord(wp,ws,Rp,As);b,a=butter(N,wn);h,f=freqz(b,a,256);subplot(2,3,2);plot(f/pi,abs(h);title('As=30带通滤波器幅频特性');subplot(2,3,5);plot(f/pi,angle(h);title('As=30带通滤波器相频特性');As=40;N,wn=buttord(wp,ws,Rp,As);b,a=butter(N,wn);h,f=freqz(b,a,256);subplot(2,3,3);plot(f/pi,abs(h);title('As=40带通滤波器