计算机编程与数字信号处理实习

1、计算机编程与数字信号处理实习能够利用Matlab熟悉地画图，内容包括：X、Y坐标轴上的label，每幅图上的title，绘画多条曲线时的legend，对图形进行适当的标注等。 （1）在一副图上画出多幅小图；clc;clear;close;x=linspace(0,5,200);y1=sin(3*x);y2=sin(x.3);y3=exp(x-2)y4=sin(2*x)+cos(3*x);for k=1:4; if k=1 subplot(2,2,k) plot(x,y1) hold on; title('x与y1') elseif k=2 subplot(2,2,k) plot

2、(x,y2) hold on; title('x与y2') elseif k=3 subplot(2,2,k) plot(x,y3) hold on; title('x与y3') elseif k=4 subplot(2,2,k) plot(x,y4) hold on; title('x与y4')endendprint -djpeg -r0 xinhao1.jpeg （2）画出一组二维图形；（以分段函数为例）clc;clear;close;x=-3:0.01:3; y1=zeros(size(x); y2=zeros(size(x); y3=ze

3、ros(size(x); N=length(x); for k=1:N if x(k)<-1&x(k)>=-3; y1(k)=(-x(k).2-4*x(k)-3)/2; elseif x(k)>=-1&x(k)<1 ; y2(k)=-x(k).2+1; else x(k)<=3&x(k)>=1 ; y3(k)=(-x(k).2+4*x(k)-3)/2; end end y=y1+y2+y3; plot(x,y)xlabel('x');ylabel('y');print -djpeg -r0 xinhao

4、2.jpeg （3）画出一组三维图形；clearx=0:0.2:3*pi;y=0:0.25:5*pi;xx,yy=meshgrid(x,y);z1=sin(xx).*sin(yy);x=-3:0.25:3;y=x;xx,yy=meshgrid(x,y);z2=xx-0.5*xx.3+0.2*yy.2+1;x=-8:0.5:8;y=x;xx,yy=meshgrid(x,y);r=sqrt(xx.2+yy.2)+eps;z3=sin(r)./r;subplot(2,2,1),mesh(z1);title('z1')subplot(2,2,2),mesh(z2);title('

5、;z2')subplot(2,2,3),waterfall(z2);title('z2')subplot(2,2,4),mesh(z3);title('z3')print -djpeg -r0 xinhao3.jpeg说明：第一个图和第二个图中对x1赋值方式不同，结果相同。可以看出，无论是第一段中对X1的定义，还是第二段中对X1的定义。只要X1的维度为10，就符合surf的输入条件。且第一段中X1的范围就是常数9，X2中同样不变，所以画出的图像当然是一样的。 （4）画出复数的实部与虚部。clc;clear;close;x=20;y=rand(x,1)+r

6、and(x,1)*i;m=real(y);n=imag(y); subplot(2,2,1)plot(y,'.');title('复数');hold on;grid on;subplot(2,2,2)plot(m,n,'.');title('复数（实部为横，虚部为纵)');xlabel('实部');ylabel('虚部');hold on;grid on;subplot(2,2,3)plot(m,'.');title('复数实部');hold on;grid on;s

7、ubplot(2,2,4)plot(n,'.');title('复数虚部');hold on;grid on;print -djpeg -r0 xinhao4.jpeg （5）完成对一个源程序进行详细注释。（仍以分段函数为例）clc;clear x=-3:0.01:3; %从-3到+3以0.01为步长给x赋值y1=zeros(size(x); %初始化y1的值y2=zeros(size(x); %初始化y1的值y3=zeros(size(x); %初始化y1的值N=length(x); %以x的维度给N赋值for k=1:N %k从1循环到Nif x(k)<

8、;-1&x(k)>=-3; %分段赋值y1(k)=(-x(k).2-4*x(k)-3)/2; elseif x(k)>=-1&x(k)<1 ; y2(k)=-x(k).2+1; else x(k)<=3&x(k)>=1 ; y3(k)=(-x(k).2+4*x(k)-3)/2; end end y=y1+y2+y3; %各段叠加成分段函数yplot(x,y)%以x为横坐标，y为纵坐标作图二、计算普通褶积与循环褶积，分别使用时间域与频率域两种方法进行正、反演计算，指出循环褶积计算时所存在的边界效应现象；编写一个做相关分析的源程序。（1）作线性

9、褶积，时间域正演，频率域反演（无）clc;clear;close;x=0 1 2'y=0 1 2 3'xy0=conv(x,y);%利用公式计算线性褶积yx0=conv(y,x);%利用公式计算线性褶积%编程计算%N1 =length(x); %N2 =length(y); %NN =N1+N2-1; %xy =zeros(1,NN); %yx =zeros(1,NN); %A1 =zeros(NN,N2); %A2 =zeros(NN,N1); %for j=1:N2 % for i=j:j+N1-1 % A1(i,j)=x(i-j+1); % end %end %for j

10、=1:N1 % for i=j:j+N2-1 % A2(i,j)=y(i-j+1); % end %end %xy1=A1*y; %yx1=A2*x; %结果验证%xy0yx0xy1yx1max(abs(xy0-yx0)%说明褶积的可交换性max(abs(xy1-yx1)%说明褶积的可交换性max(abs(xy0-xy1)%说明公式计算和编程计算的结果一致性%作图%subplot(2,2,1)stem(x)grid on;hold on;subplot(2,2,2)stem(y)grid on;hold on;subplot(2,2,3)stem(xy0)grid on;hold on;sub

11、plot(2,2,4)stem(yx0)grid on;hold on;print -djpeg -r0 xinhao.jpeg运行结果：xy0 = 0 0 1 4 7 6yx0 = 0 0 1 4 7 6xy1 = 0 0 1 4 7 6yx1 = 0 0 1 4 7 6ans = 0ans = 0ans = 0（2）作循环褶积，时间域正演与频率域反演clc;clear;close;x=0 1 2'y=0 1 2 3' xy0=conv(x,y);%利用公式计算线性褶积yx0=conv(y,x);%利用公式计算线性褶积%自带公式频率域反演计算%N1 =length(x); %

12、N2 =length(y); %N =N1+N2-1; %NN =N1+N2-1; %w =exp(-sqrt(-1)*2*pi/N); %x0=zeros(N,1); % 需要将原始数组转化为相同长度的两个数组（补零）。y0=zeros(N,1); % 需要将原始数组转化为相同长度的两个数组（补零）。x0(1:N1)=x; %y0(1:N2)=y; %X0=fft(x0); %Y0=fft(y0); %xy2=ifft(X0.*Y0); %yx2=ifft(Y0.*X0); %频率域反演%F=zeros(N); %for i=1:N % for j=1:N % F(i,j)=power(w,

13、(i-1)*(j-1); % end %end %X1=F*x0; %Y1=F*y0; %xy3=1/N*conj(F)*(X1.*Y1); % inv(F)=conj(F)/Nyx3=1/N*conj(F)*(Y1.*X1); %时间域直接计算的方法%C1 =zeros(N); %C2 =zeros(N); %for j=1:N % for i=1:N % k=i+j-1; % if k>N % k=k-N; % end % C1(k,j)=x0(i); % C2(k,j)=y0(i); % end %end %xy4=C1*y0; %yx4=C2*x0; %结果验证%xy2yx2xy

14、3yx3xy4yx4max(abs(xy2-yx2)max(abs(xy3-yx3)max(abs(xy4-yx4)max(abs(xy2-xy3)max(abs(xy2-xy4)max(abs(xy2(1:NN)-xy0) % 验证结果：可以利用循环褶积的方法计算线性褶积。 %作图%subplot(2,2,1)stem(x)grid on;hold on;subplot(2,2,2)stem(y)grid on;hold on;subplot(2,2,3)stem(xy2)grid on;hold on;subplot(2,2,4)stem(yx2)grid on;hold on;print

15、 -djpeg -r0 xinhao.jpeg运行结果：xy2 = 0 0 1 4 7 6yx2 = 0 0 1 4 7 6xy3 = 0.0000 - 0.0000i -0.0000 - 0.0000i 1.0000 - 0.0000i 4.0000 + 0.0000i 7.0000 + 0.0000i 6.0000 - 0.0000iyx3 = 0.0000 - 0.0000i -0.0000 - 0.0000i 1.0000 - 0.0000i 4.0000 + 0.0000i 7.0000 + 0.0000i 6.0000 - 0.0000ixy4 = 0 0 1 4 7 6yx4 =

16、 0 0 1 4 7 6ans = 0ans = 0ans = 0ans = 7.8221e-015ans = 0ans = 0这些ans=0说明了各种方法计算得到的褶积结果是相等的，并且说明了褶积的可交换性。（3）循环褶积计算时所存在的边界效应现象变量赋值如上所示，程序不再累述分别以xy1，xy4代表线性褶积和循环褶积当N=4时结果如下：xy1 = 0 0 1 4 7 6xy4 = 7 6 1 4当N=5时结果如下：xy1 = 0 0 1 4 7 6xy4 = 6 0 1 4 7当N=6时结果如下：xy1 = 0 0 1 4 7 6xy4 = 0 0 1 4 7 6当N=7时结果如下：xy1

17、 = 0 0 1 4 7 6xy4 = 0 0 1 4 7 6 0由此可见，当NN1+N2-1时，可以利用循环褶积计算线性褶积，N1+N2-1可以看做一个边界值，这种现象称之为循环褶积计算时所存在的边界效应现象。（4）相关分析的示例clc;clear;close;x=0 1 2'y=0 1 2 3' N1 =length(x); %N2 =length(y); %NN =N1+N2-1; % y0=zeros(NN,1);%构造相关乘积用的列向量y0(N1:NN)=y;%把y的值赋给y0A1 =zeros(NN,N2); %线性褶积矩阵A2 =zeros(NN,N1); %线性

18、相关矩阵for j=1:N2 % for i=j:j+N1-1 % A1(i,j)=x(i-j+1); % end %end %for j=1:NN for i=NN-j+1:min(NN-j+1+2,NN) % A2(i,j)=x(i-(NN-j); % end end %xyZhe=A1*y; %线性褶积xyXiang=A2*y0; %线性相关xyXiangN=xcorr(x,y) %由matlab自带的函数计算相关%作图%subplot(2,2,1)stem(x,'r','.');xlabel('n');ylabel('x'

19、);title('x向量');grid on;subplot(2,2,2)stem(y,'r','.');xlabel('n');ylabel('y');title('y向量');grid on;subplot(2,2,3)stem(xyXiang,'r','.');xlabel('n');ylabel('xyXiang');title('自己计算的线性相关');grid on;subplot(2,2,4)stem(xy

20、XiangN,'r','.');xlabel('n');ylabel('xyXiangN');title('xcorr计算得到的相关');grid on;print -djpeg -r0 xiangguang.jpeg 运行结果：由此可见，编程用循环得到的相关结果和Matlab提供的xcorr函数计算的结果一致，只不过，xcorr计算得到的结果是奇数项表示的。三、设计一个病态（矩阵）系统，分析其病态程度；找出对应的解决方法（提示：添加白噪因子）。Part One病态矩阵和及其病态程度的判断：组别的变量有冗余的矩阵，

21、或者说特征值有0的矩阵就是病态矩阵。但是在实际工作中，判断一个矩阵是否为病态矩阵是相对的。通常通过用小扰动后解的改变程度来判断。因此定义病态矩阵是求解方程组时对数据的小扰动很敏感的矩阵。解线性方程组Ax=b时，若对于系数矩阵A及右端项b的小扰动，方程组的解x与原方程组Ax=b的解差别很大，则称矩阵A为病态矩阵。方程组的近似解x一般都不可能恰好使剩余r=b-Ax为零，这时x亦可看作小扰动问题Ax=b-r(即)的解，所以当A为病态时，即使剩余r很小，仍可能得到一个与真解相差很大的近似解。例如，取 如果取作为近似解，剩余已经很小，但x与真解相差还是很大。用来衡量矩阵的病态程度，K(A)称为A的条件数

22、，它很大时，称A为病态，否则称良态；K(A)愈大，A的病态程度就愈严重。 计算得，这说明的病态程度很大。我们对，都做一个小小的扰动，那么的解与原来的真解差别很大，由此说明是病态矩阵。如果用单位矩阵来做比较的话：它的解就等于我们再对，做一个小小的扰动，那么的解与真解完全相等，这说明单位矩阵不是病态的，再用cond函数计算得，结论就是E不是病态矩阵。因此解决病态问题的关键是在其扰动的条件下，仍然保持解的不变，这就需要添加白噪因子（可以看做第一次扰动）。Part Two解决病态矩阵方程组AA*X1=C0问题clc;clear;close;A=1 2 2 1'D=1 3 5 7'B=1

23、 2 3 4'N=length(A);F=zeros(N); % w(i*j) AA=zeros(N);DD=zeros(N);w=exp(-sqrt(-1)*2*pi/N);for i=1:N for j=1:N AA(i,j)=A(mod(i-j,N)+1); % <i-j>N DD(i,j)=D(mod(i-j,N)+1); F(i,j)=power(w,(i-1)*(j-1); % DFT. % w(i-1)*(j-1) endend FA=F*A;FB=F*B; FAA=F*AA; %定义FAAFDD=F*DD; %定义FDD C0=AA*B; % = ifft(

24、fft(A).*fft(B);%C1=real(conj(F)/N*(FA.*FB); %inv(F)=conj(F)/N; %计算病态矩阵方程式AA*X1=C0%以及良态矩阵方程式DD*X2=C0%AA*X1=C0%F*AA*inv(F)*F*X1=F*C0%F*AA*(1/N)*conj(F)*F*X1=F*C0%(1/N)*conj(F)*AA*F*AA*conj(F)*F*X1=conj(F)*AA*F*C0%(1/N)*C_FA*F*AA*conj(F)*F*X1=C_FA*FC%(1/N)*C_FA*FAA*conj(F)*F*X1=FC2%(1/N)*(C_FA*FAA+0.00

25、1*eye(N)*conj(F)*F*X1=FC2%添加白噪因子%(1/N)*FA2*conj(F)*F*X1=FC2%FA2*X1=FC2%X=inv(FA2)*FC2FC=F*C0;C_FA=conj(FAA);FC2=C_FA*FC;FA2=C_FA*FAA+0.001*eye(N);X1=inv(FA2)*FC2 REALX1=inv(AA)*C0%BX=inv(AA+0.001*eye(N)*C0; % %DD*X2=C0%F*DD*inv(F)*F*X2=F*C0%F*DD*(1/N)*conj(F)*F*X2=F*C0%(1/N)*conj(F)*DD*F*DD*conj(F)*

26、F*X2=conj(F)*DD*F*C0%(1/N)*C_FD*F*DD*conj(F)*F*X2=C_FD*FC%(1/N)*C_FD*FDD*conj(F)*F*X2=FC2%(1/N)*(C_FD*FDD+0.001*eye(N)*conj(F)*F*X2=FC2%添加白噪因子%(1/N)*FD2*conj(F)*F*X2=FC2%FD2*X2=FC2%X=inv(FD2)*FC2FC=F*C0;C_FD=conj(FDD);FC2=C_FD*FC;FD2=C_FD*FDD+0.001*eye(N);X2=inv(FD2)*FC2 REALX2=inv(DD)*C0运行结果：X1 = -

27、1.9995 + 0.0006i 5.0004 - 0.0008i 0.0000 + 0.0000i 6.9996 + 0.0001iREALX1 = Inf Inf Inf InfX2 = 0.6875 + 0.0000i 1.1875 + 0.0000i 1.1875 - 0.0000i 0.6875 - 0.0000iREALX2 = 0.6875 1.1875 1.1875 0.6875结果分析：病态矩阵AA添加白噪因子前解为Inf，添加白噪因子后方可求解。良态矩阵DD添加白噪因子前解与添加白噪因子后的解是相等的。四、设计一个一维滤波处理程序（1、分别做低通、高通、带通、带阻等理想滤波

28、器进行处理；2、窗函数法或者镶边法设计滤波器，以消除Gibbs现象）。直接设计的理想滤波器（使用时会出现Gibbs现象）：clc;clear;close;f1=15,f2=35;%设定滤波器的两个频率d=0.01;%采样间隔t=0:d:1;%建立时间向量N=length(t);%时间向量长度x=5*rand(size(t)-2;%随机信号，含中高低频fc=1/(2*d);%截止频率f=0:(1/d)/N:(1/d);%求出频谱对应的频率序列FX=fft(x);%求原始信号的频谱 %理想低通滤波器设计H1=zeros(1,N);for i=1:N if f(i)<=f1|f(i)>=

29、(1/d)-f1 H1(i)=1; else H1(i)=0; endendFX1=FX.*H1; %频率滤波X1=ifft(FX1); %频率滤波后的反演信号 %理想高通滤波器设计H2=zeros(1,N);for i=1:N if (f(i)>=f2)&(f(i)<(1/d)-f2) H2(i)=1; else H2(i)=0; endendFX2=FX.*H2;X2=ifft(FX2); %理想带通滤波器设计H3=zeros(1,N);for i=1:N if (f(i)<=f2&f(i)>=f1)|(f(i)>2*fc-f2&f(i

30、)<(1/d)-f1); H3(i)=1; else H3(i)=0; endendFX3=FX.*H3;X3=ifft(FX3); %理想带阻滤波器设计H4=zeros(1,N);for i=1:N if (f(i)<=f2&f(i)>=f1)|(f(i)>2*fc-f2&f(i)<2*fc-f1); H4(i)=0; else H4(i)=1; endendFX4=FX.*H4;X4=ifft(FX4);%作图% %画出各滤波器的频谱figure(1);subplot(2,2,1),plot(H1);grid on;axis(0 100 -0.

31、5 1.5);title('理想低通H1(f)');subplot(2,2,2),plot(H2);grid on;axis(0 100 -0.5 1.5);title('理想高通H2(f)');subplot(2,2,3),plot(H3);grid on;axis(0 100 -0.5 1.5);title('理想带通H3(f)');subplot(2,2,4),plot(H4);grid on;axis(0 100 -0.5 1.5);title('理想带阻H4(f)');print -djpeg -r0 滤波器频谱图.jp

32、eg %低通滤波figure(2);subplot(2,2,1),plot(t,x);grid on;axis(0 1 -4 4);title('原始信号');subplot(2,2,2),stem(abs(FX),'.');grid on;axis(0 100 0 60);title('原始信号振幅谱');subplot(2,2,3),plot(t,X1);grid on;axis(0 1 -4 4);title('低通滤波后的低频信号');subplot(2,2,4),stem(abs(FX1),'.');gr

33、id on;axis(0 100 0 60);title('低通滤波后的低频信号的振幅谱');print -djpeg -r0 低通滤波.jpeg %高通滤波figure(3);subplot(2,2,1),plot(t,x);grid on;axis(0 1 -4 4);title('原始信号');subplot(2,2,2),stem(abs(FX),'.');grid on;axis(0 100 0 60);title('原始信号振幅谱');subplot(2,2,3),plot(t,X2);grid on;axis(0 1

34、 -4 4);title('高通滤波后的高频信号');subplot(2,2,4),stem(abs(FX2),'.');grid on;axis(0 100 0 60);title('高通滤波后的高频信号的振幅谱');print -djpeg -r0 高通滤波.jpeg %带通滤波figure(4);subplot(2,2,1),plot(t,x);grid on;axis(0 1 -4 4);title('原始信号');subplot(2,2,2),stem(abs(FX),'.');grid on;axis(

35、0 100 0 60);title('原始信号振幅谱');subplot(2,2,3),plot(t,X3);grid on;axis(0 1 -4 4);title('带通滤波后的中频信号');subplot(2,2,4),stem(abs(FX3),'.');grid on;axis(0 100 0 60);title('带通滤波后的中频信号的振幅谱');print -djpeg -r0 带通滤波.jpeg %带阻滤波figure(5);subplot(2,2,1),plot(t,x);grid on;axis(0 1 -4

36、4);title('原始信号');subplot(2,2,2),stem(abs(FX),'.');grid on;axis(0 100 0 60);title('原始信号振幅谱');subplot(2,2,3),plot(t,X4);grid on;axis(0 1 -4 4);title('带阻滤波后的含有高低频信号');subplot(2,2,4),stem(abs(FX4),'.');grid on;axis(0 100 0 60);title('带阻滤波后的高低频信号的振幅谱');prin

37、t -djpeg -r0 带阻滤波.jpeg各滤波器频谱图理想高通滤波理想低通滤波理想带通滤波理想带阻滤波镶边法设计的滤波器：clc;clear;close;f1=15,f2=35;%设定滤波器的两个频率，f1为低通滤波器的高通频率、带通滤波器的低通频率，f2为高通滤波器的高截频率、带通滤波器的高通频率d=0.01;%采样间隔t=0:d:1;%建立时间向量N=length(t);%时间向量长度%x=sin(2*pi*5*t)+sin(2*pi*20*t)+sin(2*pi*40*t);%由时间向量构成的含有低中高频的原始信号x=5*rand(size(t)-2;%随机信号，含中高低频fc=1/

38、(2*d);%截止频率f=0:(1/d)/N:(1/d);%求出频谱对应的频率序列FX=fft(x);%求原始信号的频谱%镶边函数%HH1=zeros(1,N);for i=1:N if (f(i)<=45)&(f(i)>=15) HH1(i)=0.5*sin(1/20)*pi*i+5)+0.5; else HH1(i)=1; endendHH2=zeros(1,N);for i=1:N if (f(i)<=85)&(f(i)>=55) HH2(i)=0.5*sin(1/20)*pi*i+7)+0.5; else HH2(i)=1; endendHH3=

39、zeros(1,N);for i=1:N if (f(i)<=45)&(f(i)>=30) HH3(i)=0.5*sin(1/15)*pi*i+11)+0.5; else HH3(i)=1; endendHH4=zeros(1,N);for i=1:N if (f(i)<=70)&(f(i)>=55) HH4(i)=0.5*sin(1/15)*pi*i+9)+0.5; else HH4(i)=1; endendHH5=zeros(1,N);for i=1:N if (f(i)<=22)&(f(i)>=10) HH5(i)=0.5*si

40、n(1/10)*pi*i+7)+0.5; else HH5(i)=1; endendHH6=zeros(1,N);for i=1:N if (f(i)<=45)&(f(i)>=35) HH6(i)=0.5*sin(1/10)*pi*i+15)+0.5; else HH6(i)=1; endendHH7=zeros(1,N);for i=1:N if (f(i)<=65)&(f(i)>=55) HH7(i)=0.5*sin(1/10)*pi*i+6)+0.5; else HH7(i)=1; endendHH8=zeros(1,N);for i=1:N if

41、 (f(i)<=90)&(f(i)>=78) HH8(i)=0.5*sin(1/10)*pi*i+2)+0.5; else HH8(i)=1; endendHH9=zeros(1,N);for i=1:N if (f(i)<=20)&(f(i)>=10) HH9(i)=0.5*sin(1/10)*pi*i+11)+0.5; else HH9(i)=1; endendHH10=zeros(1,N);for i=1:N if (f(i)<=40)&(f(i)>=30) HH10(i)=0.5*sin(1/10)*pi*i+1)+0.5;

42、else HH10(i)=1; endendHH11=zeros(1,N);for i=1:N if (f(i)<=70)&(f(i)>=58) HH11(i)=0.5*sin(1/10)*pi*i+8)+0.5; else HH11(i)=1; endendHH12=zeros(1,N);for i=1:N if (f(i)<=88)&(f(i)>=70) HH12(i)=0.5*sin(1/10)*pi*i+5)+0.5; else HH12(i)=1; endend %低通滤波器设计H1=zeros(1,N);for i=1:N if f(i)&l

43、t;=36|f(i)>=64 H1(i)=1; else H1(i)=0; endendH1=HH1.*H1;%镶边H1=HH2.*H1;%镶边FX1=FX.*H1;%频率滤波X1=ifft(FX1);%频率滤波后的反演信号 %高通滤波器设计H2=zeros(1,N);for i=1:N if (f(i)>=30)&(f(i)<=70) H2(i)=1; else H2(i)=0; endendH2=HH3.*H2;%镶边H2=HH4.*H2;%镶边FX2=FX.*H2;%频率滤波X2=ifft(FX2);%频率滤波后的反演信号 %带通滤波器设计H3=zeros(1,

44、N);for i=1:N if (f(i)<=45&f(i)>=10)|(f(i)>55&f(i)<88); H3(i)=1; else H3(i)=0; endendH3=HH5.*H3;%镶边H3=HH6.*H3;%镶边H3=HH7.*H3;%镶边H3=HH8.*H3;%镶边FX3=FX.*H3;X3=ifft(FX3); %带阻滤波器设计H4=zeros(1,N);for i=1:N if (f(i)<30&f(i)>20)|(f(i)>69&f(i)<80); H4(i)=0; else H4(i)=1;

45、endendH4=HH9.*H4;%镶边H4=HH10.*H4;%镶边H4=HH11.*H4;%镶边H4=HH12.*H4;%镶边FX4=FX.*H4;X4=ifft(FX4);%作图% %画出各滤波器的频谱figure(1);subplot(2,2,1),plot(H1);grid on;axis(0 100 -0.5 1.5);title('镶边低通H1(f)');%hold on;plot(HH2,'r');subplot(2,2,2),plot(H2);grid on;axis(0 100 -0.5 1.5);title('镶边高通H2(f)&#

46、39;);%hold on;plot(HH4);subplot(2,2,3),plot(H3);grid on;axis(0 100 -0.5 1.5);title('镶边带通H3(f)');%hold on;plot(HH8);subplot(2,2,4),plot(H4);grid on;axis(0 100 -0.5 1.5);title('镶边带阻H4(f)');%hold on;plot(HH12,'r');print -djpeg -r0 滤波器频谱图.jpeg %低通滤波figure(2);subplot(2,2,1),plot(t,x);grid on;axis(0 1 -4 4);title('原始信号');subplot(2,2,2),stem(abs(FX),'.');grid on;axis(0 100 0 60);title('原始信号振幅谱');subplot(2,2,3),plot(t,X1);grid on;axis(0 1 -4 4);title('低