傅里叶变换到计算机实现.doc

傅里叶变换到计算机实现 2013/8/16 Guan Jun 就拿我自身的例子来说，开始接触FFT（快速傅里叶变换）的时候并不是很熟悉，但是这种计算方法的确实很好用。那么，这个doc我想说的就是，如何从三角变换到FFT。 ,这是说一个周期性函数(T1)可以分解为不同频率的三角函数的叠加，带到原函数中，经过整理，令，再把的表达式（高数书或者信号与系统说的很清楚）带入中，我们就可以得到。以上是周期性函数的傅里叶变换，注意的是画出来的图是：在轴上频率的坐标为,即一系列间隔为的点，另外也就是说，周期函数的傅里叶变换为频域之后，是分立的频谱，不是连续的。举个栗子，函数是周期性函数吧，其频率（角频率）为，也可写成1，也就是在会有值，其余地方就没有。其实到这里，真的不难，因为求也就是带入公式的事么，不借助软件我们都能算好。但是，偏偏有那么一些人没事干非要去研究非周期性函数的傅里叶函数，然后搞出一大堆理论，让我们去学 废话不多说，如果是非周期性，是不是可以理解为周期无限大？这里的非周期函数也可由周期函数组成，例如在上，其余等0.这是不是非周期性函数？答案很显然.如果非周期性，那么公式不再适用，为什么？这得问数学系的人了。怎么办，把公式变变，移到左边，(此时频谱是连续的了，为什么，我也不晓得)那么我们就将看到最为熟悉的函数：，(也有书本写成：，).就是把，而自此，我们就开始学习一大堆公式，性质啊,我觉得这些性质不是不重要，而是没有实际的应用！为什么我这么说，因为我们用傅里叶变换，是为了什么？服务于我们的数据，没错，是数据！一堆数据给你，你能看出这函数包含的频率？你能提炼出原函数吗？Okay，你什么都没有，怎么办，望洋兴叹。 最近写的论文中，我就用到了FFT，我有图像的曲线，有曲线的数据，而且曲线明显是正余弦函数(只相差相位).大概的频率我也能看出来，但是！这个曲线并不完美，有瑕疵，但是我束手无策，这时计算机粉墨登场了，经过分析我也看出原来还是有很小的其他频率成分包含在里面。也许对傅里叶变换感兴趣的童鞋看过不少人的介绍，说时间连续，时间不连续，频谱连续，频谱不连续。22=4，这4种绕来绕去足以崩溃你（这里崩溃作动词).其实，时间连续，就是我上面讲的两种，但一个是周期性函数，一个是非周期，对应的频谱就是分立，连续。那么时间（有时候不一定是时间，也可能是位置）不连续怎么办，其实大多数应用的就是这种方法，就是我们说的采谱，说简单点就是每隔一段时间（距离）采一个点，采点间隔相同，一个点一个值.有时候觉得，数学应该是自然学科中最难的一门了.感兴趣的童鞋应该去找找资料，连续时间的傅里叶变换到离散时间傅里叶变换的过程（全是一个个公式的推导，很值得细究,我是没怎么看懂).这里，频谱连续我不想说了,说了我觉得也没用。因为,计算机绘图是画一系列点，然后把这些点连起来，点越来图像越正确，那么连续谱也可以用离散的谱画出来，只是画点的时候取密一点，画的才真,至于取点的问题最后说.先说明几个量,指的是采点的周期和频率，共取点N个,令第一个取点的时间为0，第二个就是，第三个就是，最后一个就是，那你说我取N+1个点行不行，也就是到，都可以,如果你把第一个点变成都行.接下来，说明一点，画出的频谱图是分立的点，点与点之间的间距是，而且(为什么？自己查证).在这里我需要说明的就是是你自己选取的，尽量大一点吧.时间是0,2,3,(N-1),即，对应的值我们就可以写成,当然了，我们只是写成这种形式而已,算的时候值往那个位置填即可.画频谱图的时候也是从,即，我们有公式,细心的同学肯定会发现书上跟我的公式不同，但是我这是最根本的公式，书上的公式是把乘上后的结果.关于取点呢，我想说的是,在不变的情况下，如果减小，就要变大N，就是增大取点的时间或者增大取点的位置长度.而值取的越小，就能更好地反应原函数的真实频率，举个栗子，原函数有一个频率值位于0.21，如果等于0.02，也许就不准确，而0.01()就会非常准确，当然这些都得看实际情况了.本文写的不足之处，多多批评，最后我把我老师的一个实例函数放到这里供大家更好理解,这个函数最后我加了用还原原函数的程序（Matlab软件编的，而我多次提到的FFT就是在编程的过程中用一系列方法去简化步骤，让计算机更快的算出结果，有兴趣的同学可以去看看源程序，很复杂，我是没打算掌握，因为我懒）%这个函数主要是进行FFT变换和傅里叶变%换之间(上文的蓝色公式)的验证，原函数%为f(t)=2*e-3t，傅里叶变换后的函数为%F(w)=2/(3+w*j)clc,clearN=1000;t=linspace(0,10,N);Ts=t(2)-t(1);fs=Ts-1;for k=0:N-1 y=0; for n=0:N-1 y=y+2*exp(-3*n*Ts)*exp(-j*2*pi*n*k/N); end F1(k+1)=y/fs;endf1=fs/N;w1=2*pi*f1;w=(0:N/2)*w1;F2=2./(w*j+3);abs_F2=abs(F2);abs_F1=abs(F1);plot(w,abs_F1(1:N/2+1),b)hold on, plot(w,abs_F2(1:N/2+1),r-)figure, plot(t,2*exp(-