傅里葉變換的應(yīng)用,matlab程序,C語言程序

1、1 利用FFT計(jì)算連續(xù)時(shí)間信號(hào)的傅里葉變換設(shè)是連續(xù)時(shí)間信號(hào)，并假設(shè)時(shí)，則其傅里葉變換由下式給出令是一個(gè)固定的正實(shí)數(shù)，是一個(gè)固定的正整數(shù)。當(dāng)時(shí)，利用FFT算法可計(jì)算。已知一個(gè)固定的時(shí)間間隔，選擇足夠小，使得每一個(gè)秒的間隔內(nèi)，的變化很小，則式中積分可近似為 （27）假設(shè)足夠大，對(duì)于所有的整數(shù)，幅值很小，則式（27）變?yōu)?（28）當(dāng)時(shí)，式（28）兩邊的值為 （29）其中代表抽樣信號(hào)的點(diǎn)。最后令，則上式變?yōu)?（30）首先用FFT算法求出，然后可用上式求出時(shí)的。應(yīng)該強(qiáng)調(diào)的是，式（28）只是一個(gè)近似表示，計(jì)算得到的只是一個(gè)近似值。通過取更小的抽樣間隔，或者增加點(diǎn)數(shù)，可以得到更精確的值。如果時(shí)，幅度譜很小，

2、對(duì)應(yīng)于奈奎斯特抽樣頻率，抽樣間隔選擇比較合適。如果已知信號(hào)只在時(shí)間區(qū)間內(nèi)存在，可以通過對(duì)時(shí)的抽樣信號(hào)補(bǔ)零，使足夠大。例1 利用FFT計(jì)算傅里葉變換如圖12所示的信號(hào)其傅里葉變換為：利用下面的命令，可得到的近似值和準(zhǔn)確值。 圖12 連續(xù)時(shí)間信號(hào)x(t) N=input('Input N:');T=input('Input T:');%計(jì)算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;%計(jì)算真實(shí)值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('近似值','真實(shí)值');xlabel('頻率(rad/s)');ylabel('|X|')運(yùn)行程序后輸入N=128，T=0.1，此時(shí)，得到實(shí)際的和近似的傅里葉變換的幅度譜如圖13所示，此時(shí)近似值已經(jīng)相當(dāng)準(zhǔn)確。通過增加NT可以增加更多的細(xì)節(jié)，減少T使得到

4、的值更精確。再次運(yùn)行程序后輸入N=512，T=0.05，此時(shí)，得到實(shí)際的和近似的傅里葉變換的幅度譜如圖14所示。圖13 N=128，T=0.1時(shí)的幅度譜圖14 N=512，T=0.05時(shí)的幅度譜2 利用FFT計(jì)算離散信號(hào)的線性卷積已知兩個(gè)離散時(shí)間信號(hào)與，取，對(duì)和右端補(bǔ)零，使得 （31）利用FFT算法可以求得和的L點(diǎn)DFT，分別是和，利用DTFT卷積性質(zhì)，卷積等于乘積的L點(diǎn)DFT反變換，這也可以通過FFT 算法得到。例2 利用FFT計(jì)算線性卷積已知，其中為單位階躍序列，信號(hào)如圖15所示。由于當(dāng)時(shí)，很小，故可以取為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 信號(hào)xn、yn及其卷積zn=xn*yn利用下面的Matlab命令，可得到信號(hào)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 信號(hào)xn、yn的幅度譜與相位譜3 利用FFT進(jìn)行離散信號(hào)壓縮利用FFT算法對(duì)離散信號(hào)進(jìn)行壓縮的步驟如下：1）通過采樣將信號(hào)離散化；2）對(duì)離散化信號(hào)進(jìn)行傅里葉變換；3）對(duì)變換后的系數(shù)進(jìn)行處理，將絕對(duì)值小于某一閾

8、值的系數(shù)置為0，保留剩余的系數(shù)；4）利用IFFT算法對(duì)處理后的信號(hào)進(jìn)行逆傅里葉變換。例3 對(duì)單位區(qū)間上的下列連續(xù)信號(hào)以采樣頻率進(jìn)行采樣，將其離散化為個(gè)采樣值用FFT分解信號(hào)，對(duì)信號(hào)進(jìn)行小波壓縮，然后重構(gòu)信號(hào)。令絕對(duì)值最小的80%系數(shù)為0，得到重構(gòu)信號(hào)圖形如圖17 a)所示，均方差為0.0429，相對(duì)誤差為0.0449；令絕對(duì)值最小的90%系數(shù)為0，得到重構(gòu)信號(hào)圖形如圖17 b)所示，均方差為0.0610，相對(duì)誤差為0.0638。 a) 絕對(duì)值最小的80%系數(shù)為0的重構(gòu)信號(hào)（FFT） b) 絕對(duì)值最小的90%系數(shù)為0的重構(gòu)信號(hào)（FFT）圖17 用FFT壓縮后的重構(gòu)信號(hào)相關(guān)Matlab程序如下fu

