傅里叶变换的应用matlab程序C语言程序

1、1 利用FFT计算连续时间信号的傅里叶变换设是连续时间信号，并假设时，则其傅里叶变换由下式给出令是一个固定的正实数，是一个固定的正整数。当时，利用FFT算法可计算。已知一个固定的时间间隔，选择足够小，使得每一个秒的间隔内，的变化很小，则式中积分可近似为 （27）假设足够大，对于所有的整数，幅值很小，则式（27）变为 （28）当时，式（28）两边的值为 （29）其中代表抽样信号的点。最后令，则上式变为 （30）首先用FFT算法求出，然后可用上式求出时的。应该强调的是，式（28）只是一个近似表示，计算得到的只是一个近似值。通过取更小的抽样间隔，或者增加点数，可以得到更精确的值。如果时，幅度谱很小，

2、对应于奈奎斯特抽样频率，抽样间隔选择比较合适。如果已知信号只在时间区间内存在，可以通过对时的抽样信号补零，使足够大。例1 利用FFT计算傅里叶变换如图12所示的信号其傅里叶变换为：利用下面的命令，可得到的近似值和准确值。 图12 连续时间信号x(t) N=input('Input N:');T=input('Input T:');%计算X(w)近似值t=0:T:2;x=t-1 zeros(1,N-length(t);X=fft(x);gamma=2*pi/(N*T);k=0:10/gamma;Xapp=(1-exp(-i*k*gamma*T)/(i*k*gamm

3、a)*X;%计算真实值X(w)w=0.05:0.05:10;Xact=exp(-i*w)*2*i.*(w.*cos(w)-sin(w)./(w.*w);plot(k*gamma,abs(Xapp(1:length(k),'o',w,abs(Xact);legend('近似值','真实值');xlabel('频率(rad/s)');ylabel('|X|')运行程序后输入N=128，T=0.1，此时，得到实际的和近似的傅里叶变换的幅度谱如图13所示，此时近似值已经相当准确。通过增加NT可以增加更多的细节，减少T使得到

4、的值更精确。再次运行程序后输入N=512，T=0.05，此时，得到实际的和近似的傅里叶变换的幅度谱如图14所示。图13 N=128，T=0.1时的幅度谱图14 N=512，T=0.05时的幅度谱2 利用FFT计算离散信号的线性卷积已知两个离散时间信号与，取，对和右端补零，使得 （31）利用FFT算法可以求得和的L点DFT，分别是和，利用DTFT卷积性质，卷积等于乘积的L点DFT反变换，这也可以通过FFT 算法得到。例2 利用FFT计算线性卷积已知，其中为单位阶跃序列，信号如图15所示。由于当时，很小，故可以取为17；N取10，。利用下面的Matlab命令，可得到、的卷积图形如图15所

5、示。subplot(3,1,1);n=0:16;x=0.8.n;stem(n,x);xlabel('n');ylabel('xn');subplot(3,1,2);n=0:15;y=ones(1,10) zeros(1,6);stem(n,y);xlabel('n');ylabel('yn')subplot(3,1,3);L=26;n=0:L-1;X=fft(x,L);Y=fft(y,L);Z=X.*Y;z=ifft(Z,L);stem(n,z);xlabel('n');ylabel('zn')图1

6、5 信号xn、yn及其卷积zn=xn*yn利用下面的Matlab命令，可得到信号xn、yn的幅度谱与相位谱如图16所示。subplot(2,2,1);L=26;k=0:L-1;n=0:16;x=0.8.n;X=fft(x,L);stem(k,abs(X);axis(0 25 0 5);xlabel('k');ylabel('|Xk|')subplot(2,2,2);stem(k,angle(X);axis(0 25 -1 1);xlabel('k');ylabel('Angle(Xk)(弧度)')subplot(2,2,3);y=

7、ones(1,10);Y=fft(y,L);stem(k,abs(Y);axis(0 25 0 10);xlabel('k');ylabel('|Yk|')subplot(2,2,4);stem(k,angle(Y);axis(0 25 -3 3);xlabel('k');ylabel('Angle(Yk)(弧度)')图16 信号xn、yn的幅度谱与相位谱3 利用FFT进行离散信号压缩利用FFT算法对离散信号进行压缩的步骤如下：1）通过采样将信号离散化；2）对离散化信号进行傅里叶变换；3）对变换后的系数进行处理，将绝对值小于某一阈

8、值的系数置为0，保留剩余的系数；4）利用IFFT算法对处理后的信号进行逆傅里叶变换。例3 对单位区间上的下列连续信号以采样频率进行采样，将其离散化为个采样值用FFT分解信号，对信号进行小波压缩，然后重构信号。令绝对值最小的80%系数为0，得到重构信号图形如图17 a)所示，均方差为0.0429，相对误差为0.0449；令绝对值最小的90%系数为0，得到重构信号图形如图17 b)所示，均方差为0.0610，相对误差为0.0638。 a) 绝对值最小的80%系数为0的重构信号（FFT） b) 绝对值最小的90%系数为0的重构信号（FFT）图17 用FFT压缩后的重构信号相关Matlab程序如下fu

9、nction wc=compress(w,r)%压缩函数compress.m%输入信号数据w,压缩率r%输出压缩后的信号数据if(r<0)|(r>1) error('r 应该介于0和1之间!');end;N=length(w);Nr=floor(N*r);ww=sort(abs(w);tol=abs(ww(Nr+1);wc=(abs(w)>=tol).*w;function unbiased_variance,error=fftcomp(t,y,r)%利用FFT做离散信号压缩%输入时间t,原信号y,以及压缩率r%输出原信号和压缩后重构信号的图像,以及重构均方差

10、和相对l2误差if(r<0)|(r>1) error('r 应该介于0和1之间!');end;fy=fft(y);fyc=compress(fy,r); %调用压缩函数compress.myc=ifft(fyc);plot(t,y,'r',t,yc,'b');legend('原信号','重构信号');unbiased_variance=norm(y-yc)/sqrt(length(t);error=norm(y-yc)/norm(y);输入以下Matlab命令：t=(0:255)/256;f=t+cos

11、(4*pi*t)+1/2*sin(8*pi*t);unbiased_variance,error=fftcomp(t,f,0.8)unbiased_variance = 0.0429error =0.0449如果用Harr尺度函数和Harr小波分解信号，对信号进行小波压缩，然后重构信号。令绝对值最小的80%系数为0，得到重构信号图形如图18 a)所示，均方差为0.0584，相对误差为0.0611；令绝对值最小的90%系数为0，得到重构信号图形如图18 b)所示，均方差为0.1136，相对误差为0.1190。 a) 绝对值最小的80%系数为0的重构信号（Harr） b) 绝对值最小的90%系数为

12、0的重构信号（Harr）图18 用Harr小波压缩后的重构信号相关Matlab程序如下function unbiased_variance,error=daubcomp(t,y,n,r)%利用Daubechies系列小波做离散信号压缩%输入时间t,原信号y,分解层数n,以及压缩率r%输出原信号和压缩后重构信号的图像,以及重构均方差和相对l2误差if(r<0)|(r>1) error('r应该介于0和1之间!');end;c,l=wavedec(y,n,'db1');cc=compress(c,r); %调用压缩函数compress.myc=waver

13、ec(cc,l,'db1');plot(t,y,'r',t,yc,'b');legend('原信号','重构信号');unbiased_variance=norm(y-yc)/sqrt(length(t);error=norm(y-yc)/norm(y);输入以下Matlab命令：t=(0:255)/256;f=t+cos(4*pi*t)+1/2*sin(8*pi*t);unbiased_variance,error=daubcomp(t,f,8,0.8)unbiased_variance = 0.0584erro

14、r = 0.0611结论：在信号没有突变、快变化或者大致上具有周期性的信号，用FFT可以处理得很好（甚至比小波还要好）。4 利用FFT对离散信号进行滤波利用FFT算法对信号进行滤波的步骤如下：1）通过采样将信号离散化；2）对离散化信号进行傅里叶变换；3）对变换后的系数进行处理，将绝对值大于某一阈值的系数置为0，保留剩余的系数；4）利用IFFT算法对处理后的信号进行逆傅里叶变换。例4 股票价格分析首先进入网址，点击网页底部位置的Download To Spreadsheet按钮，即可把以Excel表格格式存储的价格数据下载到本地计算机。表格从1列至第6列分别给出了从1996年4月12日至2007

15、年5月30日的交易期里每天的开盘价、最高价、最低价、收盘价、成交量以及趋势。数据下载完成后，需要颠倒顺序，使得最早时间的数据首先显示。然后另存到Matlab所在的目录中，并重新命名为“yhoodata.csv”。 本次分析选择开盘价，时间是从2007年1月1日至2007年5月30日的=102个交易日期。令代表一支股票的开盘价。为了便于分析，可以先从中减去跃变，得到，即 （32）输入以下命令，得到的频谱如图19所示。o=csvread('yhoodata.csv',2700,1,2700 1 2801 1)'N=102;for n=1:N x(n)=o(n)-o(1)-(

16、o(N)-o(1)/(N-1)*(n-1);endX=fft(x);k=0:N-1;stem(k,abs(X);axis(0 101 0 300);xlabel('k');ylabel('|Xk|')图19 xn的幅度谱可以看出上图中有5个较强谱分量，频率（）分别对应和。保留这5个频率分量的系数，将其他频率分量的系数置为0，然后再进行逆傅里叶变换，得到滤波后的近似值。输入如下Matlab程序，得到真实值与滤波后的近似值，如图20所示。plot(x);hold on;fliterX=X(1:2) 0 0 X(5) zeros(1,102-9) X(99) 0 0

17、X(102);fliterx=ifft(fliterX);plot(real(fliterx),'r');axis(1 102 0 7);xlabel('n');ylabel('xn的值和它的近似值');legend('xn真实值','xn近似值')图20 xn的真实值与滤波后的近似值从上图可以看出，滤波后的近似值既大致上保留了真实值的变化趋势，而且与其十分接近。与滤波前比较，滤波后的图形要比滤波前平滑得多。再由式（32）即可求得 （33）输入如下Matlab程序，画出真实开盘价与近似开盘价的图形。如图21所示，可

18、以看出是近似基础上的平滑值。plot(o);hold on;for n=1:N oapp(n)=fliterx(n)+o(1)+(o(N)-o(1)/(N-1)*(n-1);endplot(oapp,'r');axis(1 102 25 34);xlabel('n');ylabel('on的值和它的近似值');legend('on真实值','on近似值')图21 on的真实值与滤波后的近似值5 利用FFT提取离散信号中的最强正弦分量这里最强是指在信号中某个正弦分量的振幅远大于其它正弦分量的振幅。可以对求点来确定信号

19、中是否有最强的正弦分量。如果信号是连续时间形式的，首先还要对其进行抽样，得到离散时间形式的信号，根据Nyquist定理，抽样间隔应满，其中是中的最大频率分量。要判断信号中是否包含最强正弦分量，采样数据至少要包含该分量一个完整周期的数据，例5 太阳耀斑数据的分析太阳耀斑活动的周期是11年，这个事实可以通过提取耀斑数据的最强正弦分量加以证实。耀斑数据可以从比利时皇家天文台(Royal Observatory of Belgium)太阳耀斑数据索引中心（Sunspot Index Data Center, SIDC）网站下载。网址是：下载后的数据存放在文件“monthssn.dat”中，里面有四列数

20、据，第一年是日期，第三列是太阳耀斑的平均数，第四列平滑后太阳耀斑的平均数，可以得到从1749年到当前年月（2006年4月）的耀斑数据。本次分析选择第三列1981年1月作为开始日期，2005年12月作为结束日期，共25年300个月份的数据。为此先把相关数据复制到Excel表格的第一列中，然后保存到Matlab所在目录下，并命名为“sunspotdata.csv”。然后输入以下命令，得到耀斑曲线如图22所示。spd=csvread('sunspotdata.csv',0,0,0 0 299 0);plot(spd);grid;xlabel('月数');ylabel(

21、'耀斑平均数');axis(0 300 0 200)图22 1981年1月至2005年12月太阳耀斑的平均数由上图可见，太阳耀斑的活动确实具有周期性，但周期的准确值不明显。可以通过数第一个峰值和第二个峰值之间的月份来估计周期的值。查验表中的数据得，第一个峰值为200.3，出现在第116个月（1990年8月），第二个峰值为170.1，出现在第235个月（2000年7月），所以周期是235-116=119月，和实际值132月比较接近。下面利用来分析。首先，从图22中可以看出从到整个区间的平均值不为0，为了更方便地分析，需要减去该均值，得到结果如下： （33）然后对进行傅里叶变换，得到。在Matlab中输入以下命令，得到的幅度谱的图形如图23所示。x=spd-sum(spd)/300;X=fft(x);k=0:299;plot(k,abs(X