实验二 傅里叶分析及应用
一、实验目的
（一）掌握使用Matlab进行周期信号傅里叶级数展开和频谱分析
1、学会使用Matlab分析傅里叶级数展开，深入理解傅里叶级数的物理含义 2、学会使用Matlab分析周期信号的频谱特性
（二）掌握使用Matlab求解信号的傅里叶变换并分析傅里叶变换的性质
1、学会运用Matlab求连续时间信号的傅里叶变换 2、学会运用Matlab求连续时间信号的频谱图
3、学会运用Matlab分析连续时间信号的傅里叶变换的性质 （三） 掌握使用Matlab完成信号抽样并验证抽样定理
1、学会运用MATLAB完成信号抽样以及对抽样信号的频谱进行分析 2、学会运用MATLAB改变抽样时间间隔，观察抽样后信号的频谱变化 3、学会运用MATLAB对抽样后的信号进行重建
二、实验条件
安装winXP系统的电脑一台、matlab 7.0软件
三、实验内容
1、已知周期三角信号如下图所示[注：图中时间单位为：毫秒(ms)]：
（1）试求出该信号的傅里叶级数[自己求或参见课本P112或P394]，利用Matlab编程实现其各次谐波[如1、3、5、13、49]的叠加，并验证其收敛性；
1
解： 命令文件： clear all;close all;clc; t=-10:0.01:10; omega=pi;
y=abs(sawtooth(pi*0.5*t,0.5)); plot(t,y),grid on; axis([-10,10,0,3]); n_max=[1,3,5,13,49]; N=length(n_max); for k=1:N
n=1:2:n_max(k); b=4./((pi*n).^2); x=b*cos(omega*n'*t); figure; plot(t,y); hold on; x=x+1/2; plot(t,x); hold off; axis([-10,10,0,3]);
title(['最大谐波数=',num2str(n_max(k))]); end 图像：
2
32.521.510.50-10-8-6-4-20246810
最大谐波数=132.521.510.50-10-8-6-4-20246810
3
最大谐波数=332.521.510.50-10-8-6-4-20246810
最大谐波数=532.521.510.50-10-8-6-4-20246810
4
最大谐波数=1332.521.510.50-10-8-6-4-20246810
最大谐波数=4932.521.510.50-10-8-6-4-20246810
（2）用Matlab分析该周期三角信号的频谱[三角形式或指数形式均可]。当周期三角信号的周期（如由2ms1ms或由2ms4ms）和宽度（如2ms1ms）分别变化时，试观察
5
分析其频谱的变化。
解： 周期为2ms： 命令文件： clear all;close all;clc; dt=0.01; t=-4:dt:4;
ft=(t>=-1&t<0).*(t+1)+(t>0&t<=1).*(1-t); %subplot(2,1,1) %plot(t,ft);grid on n=2000; k=-n:n; w=pi*k/(n*dt); f=dt*ft*exp(-i*t'*w); f=abs(f); %subplot(2,1,2) plot(w,f);
axis([-20 20 0 1.1]),grid on; 图像：
6
10.80.60.40.20-20-15-10-505101520
当周期为1ms时： 命令文件： clear all;close all;clc; dt=0.01; t=-4:dt:4;
ft=(t>=-0.5&t<0).*(t+1)+(t>0&t<=0.5).*(1-t); %subplot(2,1,1) %plot(t,ft);grid on n=2000; k=-n:n; w=pi*k/(n*dt); f=dt*ft*exp(-i*t'*w); f=abs(f); %subplot(2,1,2) plot(w,f);
axis([-20 20 0 1.1]),grid on;
7
图像：
10.80.60.40.20-20-15-10-505101520
宽度与周期均变为1ms，所以只有两个图
2、分别利用Matlab符号运算求解法和数值计算法求下图所示信号的FT，并画出其频谱图（包括幅度谱和相位谱）[注：图中时间单位为：毫秒(ms)]。
解：
符号运算求解法： 命令文件： %符号运算求解法 clear all;close all;clc;
8
ft=sym('(t+2)*(Heaviside(t+2)-Heaviside(t+1))+Heaviside(t+1)-Heaviside(t-1)+(2-t)*(Heaviside(t-1)-Heaviside(t-2))');
Fw=fourier(ft); subplot(2,1,1) axis([-8,8,0,3]); ezplot(abs(Fw)),grid on title('幅度谱')
phase=atan(imag(Fw)/real(Fw)); subplot(2,1,2) ezplot(phase),grid on title('相位谱') 图像：
幅度谱3210-6-4-20w相位谱24610.50-0.5-1-6-4-20w246
9
数值计算法： 命令文件： clear all;close all;clc; dt=0.01; t=-4:dt:4;
ft=(t>=-2&t<=-1).*(t+2)+(t>=1&t<=2).*(2-t)+(t>-1&t<1); n=2000; k=-n:n; w=pi*k/(n*dt); f=dt*ft*exp(-i*t'*w); f=abs(f); plot(w,f);
axis([-10 10 0 3.2]),grid on; 图像：
10
32.521.510.50-10-8-6-4-20246810
3、试用Matlab命令求F(j)104的傅里叶反变换，并绘出其时域信号图。 -3j5j [注意：(1)写代码时ji] 解：命令文件： clear all;close all;clc; t=sym('t');
Fw=sym('10/(3+i*w)-4/(5+i*w)'); ft=ifourier(Fw); ezplot(ft),grid on axis([-1 3 -1 7]); 图像：
11
2 heaviside(x) (-2 exp(-5 x)+5 exp(-3 x))76543210-1-1-0.500.51x1.522.53
4、已知门函数自身卷积为三角波信号，试用Matlab命令验证FT的时域卷积定理。 [注：即验证门函数FT的平方与相应三角波信号的FT后结果是否一致，可结合频谱图观察分析]
解： 函数文件： function f=uCT(t) f=(t>=0) 命令文件：
%将门函数先进行时域卷积运算，再将卷积后的结果做傅里叶变换 clear all;close all;clc; dt=0.01; t=-2:dt:2.5;
12
f1=uCT(t+0.5)-uCT(t-0.5); f=conv(f1,f1)*dt; ft=sym('f'); Fw=fourier(ft) 结果为： Fw =
2*i*pi*dirac(1,w)
%将一个门函数先进行傅里叶变换，再将结果与自身相乘 clear all;close all;clc; dt=0.01; t=-2:dt:2.5;
f1=uCT(t+0.5)-uCT(t-0.5); ft=sym('f1'); Fw=fourier(ft) Fw=Fw*Fw 结果为： Fw =
2*i*pi*dirac(1,w)
13
Fw =
-4*pi^2*dirac(1,w)^2
Hz，5、设有两个不同频率的余弦信号，频率分别为f1100Hz，f23800；现在使
Hz对这三个信号进行抽样，使用MATLAB命令画出各抽样信号的波用抽样频率fs4000形和频谱，并分析其频率混叠现象[建议：抽样信号的频谱图横坐标范围不小于-10000Hz~10000Hz或-20000*pi~20000*pi rad/s]。
解：
100HZ命令文件： clear all;close all;clc; time=2.5*10^(-4); dt=0.000001; t1=0:dt:0.01;
ft=sin(2*pi*100*t1).*(t1>=0); subplot(2,2,1); plot(t1,ft),grid on axis([0,0.01,-1.1,1.1]) xlabel('time(sec)'),ylabel('f(t)') title('100HZ正弦信号') n=500; k=-n:n;
14
w=pi*k/(n*dt); fw=dt*ft*exp(-i*t1'*w); subplot(2,2,2); plot(w,abs(fw)),grid on %axis([-4 4 0 1.1*pi]); xlabel('omega'),ylabel('f(w)') title('100HZ正弦信号的频谱') t2=0:time:0.01;
fst=sin(2*pi*100*t2).*(t2>=0); subplot(2,2,3); plot(t2,fst,':'),hold on stem(t2,fst),grid on
title('100HZ抽样后的信号'), hold off fsw=time* fst*exp(-i*t2'*w); subplot(2,2,4); plot(w,abs(fsw)),grid on title('100HZ抽样后的频谱') 100HZ
图
15
像
：
100HZ正弦信号10.5)(tf0-0.5-100.0050.01time(sec)100HZ抽样后的信号10.50-0.5-100.0050.013800HZ命令文件： clear all;close all;clc; time=1/4000; dt=0.00001; t1=-0.0003:dt:0.0003; ft=sin(2*pi*3800*t1); subplot(2,2,1); plot(t1,ft),grid on %axis([-4 4 -1.1 1.1]) xlabel('time(sec)'),ylabel('f(t)') title('3800HZ正弦信号') n=500;
x 10100HZ-16正弦信号的频谱32)w(f10-4-2024omegax 106x 10-14100HZ抽样后的频谱10.50-4-2024x 10616
k=-n:n; w=pi*k/(n*dt); fw=dt*ft*exp(-i*t1'*w); subplot(2,2,2); plot(w,abs(fw)),grid on
axis([-100000 100000 0 3*10^(-4)]); xlabel('omega'),ylabel('f(w)') title('3800HZ正弦信号的频谱') t2=-0.0003:time:0.0003; fst=sin(2*pi*3800*t2); subplot(2,2,3); plot(t2,fst,':'),hold on stem(t2,fst),grid on
axis([-0.0004 0.0004 -1.2 1.2]) title('3800HZ抽样后的信号'), hold off fsw=time* fst*exp(-i*t2'*w); subplot(2,2,4); plot(w,abs(fsw)),grid on
axis([-200000 200000 0 1*10^(-3)]) title('3800HZ抽样后的频谱') 3800HZ图像：
17
3800HZ正弦信号10.50-0.5-1-4-20time(sec)24-432-43800HZ正弦信号的频谱x 10f(w)10-1f(t)-0.5x 103800HZ抽样后的信号10.50-0.5-1-4-202x 101x 10-3x 103800HZ抽样后的频谱0omega0.5150.54-40-2-101x 1025
6、结合抽样定理，利用MATLAB编程实现Sa(t)信号经过冲激脉冲抽样后得到的抽样信号
fst及其频谱[建议：冲激脉冲的周期分别取4*pi/3 s、pi s、2*pi/3 s三种情况对比]，并利
用fst构建Sa(t)信号。）
解：
冲激脉冲的周期取4*pi/3 s：
命令文件：
clear all;close all;clc; Ts=2/3;
t1=-5:0.01:5; ft=sinc(t1); subplot(2,2,1) plot(t1,ft),grid on axis([-6 6 -0.5 1.2]) title('Sa(t)') N=500; k=-N:N;
W=pi*k/(N*0.01);
Fw=0.01*ft*exp(-j*t1'*W); subplot(2,2,2)
plot(W,abs(Fw)),grid on axis([-30 30 -0.05 1.5])
18
title('Sa(t)频谱') t2=-5:Ts:5; fst=sinc(t2); subplot(2,2,3)
plot(t1,ft,':'),hold on stem(t2,fst),grid on axis([-6 6 -0.5 1.2]) title('抽样后信号')
Fsw=Ts*fst*exp(-j*t2'*W); subplot(2,2,4)
plot(W,abs(Fsw)),grid on axis([-50 50 -0.05 1.5]) title('抽样后信号频谱') 图像：
Sa(t)10.50-0.5-505抽样后信号10.50-0.5-505冲击脉冲的周期取pi： 命令文件：
clear all;close all;clc; Ts=1/2;
t1=-5:0.01:5; ft=sinc(t1); subplot(2,2,1) plot(t1,ft),grid on axis([-6 6 -0.5 1.2]) title('Sa(t)')
Sa(t)频谱1.510.50-20020抽样后信号频谱1.510.50-5005019
N=500; k=-N:N;
W=pi*k/(N*0.01);
Fw=0.01*ft*exp(-j*t1'*W); subplot(2,2,2)
plot(W,abs(Fw)),grid on axis([-30 30 -0.05 1.5]) title('Sa(t)频谱') t2=-5:Ts:5; fst=sinc(t2); subplot(2,2,3)
plot(t1,ft,':'),hold on stem(t2,fst),grid on axis([-6 6 -0.5 1.2]) title('抽样后信号')
Fsw=Ts*fst*exp(-j*t2'*W); subplot(2,2,4)
plot(W,abs(Fsw)),grid on axis([-50 50 -0.05 1.5]) title('抽样后信号频谱') 图像：
Sa(t)10.50-0.5-505抽样后信号10.50-0.5-505冲击脉冲的周期取2/3*pi: 命令文件：
clear all;close all;clc;
Sa(t)频谱1.510.50-20020抽样后信号频谱1.510.50-5005020
Ts=1/3;
t1=-5:0.01:5; ft=sinc(t1); subplot(2,2,1) plot(t1,ft),grid on axis([-6 6 -0.5 1.2]) title('Sa(t)') N=500; k=-N:N;
W=pi*k/(N*0.01);
Fw=0.01*ft*exp(-j*t1'*W); subplot(2,2,2)
plot(W,abs(Fw)),grid on axis([-30 30 -0.05 1.5]) title('Sa(t)频谱') t2=-5:Ts:5; fst=sinc(t2); subplot(2,2,3)
plot(t1,ft,':'),hold on stem(t2,fst),grid on axis([-6 6 -0.5 1.2]) title('抽样后信号')
Fsw=Ts*fst*exp(-j*t2'*W); subplot(2,2,4)
plot(W,abs(Fsw)),grid on axis([-50 50 -0.05 1.5]) title('抽样后信号频谱') 图像:
Sa(t)1.510.50-0.510.50-50抽样后信号1.510.50-0.510.50-505-20020Sa(t)频谱抽样后信号频谱-505050
四、实验结论和讨论
本次试验帮助我更好的理解了傅里叶变换的概念，对抽样定理也有了很好的领悟，主要不足有一些概念理解不清，导致试验很难进行。
21
五、实验思考
以后将加强概念的理解，熟悉matlab的操作！
22
