测试信号处理（第3版）程序代码 第3章

例3.1%例3.1求解周期信号傅里叶级数的MATLAB程序%先求出信号傅里叶级数的系数symsTEtnx;pi=sym('pi');%创建符号对象a0=2/T*int(-E,t,-T/2,0)+2/T*int(E,t,0,T/2),A0=a0/2;%计算系数an=2/T*int(-E*cos(2*pi*n*t/T),t,-T/2,0)+2/T*int(E*cos(2*pi*n*t/T),t,0,T/2)bn=2/T*int(-E*sin(2*pi*n*t/T),t,-T/2,0)+2/T*int(E*sin(2*pi*n*t/T),t,0,T/2)%求出E=1时各谐波分量的幅度symsEAncbn=2*E*c;E=1;c=-(cos(n*pi)-1)/(n*pi);bn=subs(bn,{sym('E'),sym('c')},{E,c})%符号变量置换%画出频谱图A=[00000000000];%存储11个谐波分量幅度的数组n=1:1:11;forn=1:11bn=2*(-cos(n*pi)+1)/(n*pi);A(n)=double(vpa(bn));%对任意精度的符号类数据进行规范An=A(n)endx1=1;x2=11;x=x1:x2;stem(x,A,'r','filled')%画出11个谐波分量的幅度谱例3.3%产生周期为2π的方波信号t=0:1/10:10;y=square(2*pi*(0.05*pi)*t);subplot(2,2,1);plot(t,y);xlabel('(a)方波波形')axis([010-1.51.5])%基波信号t=0:1/10:10;y=sin(t);subplot(2,2,2);plot(t,y);xlabel('(b)基波波形')%基波+3次谐波合成的波形t=0:1/10:10;y=sin(t)+sin(3*t)/3;subplot(2,2,3);plot(t,y);xlabel('(c)基波+3次谐波')%基波+3次+5次+7次+9次谐波合成的波形t=0:1/10:10;y=sin(t)+sin(3*t)/3+sin(5*t)/5+sin(7*t)/7+sin(9*t)/9;subplot(2,2,4);plot(t,y);xlabel('(d)基波+3次+5次+7次+9次谐波') %最高谐波阶次为19时的吉伯斯现象t=0:31/1000:5;y=zeros(10,max(size(t)));x=zeros(size(t));fork=1:2:19,x=x+sin(k*t)/k;y((k+1)/2,:)=x;end;pause,plot(y(1:9,:)),pause,mesh(y),pauseclc%最高谐波阶次为45时的吉伯斯现象y=zeros(10,max(size(t)));x=zeros(size(t));fork=1:2:45,x=x+sin(k*t)/k;y((k+1)/2,:)=x;end;pause,plot(y(1:9,:)),pause,mesh(y),pauseclc例3.10%例3.10中傅里叶变换的MATLAB程序symsEtaw;f=E*Heaviside(t+a)-E*Heaviside(t-a)F=fourier(f);F1=simplify(F)%得到F的简化形式E=1;a=1;F2=subs(F1,{sym('E'),sym('a')},{Ea})%将F1中的变量替换为新变量ezplot(abs(F2),[-3*pi,3*pi])title('F(Ω)')xlabel('Ω')f2=ifourier(F2)运行结果为F=(E*(sin(a*w)+cos(a*w)*1i))/w-(E*(cos(a*w)*1i-sin(a*w)))/w例3.14%例3.14高斯信号的频谱symstwsigma=sym('sigma','positive');f=exp(-t^2/(2*sigma^2))/(sqrt(2*pi)*sigma);F=fourier(f,t,w);F1=simplify(F);digits(4);F2=vpa(F1)sigma=1;F3=subs(F2,'sigma',sigma)ezplot(F3,[-2*pi,2*pi]),gridontitle('F(Ω)');xlabel('Ω')例3.15%信号最高频率wm＝2ws（ws:采样频率）临界采样的情况wm=1;%信号最高频率（带宽）wc=wm;%预处理用低通滤波器截止频率Ts=pi/wm;%采样间隔ws=2*pi/Ts;%变为采样角频率n=-100:100;%采样点数nTs=n*Ts;f=sinc(nTs/pi);%对信号Sa(t)的抽样信号Dt=0.005;t=-10:Dt:10;fa=f*Ts*wc/pi*sinc((wc/pi)*(ones(length(nTs),1)*t-nTs'*ones(1,length(t))));%信号的重构t1=-10:0.5:10;f1=sinc(t1/pi);%Sa(t)的连续波形表示%画出抽样信号与抽样信号的重构波形subplot(2,1,1);stem(t1,f1);xlabel('kTs');ylabel('f(kTs)');title('sa(t)的临界抽样');subplot(2,1,2);plot(t,fa)xlabel('t');ylabel('fa(t)');title('由抽样信号重构sa(t)');grid%欠抽样情况，抽样频率ws与信号最高频率wm之间，wm>2ws.wm=1;wc=1.1*wm;Ts=2*pi/wm;ws=2*pi/Ts;n=-100:100;nTs=n*Ts;f=sinc(nTs/pi);Dt=0.005;t=-10:Dt:10;fa=f*Ts*wc/pi*sinc((wc/pi)*(ones(length(nTs),1)*t-nTs'*ones(1,length(t))));error=abs(fa-sinc(t/pi));%求解重构信号与原信号的误差t1=-10:0.5:10;f1=sinc(t1/pi);subplot(311);stem(t1,f1);xlabel('kTs');ylabel('f(kTs)');title('sa(t)的欠抽样信号');subplot(312);plot(t,fa)xlabel('t');ylabel('fa(t)');title('由欠抽样信号重构的sa(t)');grid;subplot(313);plot(t,error);xlabel('t');ylabel('error(t)');title('欠抽样信号与原信号的误差error(t)')%过抽样情况，wm<2ws.wm=1;wc=1.1*wm;Ts=0.5*pi/wm;ws=2*pi/Ts;n=-100:100;nTs=n*Ts;f=sinc(nTs/pi);Dt=0.005;t=-10:Dt:10;fa=f*Ts*wc/pi*sinc((wc/pi)*(ones(length(nTs),1)*t-nTs'*ones(1,length(t))));error=abs(fa-sinc(t/pi));t1=-10:0.5:10;f1=sinc(t1/pi);subplot(311);stem(t1,f1);xlabel('kTs');ylabel('f(kTs)')