【正文】 查驗表中的數(shù)據(jù)得，第一個峰值為 ，出現(xiàn)在第 116 個月（ 1990 年 8 月），第二個峰值為 ，出現(xiàn)在第 235個月（ 2020 年 7月），所以周期是 235116=119 月，和實際值 132月比較接近。 axis([0 300 0 200]) 圖 49 1981年 1月至 2020年 12月太陽耀斑的平均數(shù) 由上圖可見，太陽耀斑的活動確實具有周期性，但周期的準確值不明顯。耀斑平均數(shù) 39。)。 xlabel(39。 plot(spd)。39。然后輸入以下命令，得到耀斑曲線如圖 9所示 [16]。本次分析選擇第 三列 1981年 1 月作為開始日期， 2020 年 12 月作為結(jié)束日期，共 25 年 300 個月份的數(shù)據(jù)。耀斑數(shù)據(jù)可以從比利時皇家天文臺 (Royal Observatory of Belgium)太陽耀斑數(shù)據(jù)索引中心（ Sunspot Index Data Center, SIDC）網(wǎng)站下載。要判斷信號 []xn 中是否包含最強正弦分量，采樣數(shù)據(jù)至少要包含該分量一個完整周期的數(shù)據(jù) [15]?？梢詫?[]xn 求 N 點 DFT 來確定信號中是否有最強的正弦分量。 grid。N=51239。相位 39。)。%畫出相頻響應(yīng) 武漢工程大學(xué)畢業(yè)設(shè)計 (論文 )說明書 24 xlabel(39。 subplot(3,1,3)。)。 title(39。幅度 39。)。%畫出幅頻響應(yīng) xlabel(39。%頻率序列 subplot(3,1,2)。%求得 FFT 變換后的幅度 angY=angle(Y)*180/pi。 Y=fft(x,N)。時域信號 39。)。 ylabel(39。t39。 plot(t,x)。%信號頻率 x=2+3*cos(100*pi*tpi/6)+(3/2)*cos(150*pi*t+pi/2)。 t=0:1/fs:(N1)/fs。%數(shù)據(jù)點數(shù) n=0:N1。 相關(guān) Matlab 程序如下 fs=1000。 例 44 對被白噪聲污染的信 ? ? ? ?2/150c o ???? ????? tt進行頻譜分析，從中鑒別出有用的信號。 [unbiased_variance,error]=daubp(t,f,8,) unbiased_variance = error = 結(jié)論：在信號沒有突變、快變化或者大致上具有周期性的信號，用 FFT 可以處理得很好（甚至比小波還要好）。 輸入以下 Matlab 命令： =(0:255)/256。 unbiased_variance=norm(yyc)/sqrt(length(t))。重構(gòu)信號 39。原信號 39。)。,t,yc,39。 plot(t,y,39。db139。 cc=press(c,r)。db139。 end。r 應(yīng)該介于 0和 1 之間 !39。令絕對值最小的 80%系數(shù)為 0，得到重構(gòu)信號圖形如圖 47a)所示，均方差為 ，相對誤差為 ；令絕對值最小的 90%系數(shù)為 0，得到重構(gòu) 信號圖形如圖 47b)所示，均方差為 ，相對誤差為 。 f=t+cos(4*pi*t)+1/2*sin(8*pi*t)。 error=norm(yyc)/norm(y)。)。,39。 legend(39。b39。r39。 %調(diào)用壓縮函數(shù) yc=ifft(fyc)。 fy=fft(y)。)。 function [unbiased_variance,error]=fftp(t,y,r) %利用 FFT 做離散信號壓縮 %輸入 時間 t,原信號 y,以及壓縮率 r if(r0)|(r1) error(39。 tol=abs(ww(Nr+1))。 Nr=floor(N*r)。 end。r 應(yīng)該介于 0和 1之間 !39。令絕對值最小的80%系數(shù)為 0，得到重構(gòu)信號圖形如圖 46a)所示，均方差為 ，相對誤差武漢工程大學(xué)畢業(yè)設(shè)計 (論文 )說明書 20 為 ；令絕對值最小的 90%系數(shù)為 0，得到重構(gòu)信號圖形如圖 46b)所示，均方差為 ，相對誤差為 。) 圖 45 信號 x[n]、 y[n]的幅度譜與相位譜 FFT 進行離散信號壓縮 利用 FFT算法對離散信號進行壓縮的步驟如下： 1）通過采樣將信號離散化；2）對離散 化信號進行傅里葉變換； 3）對變換后的系數(shù)進行處理，將絕對值小于某一閾值的系數(shù)置為 0，保留剩余的系數(shù)； 4）利用 IFFT 算法對處理后的信號進行逆傅里葉變換 [13]。ylabel(39。k39。 axis([0 25 3 3])。) subplot(2,2,4)。ylabel(39。k39。 axis([0 25 0 10])。Y=fft(y,L)。) subplot(2,2,3)。ylabel(39。k39。 axis([0 25 1 1])。) subplot(2,2,2)。ylabel(39。k39。 axis([0 25 0 5])。X=fft(x,L)。 n=0:16。 L=26。) 圖 44 信號 x[n]、 y[n]及其卷積 z[n]=x[n]*y[n] 利用下面的 Matlab 命令，可得到信號 x[n]、 y[n]的幅度譜與相位譜如圖 45所示。ylabel(39。n39。 武漢工程大學(xué)畢業(yè)設(shè)計 (論文 )說明書 18 stem(n,z)。 Z=X.*Y。 X=fft(x,L)。 L=26。y[n]39。)。 xlabel(39。 y=[ones(1,10) zeros(1,6)]。 subplot(3,1,2)。x[n]39。)。 xlabel(39。 x=.^n。 subplot(3,1,1)。由于當 16n? 時， []xn 很小，故 M 可以取為 17； N取 10， 1 26L M N? ? ? ?。 武漢工程大學(xué)畢業(yè)設(shè)計 (論文 )說明書 16 圖 42 N=128， T= 圖 43 N=512， T= 武漢工程大學(xué)畢業(yè)設(shè)計 (論文 )說明書 17 FFT 計算離散信號的線性卷積 已知兩個離散時間信號 [ ] ( 0 ,1, 2 , 1)x n n M??與 [ ] ( 0 ,1, 2 , 1)y n n N??，取L? 1MN??，對 []xn 和 []yn 右端補零，使得 [ ] 0 , 1 , 2 , , 1x n n M M L? ? ? ? ? [ ] 0 , 1 , 2 , , 1y n n N N L? ? ? ? ? （ 49） 利用 FFT算法可以求得 []xn 和 []yn 的 L 點 DFT，分別是 []Xk 和 []Yk，利用DTFT 卷積性質(zhì)，卷積 [ ]* [ ]x n y n 等于乘積 [ ] [ ]X k Y k 的 L 點 DFT 反變換，這也可以通過 FFT 算法得到。通過增加 NT 可以增加更多的細節(jié)，減少 T使得到的值更精確。|X|39。)。 xlabel(39。真實值 39。近似值 39。,w,abs(Xact))。 plot(k*gamma,abs(Xapp(1:length(k))),39。 %計算真實值 X(w) w=::10。 k=0:10/gamma。 X=fft(x)。 %計算 X(w)近似值 t=0:T:2。Input T:39。)。 圖 41 連續(xù)時間信號 x(t) N=input(39。為了將存入 X 中的 )( kjX ? 的樣本排列成使 )1( ?KX 就是對于 10 ??? Nk ，在 )2( ???? Nk?? 上求得武漢工程大學(xué)畢業(yè)設(shè)計 (論文 )說明書 15 的 CTFT,可用 X=fftshift(tau*fft(x))。為了計算高效， fft 在負的頻率樣本之前先產(chǎn)生正頻率樣本?？梢岳煤瘮?shù) fft 對一組離散的頻率 k? 計算上式中的和式。所用的近似式是根據(jù)積分的定義得到的，即 ?? ????????? ? ? n njtj enxdtetx ?? ???? )(lim)( 0 （ 46） 對于一般信號，在足夠小的τ下，上式右邊的和式是對于 CTFT 積分的一個好的近似。利用 MATLAB 可以計算（ CTFT）積分的數(shù)值近似。 將連續(xù)時間傅立葉變換進行數(shù)字近似，用函數(shù) fft（快速傅立葉算法）高效地計算這個近似值。如果B?? 時，幅度譜 ()X? 很小，對應(yīng)于奈奎斯特抽樣頻率 2s B?? ，抽樣間隔 T 選擇 /B? 比較合適。 應(yīng)該強調(diào)的是，式（ 43）只是一個近似表示，計算得到的 ()X? 只是一個近似值。 已知一個固定的時間間隔 T ，選擇讓 T 足夠小，使得每一個 T秒的間隔( 1)nT t n T? ? ? 內(nèi)， ()xt 的變化很小，則式中積分可近似為 ( 1 )0( ) ( ) ( )nT i w tnTnX e d t x n T?? ? ??? ? ? ( 1 )01[ ] ( )i t t n Tt n Tn e x n Ti??? ? ? ????? ? 01 ()iT i n Tne e x n Ti ? ??? ? ???? ? （ 42） 假設(shè) N 足夠大，對于所有 nN? 的整數(shù)，幅值 ()xnT 很小，則式（ 42）變?yōu)? 101( ) ( )iT N i n TneX e x n Ti ? ?? ?? ? ???? ? （ 43） 當 2/k NT??? 時，式（ 42）兩邊的值為 2 / 2 /1 2/02 1 1( ) ( ) [ ]2 / 2 /i k N i k NN i n k Nnk e eX e x n T X kN T i k N T i k N T???? ????? ??????? （ 44） 其中 []Xk 代表抽樣信號 [ ] ( )x n x nT? 的 N 點 DFT 。 武漢工程大學(xué)畢業(yè)設(shè)計 (論文 )說明書 13 在信號處理 中 的理論應(yīng)用 連續(xù)時間信號的快速傅里葉變換 設(shè) ()xt 是連續(xù)時間信號，并假設(shè) 0t? 時 () 0xt? ，則其傅里葉變換由下式給出： 0( ) ( ) itX x t e d t?? ? ?? ? （ 41） 令 ? 是一固定的正實數(shù)， N 是一個固定的正整數(shù)。因此，若用Goertzel 算法計算離散付里葉變換的所有 N 個點的值，需要的乘法次數(shù)近似為2N ，加法次數(shù)近似為 22N 。容易證明，圖 33的網(wǎng)絡(luò)形式在計算 )( kNX ? 時和計算 )(kX 時具有完全相同的極點，但前者的零點系數(shù)與后者的零點系數(shù)成復(fù)共軛關(guān)系。 還要說明圖 33所示的算法的另一個優(yōu)點。這種算法要求的乘法次數(shù)約為直接算法的一半。所需的計算量是 4 次實數(shù)乘法和 N4 次實數(shù)加法 [11]。計算狀態(tài)變量 u 和 v 。 )(nx 是一個 N 點序列， 1,...,2,1,0 ?? Nn 。作 N 次迭代所需的計算量是 N2 次實數(shù)乘法和 N4 次實數(shù)加法。按照圖 33計算 )(kX 時，步驟有二，即： 輸入點依次取 )1() , . . . ,2(),1(),0( ?nxxxx ，進行遞推運算。由式（ 324)可見，濾波器是一個二階系統(tǒng)，有一對復(fù)數(shù)共軛極點和一個復(fù)數(shù)零點。下面介紹改進的 Goertzel 算法，這種算法所需的實數(shù)乘法次數(shù)約為直接方法的一半。按照上述 Goertzel 算法，所需的實數(shù)乘法和實數(shù)加法都是 24N 次。每次遞推包含一次復(fù)數(shù)乘法。圖 32示出這個系統(tǒng)的信號流圖，相應(yīng)的差分方程為 )()1()( nxnyWny kkNk ??? ? ， 0)1( ??y (323) 按照此式進行遞推運算，到了 Nn? 時刻，即可依據(jù)式（ 321）得到 )(kX 。 )()()( nhnxny