附录
附录
注意：由于大部分程序都很相似，这里只给出五个主要的源程序。对于那些只需要该参数就可以实现的程序也给只出主要程序。
源程序一：(连续周期信号的分解与综合——谐波分析) function[A_sym,B_sym]=CTFS1 syms t n k x
T=5;tao=0.2*T;a=0.5; if nargin<4;Nf=6;end if nargin<5;Nn=32;end
x=time_fun_x(t); A0=2*int(x,t,-a,T-a)/T; As=int(2*x*cos(2*pi*n*t/T)/T,t,-a,T-a); Bs=int(2*x*sin(2*pi*n*t/T)/T,t,-a,T-a); A_sym(1)=vpa(A0,Nn); for k=1:Nf A_sym(k+1)=vpa(subs(As,n,k),Nn); B_sym(k+1)=vpa(subs(Bs,n,k),Nn); end
if nargout==0
c=A_sym;disp(c) d=B_sym;disp(d) t=-8*a:0.01:T-a;
f1=2*(0.2/2+0.1871.*cos(2*pi*1*t/5)+0.*sin(2*pi*1*t/5)); f2=2*(0.1514.*cos(2*pi*2*t/5)+0.*sin(2*pi*2*t/5)); f3=2*(0.1009.*cos(2*pi*3*t/5)+0.*sin(2*pi*3*t/5)); f4=2*(0.0468.*cos(2*pi*4*t/5)+0.*sin(2*pi*4*t/5)); f5=2*(-0.0312.*cos(2*pi*6*t/5)+0.*sin(2*pi*6*t/5)); f6=f1+f2; f7=f6+f3; f8=f7+f4+f5; subplot(2,2,1) plot(t,f1),hold on title('基波') subplot(2,2,2) plot(t,f6),hold on title('基波+2次谐波') subplot(2,2,3) plot(t,f7),hold on
title('基波+2次谐波+3次谐波') subplot(2,2,4) plot(t,f8),hold on
title('基波+2次谐波+3次谐波+6次谐波')
22
附录
end
function y=time_fun_e a=0.5;T=5;h=1; tao=0.2*T;
t=-8*a:0.01:T-a;
e1=1/2+1/2.*sign(t+tao/2); e2=1/2+1/2.*sign(t-tao/2); y=h.*(e1-e2);
function x=time_fun_x(t) h=1;
x1=sym('Heaviside(t+0.5)')*h; x=x1-sym('Heaviside(t-0.5)')*h;
源程序二：（吉布斯效应） t=-2:0.001:2;
N=input('N='); c0=0.5;
fN=c0*ones(1,length(t)); for n=1:2:N
fN=fN+cos(pi*n*t)*sinc(n/2); end figure plot(t,fN)
axis([-2 2 -0.2 1.2])
源程序三：（矩形脉冲的频谱分析） function [A_sym,B_sym]=CTFS2 syms t n y
if nargin<3;Nf=input('pleas Input 所需展开的最高谐波次数: Nf=');end T=input('pleas Input 信号的周期T='); if nargin<5;Nn=32;end y=fun_in(t);
A0=2*int(y,t,0,T)/T;
As=int(2*y*cos(2*pi*n*t/T)/T,t,0,T); Bs=int(2*y*sin(2*pi*n*t/T)/T,t,0,T); A_sym(1)=double(vpa(A0,Nn)); for k=1:Nf
A_sym(k+1)=double(vpa(subs(As,n,k),Nn)); B_sym(k+1)=double(vpa(subs(Bs,n,k),Nn));end if nargout==0
S1=fliplr(A_sym) S1(1,k+1)=A_sym(1)
23
附录
S2=fliplr(1/2*S1) S3=fliplr(1/2*B_sym) S3(1,k+1)=0 S4=fliplr(S3) S5=S2-i*S4; S6=fliplr(S5); N=Nf*2*pi/T; k2=-N:2*pi/T:N;
S7=[S6,S5(2:end)]; x=fun_mc subplot(3,1,2)
stem(k2,abs(S7)); axis([-150,150,0,0.12]) subplot(3,1,3)
stem(k2,angle(S7)); axis([-150,150,-4,4]) end
function y=fun_in(t) syms a a1
T=input('pleas Input 信号的周期T='); M=input('周期与脉冲宽度之比M='); A=1;tao=T/M;a=tao/2;
y1=sym('Heaviside(t+a1)')*A; y=y1-sym('Heaviside(t-a1)')*A; y=subs(y,a1,a); y=simple(y);
function x=fun_mc T=5;tao=T/5;n=4; t=-n*T:0.01:n*T; x=rectpuls(t,1); for i=1:n;
x=x+rectpuls(t-i*T,1)+rectpuls(t+i*T,1); end
subplot(3,1,1) plot(t,x) hold on
axis([-20,20,0,1.2]) 源程序四：（典型周期脉冲的频谱分析——方波） function [A_sym,B_sym]=fangbo syms t n k y T=10;
if nargin<4;Nf=input('pleas Input 所需展开的最高谐波次数Nf= ');end if nargin<5;Nn=32;end y=fangbo_1;
24
附录
A0=2*int(y,t,0,T)/T;
As=int(2*y*cos(2*pi*n*t/T)/T,t,0,T); Bs=int(2*y*sin(2*pi*n*t/T)/T,t,0,T); A_sym(1)=double(vpa(A0,Nn)); for k=1:Nf
A_sym(k+1)=double(vpa(subs(As,n,k),Nn)); B_sym(k+1)=double(vpa(subs(Bs,n,k),Nn));end if nargout==0
S1=fliplr(A_sym) S1(1,k+1)=A_sym(1) S2=fliplr(1/2*S1) S3=fliplr(1/2*B_sym) S3(1,k+1)=0 S4=fliplr(S3) S5=S2-i*S4; S6=fliplr(S5); N=Nf*2*pi/T;
k2=-N:2*pi/T:N; S7=[S6,S5(2:end)]; subplot(2,1,1)
x=fangbo_2 T=5;t=-2*T:0.01:2*T; plot(t,x)
title('T=5，占空比为50%的周期方波脉冲') axis([-10,10,-1.2,1.2]) subplot(2,1,2)
stem(k2,abs(S7));
title('连续时间函数周期方波脉冲的双边幅度谱') axis([-20,20,0,0.6]) end
function y=fangbo_1 syms a a1 T=5;a=T/2;
y1=sym('Heaviside(t)')*2-sym('Heaviside(t-a1)'); y=y1-sym('Heaviside(t+a1)'); y=subs(y,a1,a); y=simple(y);
function x=fangbo_2
T=5;t=-2*T:0.01:2*T;duty=50; x=square(t,duty); 源程序五：（典型周期脉冲的频谱分析——三角波） function [A_sym,B_sym]=sjb syms t n k y T=10;
25
附录
if nargin<4;Nf=input('pleas Input 所需展开的最高谐波次数Nf= ');end if nargin<5;Nn=32;end y=time_fun_s(t); A0=2*int(y,t,0,T)/T;
As=int(2*y*cos(2*pi*n*t/T)/T,t,0,T); Bs=int(2*y*sin(2*pi*n*t/T)/T,t,0,T); A_sym(1)=double(vpa(A0,Nn)); for k=1:Nf
A_sym(k+1)=double(vpa(subs(As,n,k),Nn)); B_sym(k+1)=double(vpa(subs(Bs,n,k),Nn));end if nargout==0
S1=fliplr(A_sym) S1(1,k+1)=A_sym(1) S2=fliplr(1/2*S1) S3=fliplr(1/2*B_sym) S3(1,k+1)=0 S4=fliplr(S3) S5=S2-i*S4; S6=fliplr(S5); N=Nf*2*pi/T;
k2=-N:2*pi/T:N; S7=[S6,S5(2:end)]; subplot(2,1,1)
x=sjb_timefun T=10;t=-2*T:0.01:2*T; plot(t,x)
title('T=10，脉宽为1的周期三角波脉冲') axis([-20,20,-1,1.2]) subplot(2,1,2)
stem(k2,abs(S7));
title('连续时间函数周期三角脉冲的双边幅度谱') axis([0,20,0,0.6]) end
function x=sjb_timefun T=10;t=-2*T:0.01:2*T; x=sawtooth(t-2*T/3,1);
function y=time_fun_s(t) syms a a1 T=10;a=T/2;
y=sym('Heaviside(t+a1)')*(2*t/a1)-sym('Heaviside(t-a1)')*(2*t/a1); %y=y1-sym('Heaviside(t)')*(4*t/a1); y=subs(y,a1,a); y=simple(y);
26