9、nction wc=compress(w,r)%壓縮函數(shù)compress.m%輸入信號(hào)數(shù)據(jù)w,壓縮率r%輸出壓縮后的信號(hào)數(shù)據(jù)if(r<0)|(r>1) error('r 應(yīng)該介于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做離散信號(hào)壓縮%輸入時(shí)間t,原信號(hào)y,以及壓縮率r%輸出原信號(hào)和壓縮后重構(gòu)信號(hào)的圖像,以及重構(gòu)均方差

10、和相對(duì)l2誤差if(r<0)|(r>1) error('r 應(yīng)該介于0和1之間!');end;fy=fft(y);fyc=compress(fy,r); %調(diào)用壓縮函數(shù)compress.myc=ifft(fyc);plot(t,y,'r',t,yc,'b');legend('原信號(hào)','重構(gòu)信號(hào)');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尺度函數(shù)和Harr小波分解信號(hào)，對(duì)信號(hào)進(jìn)行小波壓縮，然后重構(gòu)信號(hào)。令絕對(duì)值最小的80%系數(shù)為0，得到重構(gòu)信號(hào)圖形如圖18 a)所示，均方差為0.0584，相對(duì)誤差為0.0611；令絕對(duì)值最小的90%系數(shù)為0，得到重構(gòu)信號(hào)圖形如圖18 b)所示，均方差為0.1136，相對(duì)誤差為0.1190。 a) 絕對(duì)值最小的80%系數(shù)為0的重構(gòu)信號(hào)（Harr） b) 絕對(duì)值最小的90%系數(shù)為

12、0的重構(gòu)信號(hào)（Harr）圖18 用Harr小波壓縮后的重構(gòu)信號(hào)相關(guān)Matlab程序如下function unbiased_variance,error=daubcomp(t,y,n,r)%利用Daubechies系列小波做離散信號(hào)壓縮%輸入時(shí)間t,原信號(hào)y,分解層數(shù)n,以及壓縮率r%輸出原信號(hào)和壓縮后重構(gòu)信號(hào)的圖像,以及重構(gòu)均方差和相對(duì)l2誤差if(r<0)|(r>1) error('r應(yīng)該介于0和1之間!');end;c,l=wavedec(y,n,'db1');cc=compress(c,r); %調(diào)用壓縮函數(shù)compress.myc=waver

13、ec(cc,l,'db1');plot(t,y,'r',t,yc,'b');legend('原信號(hào)','重構(gòu)信號(hào)');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結(jié)論：在信號(hào)沒有突變、快變化或者大致上具有周期性的信號(hào)，用FFT可以處理得很好（甚至比小波還要好）。4 利用FFT對(duì)離散信號(hào)進(jìn)行濾波利用FFT算法對(duì)信號(hào)進(jìn)行濾波的步驟如下：1）通過采樣將信號(hào)離散化；2）對(duì)離散化信號(hào)進(jìn)行傅里葉變換；3）對(duì)變換后的系數(shù)進(jìn)行處理，將絕對(duì)值大于某一閾值的系數(shù)置為0，保留剩余的系數(shù)；4）利用IFFT算法對(duì)處理后的信號(hào)進(jìn)行逆傅里葉變換。例4 股票價(jià)格分析首先進(jìn)入網(wǎng)址 To Spreadsheet按鈕，即可把以Excel表格格式存儲(chǔ)的價(jià)格數(shù)據(jù)下載到本地計(jì)算機(jī)。表格從1列至第6列分別給出了從1996年4月12日至2007年5月30日的交易期里每天的開盤價(jià)、

15、最高價(jià)、最低價(jià)、收盤價(jià)、成交量以及趨勢(shì)。數(shù)據(jù)下載完成后，需要顛倒順序，使得最早時(shí)間的數(shù)據(jù)首先顯示。然后另存到Matlab所在的目錄中，并重新命名為“yhoodata.csv”。 本次分析選擇開盤價(jià)，時(shí)間是從2007年1月1日至2007年5月30日的=102個(gè)交易日期。令代表一支股票的開盤價(jià)。為了便于分析，可以先從中減去躍變，得到，即 （32）輸入以下命令，得到的頻譜如圖19所示。o=csvread('yhoodata.csv',2700,1,2700 1 2801 1)'N=102;for n=1:N x(n)=o(n)-o(1)-(o(N)-o(1)/(N-1)*(n

