实验三：周期信号的傅利叶级数分析及MATLAB实现
学院：信息学院 专业：通信工程(武警国防生) 指导教师：肖琦 姓名：梁由勇 学号：20081060094 成绩：
实验目的：1.用MATLAB实现周期信号的傅利叶级数；
2.用MATLAB实现周期信号(包括典型周期信号)的频谱分析；
3.观察利用MATLAB生成的图形及结果，与信号与系统理论知识相连系，加深对信号与系统理论知识的深入理解。
9.1 已知周期半波余弦信号和周期全波余弦信号的波形分别如图所示（图略），用MATLAB编程
求出它们的傅里叶系数，绘出其直流、一次、二次、三次、四次及五次谐波叠加后的波形图，并将其与原周期信号的时域波形进行比较，观察周期信号的分解与合成过程。 解：
% 观察周期余弦半波信号的分解和合成 % m:傅立叶级数展开的项数
display('Please input the value of m (傅立叶级数展开的项数)'); m=input('m = '); t=-2.5*pi:0.01:2.5*pi; t1=-0.5*pi:0.01:0.5*pi;
n=round(length(t)/5);
f=[cos(t1)';zeros(n-1,1);cos(t1)';zeros(n-1,1);cos(t1)']; y=zeros(m+1,max(size(t))); y(m+1,:)=f';
subplot((m+2),1,1)
plot(t/pi,y(m+1,:)); grid;
axis([-2.5 2.5 -0.5 1.5]); title('周期半波余弦信号'); xlabel('t/pi','Fontsize', 8); x=zeros(size(t)); kk='1';
%计算系数
syms tx n T=2*pi;
fx=sym('cos(tx)'); Nn=30;
an=zeros(m+1,1); bn=zeros(m+1,1);
A0=2*int(fx,tx,-T/4,T/4)/T;
An=2*int(fx*cos(2*pi*(n+eps/2)*tx/T),tx,-T/4,T/4)/T; Bn=2*int(fx*sin(2*pi*(n+eps/2)*tx/T),tx,-T/4,T/4)/T; an(1)=double(vpa(A0,Nn)); an(2)=0.5; for k=2:m
an(k+1)=double(vpa(subs(An,n,k),Nn)); bn(k+1)=double(vpa(subs(Bn,n,k),Nn)); end
%计算直流分量
pause;
x=an(1)*cos(0*t)/2; plot(t/pi,y(m+1,:)); hold on;
信号与线性系统分析实验
plot(t/pi,x); grid;
hold off;
axis([-2.5 2.5 -0.5 1.5]); title('直流分量');
xlabel('t/pi','Fontsize', 8);
%各次谐波叠加
for k=1:m pause;
x=x+an(k+1).*cos(k*t); y(k,:)=x;
subplot((m+2),1,k+1); plot(t/pi,y(m+1,:)); hold on;
plot(t/pi,y(k,:)); hold off; grid;
axis([-2.5 2.5 -0.5 1.5]); title(strcat('第',kk,'次谐波叠加')); xlabel('t/pi','Fontsize', 8); kk=strcat(kk,'、',num2str(k+1)); end pause;
subplot((m+2),1,m+2) plot(t/pi,y(1:m+1,:)); grid;
axis([-2.5 2.5 -0.5 1.5]); title('各次谐波叠加波形');
xlabel('t/pi','Fontsize', 8); % End
% 观察周期余弦全波信号的分解和合成 % m:傅立叶级数展开的项数
display('Please input the value of m (傅立叶级数展开的项数); t = -2.5*pi:0.01:2.5*pi; t1=-0.5*pi:0.01:0.5*pi-0.01;
n = round(length(t)/5);
f = [cos(t1)';cos(t1)';cos(t1)';cos(t1)';cos(t1)';0]; y = zeros(m+1,max(size(t))); y(m+1,:) = f'; subplot(m+2,1,1)
plot(t/pi,y(m+1,:)); grid on; axis([-2.5 2.5 -0.2 1.2]); title('周期全波余弦信号'); xlabel('t/pi','Fontsize', 8); x=zeros(size(t)); kk = '1';
%计算系数
2
信号与线性系统分析实验
syms tx n T=pi;
fx=sym('cos(tx)'); Nn=32;
an=zeros(m+1,1); bn=zeros(m+1,1);
A0 =2*int(fx,tx,-T/2,T/2)/T;
An=2*int(fx*cos(2*pi*(n+eps/2)*tx/T),tx,-T/2,T/2)/T; Bn=2*int(fx*sin(2*pi*(n+eps/2)*tx/T),tx,-T/2,T/2)/T; an(1) = double(vpa(A0,Nn)); for k=1:m
an(k+1)=double(vpa(subs(An,n,k),Nn)); bn(k+1)=double(vpa(subs(Bn,n,k),Nn)); end
%求直流信号 pause;
x=an(1)*cos(0*t)/2; subplot(m+2,1,1)
plot(t/pi,y(m+1,:)); hold on;
plot(t/pi,x); grid on; hold off;
axis([-2.5 2.5 -0.2 1.2]); title('周期全波余弦信号');
xlabel('t/pi','Fontsize', 8);
%各次谐波叠加 for k=1:m pause;
x=x+an(k+1).*cos(2*k*t); y(k,:) = x; subplot(m+2,1,k+1) plot(t/pi,y(m+1,:)); hold on;
plot(t/pi,y(k,:)); hold off; grid on;
axis([-2.5 2.5 -0.2 1.2]); title(strcat('第',kk,'次谐波叠加')); xlabel('t/pi','Fontsize', 8);
kk = strcat(kk,'、',num2str(k+1)); end pause;
subplot(m+2,1,m+2)
plot(t/pi,y(1:m+1,:)); grid on;
axis([-2.5 2.5 -0.2 1.2]); title('各次次谐波叠加波形');
xlabel('t/pi','Fontsize', 8); % End
9.2 试用MATLAB编程会出9.1中所示周期信号的幅度频谱，要求交互式输入信号周期，观察分
析信号周期与频谱的关系。当周期T趋于无穷大时，频谱谱线将发生什么样的变化？
3
信号与线性系统分析实验
解：
syms T n t;
display('Please input the value of NF,T,tao'); Nf=input('Nf='); T=input('T=');
tao=input('tao=');
an=2/T*int(2/T*t*cos(2*pi*n*t/T),t,0,tao); bn=2/T*int(2/T*t*sin(2*pi*n*t/T),t,0,tao); a1n=zeros(1,Nf+1); b1n=zeros(1,Nf+1); cn=zeros(1,Nf+1); cn(1)=1/2; for i=2:Nf+1
a1n(i)=subs(an,n,i); b1n(i)=subs(bn,n,i);
cn(i)=(a1n(i)^2+b1n(i)^2)^(1/2); end
k=0:Nf;
stem(k,cn); hold on
plot(k,cn);
title(strcat('幅度频谱，周期、时域信号分别为',num2str(T),'和',num2str(tao))); xlabel(strcat('谐波次数',k));
运行：
Please input the value of m (傅立叶级数展开的项数) Please input the value of m (傅立叶级数展开的项数) Please input the value of NF,T,tao Nf=30 T=2*pi
tao=pi
tao=pi/2
4
信号与线性系统分析实验
tao=pi/4
固定tao=1, Nf=30
改变T分别为2*pi ,4*pi ,8*pi 图形分别如下： T=2*pi
T=8*pi
5
信号与线性系统分析实验
T=4*pi
9.3 已知周期锯齿脉冲信号如图所示（图略），用MATLAB绘制其频谱，要求交互式设置信
号周期和时域宽度，观察信号周期T及时域宽度对信号频谱的影响。
?1?x......(0?x??)解：?在一个周期内f(t)???
?(??x?T)?0..........?a0?而：f(t)???ancos(n?t)??bncos(n?t)
2n?1n?122首先求解其傅立叶系数：a0??f(t)dt;an??f(t)cos(n?t)dt
T0T02bn??f(t)sin(n?t)dt
T0
TTT编写.m文件如下：
disp('please input the value of T ,tao and Nf'); T=input('T=');
tao=input('tao='); Nf=input('Nf='); syms x n k Nn=32;
An=zeros(Nf+1,1); Bn=zeros(Nf+1,1); f=x/tao;
6
信号与线性系统分析实验
a0=2*int(f,x,0,tao)/T;
an=2*int(f*cos(n*x),x,0,tao)/T; bn=2*int(f*sin(n*x),x,0,tao)/T; An(1)=double(vpa(a0,Nn)); for k=1:Nf
An(k+1)=double(vpa(subs(an,n,k),Nn)); Bn(k+1)=double(vpa(subs(bn,n,k),Nn)); end
cn=sqrt(An.*An+Bn.*Bn); m=0:Nf;
stem(m,cn); hold on; plot(m,cn);
xlabel('·ù?è\\omega','fontsize',8);
7
信号与线性系统分析实验
其次固定tao，改变T:
8