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

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

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

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

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

9、nction wc=compress(w,r)%壓縮函數(shù)compress.m%輸入信號數(shù)據(jù)w,壓縮率r%輸出壓縮后的信號數(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做離散信號壓縮%輸入時間t,原信號y,以及壓縮率r%輸出原信號和壓縮后重構(gòu)信號的圖像,以及重構(gòu)均方差

10、和相對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('原信號','重構(gòu)信號');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小波分解信號，對信號進(jìn)行小波壓縮，然后重構(gòu)信號。令絕對值最小的80%系數(shù)為0，得到重構(gòu)信號圖形如圖18 a)所示，均方差為0.0584，相對誤差為0.0611；令絕對值最小的90%系數(shù)為0，得到重構(gòu)信號圖形如圖18 b)所示，均方差為0.1136，相對誤差為0.1190。 a) 絕對值最小的80%系數(shù)為0的重構(gòu)信號（Harr） b) 絕對值最小的90%系數(shù)為

12、0的重構(gòu)信號（Harr）圖18 用Harr小波壓縮后的重構(gòu)信號相關(guān)Matlab程序如下function unbiased_variance,error=daubcomp(t,y,n,r)%利用Daubechies系列小波做離散信號壓縮%輸入時間t,原信號y,分解層數(shù)n,以及壓縮率r%輸出原信號和壓縮后重構(gòu)信號的圖像,以及重構(gòu)均方差和相對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('原信號','重構(gòu)信號');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é)論：在信號沒有突變、快變化或者大致上具有周期性的信號，用FFT可以處理得很好（甚至比小波還要好）。4 利用FFT對離散信號進(jìn)行濾波利用FFT算法對信號進(jìn)行濾波的步驟如下：1）通過采樣將信號離散化；2）對離散化信號進(jìn)行傅里葉變換；3）對變換后的系數(shù)進(jìn)行處理，將絕對值大于某一閾值的系數(shù)置為0，保留剩余的系數(shù)；4）利用IFFT算法對處理后的信號進(jìn)行逆傅里葉變換。例4 股票價格分析首先進(jìn)入網(wǎng)址 To Spreadsheet按鈕，即可把以Excel表格格式存儲的價格數(shù)據(jù)下載到本地計算機(jī)。表格從1列至第6列分別給出了從1996年4月12日至2007年5月30日的交易期里每天的開盤價、

15、最高價、最低價、收盤價、成交量以及趨勢。數(shù)據(jù)下載完成后，需要顛倒順序，使得最早時間的數(shù)據(jù)首先顯示。然后另存到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)-(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個較強(qiáng)譜分量，頻率（）分別對應(yīng)和。保留這5個頻率分量的系數(shù)，將其他頻率分量的系數(shù)置為0，然后再進(jìn)行逆傅里葉變換，得到濾波后的近似值。輸入如下Matlab程序，得到真實值與濾波后的近似值，如圖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真實值','xn近似值')圖20 xn的真實值與濾波后的近似值從上圖可以看出，濾波后的近似值既大致上保留了真實值的變化趨勢，而且與其十分接近。與濾波前比較，濾波后的圖形要比濾波前平滑得多。再由式（32）即可求得 （33）輸入如下Matlab程序，畫出真實開盤價與近似開盤價的圖形。如圖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真實值','on近似值')圖21 on的真實值與濾波后的近似值5 利用FFT提取離散信號中的最強(qiáng)正弦分量這里最強(qiáng)是指在信號中某個正弦分量的振幅遠(yuǎn)大于其它正弦分量的振幅?？梢詫η簏c來確定信號中是否有最強(qiáng)的正弦分量。如果信號是連

19、續(xù)時間形式的，首先還要對其進(jìn)行抽樣，得到離散時間形式的信號，根據(jù)Nyquist定理，抽樣間隔應(yīng)滿，其中是中的最大頻率分量。要判斷信號中是否包含最強(qiáng)正弦分量，采樣數(shù)據(jù)至少要包含該分量一個完整周期的數(shù)據(jù)，例5 太陽耀斑數(shù)據(jù)的分析太陽耀斑活動的周期是11年，這個事實可以通過提取耀斑數(shù)據(jù)的最強(qiáng)正弦分量加以證實。耀斑數(shù)據(jù)可以從比利時皇家天文臺(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個月份的數(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ù)由上圖可見，太陽耀斑的活動確實具有周期性，但周期的準(zhǔn)確值不明顯?？梢酝ㄟ^數(shù)第一個峰值和第二個峰值之間的月份來估計周期的值。查驗表中的數(shù)據(jù)得，第一個峰值為200.3，出現(xiàn)在第116個月（1990年8月），第二個峰值為170.1，出現(xiàn)在第235個月（2000年7月），所以周期是235-116=119月，和實際值132月比較接近。下面利用來分析。首先，從圖22中可以看出從到整個區(qū)間的平均值不為0，為了更方便地分析，需要減去該均值，得到結(jié)果如下： （33）然后對進(jìn)行傅里葉變換，得到。在Matlab中輸入以下命令，得到的幅度譜的圖形如圖23所示。x=spd-sum(spd)/300;X=fft(x);k=0:299;plot(k,abs(