16、-1);endX=fft(x);k=0:N-1;stem(k,abs(X);axis(0 101 0 300);xlabel('k');ylabel('|Xk|')圖19 xn的幅度譜可以看出上圖中有5個(gè)較強(qiáng)譜分量，頻率（）分別對(duì)應(yīng)和。保留這5個(gè)頻率分量的系數(shù)，將其他頻率分量的系數(shù)置為0，然后再進(jìn)行逆傅里葉變換，得到濾波后的近似值。輸入如下Matlab程序，得到真實(shí)值與濾波后的近似值，如圖20所示。plot(x);hold on;fliterX=X(1:2) 0 0 X(5) zeros(1,102-9) X(99) 0 0 X(102);fliterx=iff

17、t(fliterX);plot(real(fliterx),'r');axis(1 102 0 7);xlabel('n');ylabel('xn的值和它的近似值');legend('xn真實(shí)值','xn近似值')圖20 xn的真實(shí)值與濾波后的近似值從上圖可以看出，濾波后的近似值既大致上保留了真實(shí)值的變化趨勢(shì)，而且與其十分接近。與濾波前比較，濾波后的圖形要比濾波前平滑得多。再由式（32）即可求得 （33）輸入如下Matlab程序，畫出真實(shí)開盤價(jià)與近似開盤價(jià)的圖形。如圖21所示，可以看出是近似基礎(chǔ)上的平滑值。plot

18、(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真實(shí)值','on近似值')圖21 on的真實(shí)值與濾波后的近似值5 利用FFT提取離散信號(hào)中的最強(qiáng)正弦分量這里最強(qiáng)是指在信號(hào)中某個(gè)正弦分量的振幅遠(yuǎn)大于其它正弦分量的振幅?？梢詫?duì)求點(diǎn)來確定信號(hào)中是否有最強(qiáng)的正弦分量。如果信號(hào)是連

19、續(xù)時(shí)間形式的，首先還要對(duì)其進(jìn)行抽樣，得到離散時(shí)間形式的信號(hào)，根據(jù)Nyquist定理，抽樣間隔應(yīng)滿，其中是中的最大頻率分量。要判斷信號(hào)中是否包含最強(qiáng)正弦分量，采樣數(shù)據(jù)至少要包含該分量一個(gè)完整周期的數(shù)據(jù)，例5 太陽耀斑數(shù)據(jù)的分析太陽耀斑活動(dòng)的周期是11年，這個(gè)事實(shí)可以通過提取耀斑數(shù)據(jù)的最強(qiáng)正弦分量加以證實(shí)。耀斑數(shù)據(jù)可以從比利時(shí)皇家天文臺(tái)(Royal Observatory of Belgium)太陽耀斑數(shù)據(jù)索引中心（Sunspot Index Data Center, SIDC）網(wǎng)站下載。網(wǎng)址是：http:/sidc.oma.be/index.php3下載后的數(shù)據(jù)存放在文件“monthssn.da

20、t”中，里面有四列數(shù)據(jù)，第一年是日期，第三列是太陽耀斑的平均數(shù)，第四列平滑后太陽耀斑的平均數(shù)，可以得到從1749年到當(dāng)前年月（2006年4月）的耀斑數(shù)據(jù)。本次分析選擇第三列1981年1月作為開始日期，2005年12月作為結(jié)束日期，共25年300個(gè)月份的數(shù)據(jù)。為此先把相關(guān)數(shù)據(jù)復(fù)制到Excel表格的第一列中，然后保存到Matlab所在目錄下，并命名為“sunspotdata.csv”。然后輸入以下命令，得到耀斑曲線如圖22所示。spd=csvread('sunspotdata.csv',0,0,0 0 299 0);plot(spd);grid;xlabel('月數(shù)'

21、;);ylabel('耀斑平均數(shù)');axis(0 300 0 200)圖22 1981年1月至2005年12月太陽耀斑的平均數(shù)由上圖可見，太陽耀斑的活動(dòng)確實(shí)具有周期性，但周期的準(zhǔn)確值不明顯?？梢酝ㄟ^數(shù)第一個(gè)峰值和第二個(gè)峰值之間的月份來估計(jì)周期的值。查驗(yàn)表中的數(shù)據(jù)得，第一個(gè)峰值為200.3，出現(xiàn)在第116個(gè)月（1990年8月），第二個(gè)峰值為170.1，出現(xiàn)在第235個(gè)月（2000年7月），所以周期是235-116=119月，和實(shí)際值132月比較接近。下面利用來分析。首先，從圖22中可以看出從到整個(gè)區(qū)間的平均值不為0，為了更方便地分析，需要減去該均值，得到結(jié)果如下： （33）然后對(duì)進(jìn)行傅里葉變換，得到。在Matlab中輸入以下命令，得到的幅度譜的圖形如圖23所示。x=spd-sum(spd)/300;X=fft(x);k=0:299;plot(k,abs(