实验二_用matlab实现傅立叶变换

1、实 验 报 告2010年 4月 15日实验二 用matlab实现傅立叶变换实验目的： 1.掌握傅立叶数值实现方法（矩阵算法）2.生成连续周期信号，掌握程序优化技巧3.对于自定义函数参数有效性的检查4.复习并巩固“信号与系统”相关知识内容，学习用matlab实现问题实验环境：运行于Matlab7.6环境实验内容： 本次实验参照信号与系统“Matlab综合实验”55页课后练习习题：1. 如图4.4所示锯齿波信号，分别去一个周期的抽样数据X1（t），0<=t<=1和五个周期的数据X（t），0<=t<5,计算其傅立叶变换X1（w）和X（w），比较有和不同并解释原因。图4.4 练

2、习题2图 编程如下：方法1：%计算单位锯齿波和五个周期波形的傅立叶变换 %解法1：基本用循环实现数值的计算；对于5个周期锯齿波用内外循环来生成实现T1=1; %一个周期锯齿波N1=10000;t1=linspace(0,T1-T1/N1,N1'f1=0*t1;f1=1-2*t1;OMG=32*pi; %频率抽样区间K1=100; %频率抽样点数omg=linspace(-OMG/2,OMG/2-OMG/K1,K1'X1=0*omg;for k=1:K1 %求解五个周期函数的傅里叶变换系数for n=1:N1X1(k=X1(k+T1/N1*f1(n*exp(-j*omg(k*t1

3、(n;endendfs1=0*t1;for n=1:N1 %通过傅里叶逆变换还原原函数for k=1:K1fs1(n=fs1(n+OMG/2/pi/K1*X1(k*exp(j*omg(k*t1(n;endendT2=5; %五个周期锯齿波N2=10000;t2=linspace(0,T2-T2/N2,N2'f2=0*t2;t3=linspace(0,T2/5-T2/N2,N2/5' %先定义一个周期内的锯齿波变量抽样值f3=0*t3; %初始化一个周期的函数抽样值f3=1-2*t3; %表示出一个周期内函数抽样值for s=0:4 %将一个周期锯齿波平移到五个周期，通过循环控制

4、for a=1:N2/5f2(2000*s+a=f3(a;endend %将函数拓展表示为五个周期X2=0*omg;for k=1:K1 %求解五个周期函数的傅里叶变换系数for n=1:N2X2(k=X2(k+T2/N2*f2(n*exp(-j*omg(k*t2(n;endendfs2=0*t2;for n=1:N1 %通过傅里叶逆变换还原原函数for k=1:K1fs2(n=fs2(n+OMG/2/pi/K1*X2(k*exp(j*omg(k*t2(n;endendfigure;subplot(2,2,1;plot(omg,abs(X1,'r' %以幅度频谱画图xlabel

5、('Frequency',ylabel('Amplitude'title('单个锯齿波的幅频曲线'subplot(2,2,2;plot(t1,fs1,'r'xlabel('Second(s',ylabel('Amplitude'title('由频域还原时域函数'subplot(2,2,3;plot(omg,abs(X2,'r'xlabel('Frequency',ylabel('Amplitude'title('五个周期锯齿波的

6、幅频曲线'subplot(2,2,4;plot(t2,fs2,'r'xlabel('Second(s',ylabel('Amplitude'title('由频域还原时域函数'相关曲线：方法2：%计算单位锯齿波和五个周期波形的傅立叶变换 %解法2：数值算法用矩阵实现，大大加快了运行速度；并且直接调用“sawtooth”生成5个周期的锯齿波T1=1; %单个周期时域范围 N1=10000; %时域抽样点数 t1=linspace(0,T1-T1/N1,N1' %生成抽样时间点f1=1-2*t1; %生成抽样函数值 OM

7、G=32*pi; %频域范围K1=100; %频域抽样点数omg=linspace(-OMG/2,OMG/2-OMG/K1,K1' %生成抽样频率点X1=T1/N1*exp(-j*kron(omg,t1.'*f1; %傅里叶正变换求解傅里叶系数fs1=OMG/2/pi/K1*exp(j*kron(t1,omg.'*X1; %傅里叶逆变换还原时域函数T2=5; %五个周期时域范围N2=10000; %时域抽样点数t2=linspace(0,T2-T2/N2,N2' %生成抽样时间点fs2=0*t2;f2=sawtooth(t2*2*pi,0; %生成五个周期的锯齿

8、波 X2=T2/N2*exp(-j*kron(omg,t2.'*f2; %傅里叶正变换求解傅里叶系数fs2=fs2+OMG/2/pi/K1*exp(j*kron(t2,omg.'*X2; %傅里叶逆变换还原时域函数figure; %生成一个2*2矩阵子图subplot(2,2,1;plot(omg,abs(X1,'r' %一个周期时的频谱图xlabel('Frequency',ylabel('Amplitude'title('单个锯齿周期幅频特性曲线'subplot(2,2,2; plot(t1,fs1,'

9、r' %还原的时域函数 xlabel('Time',ylabel('Amplitude'title('Function after recovered'subplot(2,2,3;plot(omg,abs(X2,'r' %五个周期时的频谱图 xlabel('Frequency',ylabel('Amplitude'title('五个锯齿周期幅频特性曲线'subplot(2,2,4;plot(t2,fs2,'r' %还原的时域函数xlabel('Time

10、',ylabel('Function after recovered'title('Function after recovered'相关曲线：2. 请编写函数F=fsana(t,f,N,计算周期信号f的前N个指数形式的傅立叶级数系数，t表示f对应的抽样时间（均为一个周期）；再编写函数f=fssyn（F，t），由傅立叶级数系数F合成抽样时间t对应的函数。设计信号验证这两个是否正确。3. 利用fsana和fssyn计算练习题1中X1（t）的前10个傅立叶级数系数F，请用这些系数合成周期为0.5的锯齿波，0<=t<2并绘图。解：这两题验证函数直接

11、用练习1所给函数.编程如下：编写fsana函数：%函数功能：定义一个功能函数来实现周期函数的傅里叶正变换 t表示f表示的抽样时间 f为一个周期内抽样函数值 N为前几个指数形式的傅立叶级数系数function F=fsana(t,f,Nomg1=2*pi/(max(t-min(t; k=-N:N'F=1/length(t*exp(-j*kron(k*omg1,t.'*f;编写fssyn函数：%函数功能：定义功能函数fssyn由傅立叶级数系数F合成时间t对应的函数function f=fssyn(F,tomg1=2*pi/(max(t-min(t;N=floor(length(F/

12、2;k=-N:N;f=exp(j*kron(t,k*omg1*F;验证所写函数：clcclearclose allT1=1; %一个周期时域范围N1=256; %时域抽样点数 t=linspace(0,T1-T1/N1,N1' %生成抽样时间点f=1-2*t; %生成抽样函数值 plot(fN=25F1=fsana(t,f,N %调用fsana函数求解前N项傅立叶级数系数figure(1stem(abs(F1,'s'title('前N项傅立叶级数系数幅度曲线'f2=fssyn(F1,t %调用fssyn函数求原时域函数f3=fssyn1(F1,t %用所

13、求傅立叶系数合成原基频2倍的函数 调用修改后的fssyn函数（修改基频）figure(2plot(t,f,'b',t,f2,'k'xlabel('times',ylabel('Amplitude'legend('origin','recover'title('傅立叶逆变换后时域函数与原时域函数对比'figure(3plot(t,f3xlabel('times',ylabel('Amplitude'title('时域压缩 基频变为2倍后的时域函数

14、'注：练习3要生成0.5周期的锯齿波，和原锯齿波相比是基频变为2倍，同样改写fssyn函数，修改基频，调用fssyn1函数：function f=fssyn1(F,tomg1=2*2*pi/(max(t-min(t;N=floor(length(F/2;k=-N:N;f=exp(j*kron(t,k*omg1*F;相关曲线：实验小结：这个多礼拜的学习，我进一步对于matlab的数值算法有了自己的了解。很多数学运算，可以类似地通过高级语言的方法实现，比如循环，可是这样花费的时间比较久，相应的代码也比较长。Matlab软件库函数中有很多基于矩阵算法，因此，将问题通过矩阵来实现，将大大缩短程序和运行时间，通常也比较读懂（当然是在练习后）。这次的学习，我另一方面明白了对于自定义函数的调用，它的文件名称应该与函数名称一样。更重要的是，我学习到了一些基本的Matlab调试方法，在确认所使用算法和数学公式无误的情况下，我通