【正文】 xny ?????? （ 318） 可見 )(nyk 是由兩個序列卷積而得到的序列。 （ 320） 對比式 (317）和式 (318），可知：按式（ 319）進行卷積運算，當 Nn? 時，濾波器的輸出 )(Nyk 就是 )(kX ： Nnk nykX ?? |)()( （ 321） 對式 (320)進行 Z變換，可得濾波器的系統(tǒng)函數(shù) 11 1)( ???? zWzH kk （ 322） 這是一個一階系統(tǒng)。 按照式 (322)進行運算時，可先算好旋轉(zhuǎn)因子 kNW? ，儲存起來。按式 (316)直接計算 N 點離散付里葉變換 ,需要 24N 次實數(shù)乘法和)NN( 24 次實數(shù)加法。所以當 N 不大時，上述算法的效率稍差 [10]。 圖 32 用一階系統(tǒng)實現(xiàn) Goertzel 算法 圖 33 用二階系統(tǒng)實現(xiàn) Goert算法 把式（ 322）的分子 和分母都乘上因子 )1( 1?? zWkN ，就得到第 k 個濾波器的系統(tǒng)函數(shù)為 )1)(1( 1)( 111?????? ?? zWzW zWzH kNkNkNk 武漢工程大學(xué)畢業(yè)設(shè)計 (論文 )說明書 12 211 ))/2c o s ( (21 1??? ??? ?? zzkN zW kN （ 324） 與此相應(yīng)的信號流圖示于圖 33。 為了便于運算，在 圖 33 所示的流圖中，設(shè)立狀態(tài)變量 u 和 v 。每次運算中，更新狀態(tài)變量 u 和 v 。 。在 0)( ?Nx 點上。這時，按照圖 32 算出濾波器的輸出 )(Nyk ，此即 )(kX 。 綜上所述，計算一點 )(kX 需要進行 )2(2 ?N 次實數(shù)乘法和 )1(4 ?N 次實數(shù)加法。在這種較為有效的方案中，仍具有這樣的優(yōu)點，即必須計算和存儲的系數(shù)只有 ))/2cos(( kN? 和 kNW 。當輸入序列為實序列時，離散付里葉變換序列 )(kX 是對稱的，即 )()( * kNXkX ?? 。由于零點僅在最后的迭代中實現(xiàn)，所以諸極點要求的N2 次乘法和 N4 次加法可以用來計算離散付里葉變換的兩個值。然而，它同直接計算離散付里葉變換一樣，計算量仍然正比于 2N 。當 , 0 ,1, 2 , , 1k k N? ? ? ? ?時，利用 FFT算法可計算 ()X? 。最后令 2/NT??? ，則上式變?yōu)? 2/1( ) [ ] 0 , 1 , 2 , , 12/ i k NeX k X k k Ni k N T?? ??? ? ? ? （ 45） 首先用 FFT 算法求出 []Xk ，然后可用上式求出 0,1, 2, , 1kN??時的武漢工程大學(xué)畢業(yè)設(shè)計 (論文 )說明書 14 ()Xk? 。通過取更小的抽樣間隔 T ，或者增加點數(shù) N ，可以得到更精確的值。如果已知信號只在時間區(qū)間 10 tt?? 內(nèi)存在，可以通過對1nT t? 時的抽樣信號 [ ] ( )x n x nT? 補零，使 N 足夠大[12]。很多信號都能用（ 41）式連續(xù)時間傅立葉變換（ CTFT）來表示。利用在密集的等間隔 t的樣本上的求和來近似這個積分，就可以用函數(shù) fft高效地計算這個近似值。若信號 )(tx 對于 0?t 和 Tt? 為零，那么這個近似式就能寫成 ??? ???????? ?? 100 )()()(NnnjT tjtj enxdtetxdtetx ?? ???? （ 47） 式中 ?nT? ， N 為一整數(shù)。如果 N個樣本 )(?nx 是存在向量 x 內(nèi)的話，那么調(diào)用函數(shù) X=tau*fft(x)就可以計算出 )()()1( 10 kNnnj jXenxkX k ??? ?? ??? ???? （ 48） 式中 ?????????????? 12 2220 2NkNN kNkNkk??????? 以及 N假設(shè)為偶數(shù)。為了將頻率樣本置于上升的順序，能用函數(shù) fftshift。 例 41 利用 FFT計算傅里葉變換 如圖 41所示的信號 1 0 2() 0ttxt ? ? ??? ?? 其 它 其傅里葉變換為： 2c o s( ) sin ( )( ) 2 iX e i?? ? ?? ? ??? 利用下面的命令，可得到 ()xt 的近似值和準確值。Input N:39。 T=input(39。)。 x=[t1 zeros(1,Nlength(t))]。 gamma=2*pi/(N*T)。 Xapp=(1exp(i*k*gamma*T))/(i*k*gamma)*X。 Xact=exp(i*w)*2*i.*(w.*cos(w)sin(w))./(w.*w)。o39。 legend(39。,39。)。頻率 (rad/s)39。ylabel(39。) 運行程序后輸入 N=128， T=，此時 ?? ，得到實際的和近似的傅里葉變換的幅度譜如圖 42 所示，此時近似值已經(jīng)相當準確。再次運行程序后輸入 N=512， T=，此時 ?? ，得到實際的和近似的傅里葉變換的幅度譜如圖 43 所示。 例 42 利用 FFT 計算線性卷積 已知 [ ] [ ]nx n u n? ，其中 []un為單位階躍序列，信號 []yn 如圖 44所示。 利用下面的 Matlab 命令，可得到 []xn 、 []yn 的卷積圖形如圖 44所示。 n=0:16。 stem(n,x)。n39。ylabel(39。)。 n=0:15。 stem(n,y)。n39。ylabel(39。) subplot(3,1,3)。n=0:L1。Y=fft(y,L)。z=ifft(Z,L)。 xlabel(39。)。z[n]39。 subplot(2,2,1)。k=0:L1。x=.^n。 stem(k,abs(X))。 xlabel(39。)。|X[k]|39。 stem(k,angle(X))。 xlabel(39。)。Angle(X[k])(弧度 )39。 y=ones(1,10)。 武漢工程大學(xué)畢業(yè)設(shè)計 (論文 )說明書 19 stem(k,abs(Y))。 xlabel(39。)。|Y[k]|39。 stem(k,angle(Y))。 xlabel(39。)。Angle(Y[k])(弧度 )39。 例 43 對單位區(qū)間上的下列連續(xù)信號 1( ) c o s ( 4 ) s in ( 8 )2f t t t t??? ? ? 以 256sf Hz? 采樣頻率進行采樣，將其離散化為 82 個采樣值 [ ] ( ) | ( ) ( / 2 5 6 ) ,t n Tf n f t f n T f n?? ? ?0,1, 2, , 255n ? 用 FFT 分解信號，對信號進行小波壓縮，然后重構(gòu)信號。 a) 絕對值最小的 80%系數(shù)為 0 的重構(gòu)信號（ FFT） b) 絕對值最小的 90%系數(shù)為 0 的重構(gòu)信號（ FFT） 圖 46 用 FFT壓縮后的重構(gòu)信號 相關(guān) Matlab 程序如下 : function wc=press(w,r) %壓縮函數(shù) %輸入信號數(shù)據(jù) w,壓縮率 r %輸出壓縮后的信號數(shù)據(jù) if(r0)|(r1) error(39。)。 N=length(w)。 ww=sort(abs(w))。 wc=(abs(w)=tol).*w。r 應(yīng)該介于 0和 1之間 !39。 武漢工程大學(xué)畢業(yè)設(shè)計 (論文 )說明書 21 end。 fyc=press(fy,r)。 plot(t,y,39。,t,yc,39。)。原信號 39。重構(gòu)信號 39。 unbiased_variance=norm(yyc)/sqrt(length(t))。 輸入以下 Matlab 命 令： t=(0:255)/256。 [unbiased_variance,error]=fftp(t,f,) unbiased_variance = error = 如果用 Harr 尺度函數(shù)和 Harr 小波分解信號，對信號進行小波壓縮，然后重構(gòu)信號。 a) 絕對值最小的 80%系數(shù)為 0 的重構(gòu)信號（ Harr） b) 絕對值最小的 90%系數(shù)為 0 的重構(gòu)信號（ Harr） 圖 47 用 Harr小波壓縮后的重構(gòu)信號 武漢工程大學(xué)畢業(yè)設(shè)計 (論文 )說明書 22 相關(guān) Matlab 程序如下 function [unbiased_variance,error]=daubp(t,y,n,r) %利用 Daubechies 系列小波做離散信號壓縮 %輸入時間 t,原信號 y,分解層數(shù) n,以及壓縮率 r %輸出原信號和壓縮后重構(gòu)信號的 圖像 ,以及重構(gòu)均方差和相對 l^2誤差 if(r0)|(r1) error(39。)。 [c,l]=wavedec(y,n,39。)。 %調(diào)用壓縮函數(shù) yc=waverec(cc,l,39。)。r39。b39。 legend(39。,39。)。 error=norm(yyc)/norm(y)。 f=t+cos(4*pi*t)+1/2*sin(8*pi*t)。 4. 4 利用 FFT對離散信號進行濾波 利用 FFT算法對信號進行濾波的步驟如下： 1）通過采樣將信號離散化； 2）對離散化信號進行傅里葉變換； 3）對變換后的系數(shù)進行處理，將絕對值大于某一閾值的系數(shù)置為 0，保留剩余的系數(shù)； 4）利用 IFFT 算法對處理后的信號進行武漢工程大學(xué)畢業(yè)設(shè)計 (論文 )說明書 23 逆傅里葉變換 [14]。要求：將信號的幅度換算成實際的幅度，信號的頻率換算成實際的頻率。%采樣頻率 N=512。 randn(1,N)。%采樣時間序列 f0=100。 subplot(3,1,1)。 xlabel(39。)。sin(2+3*cos(100*pi*tpi/6)+(3/2)*cos(150*pi*t+pi/2))39。 title(39。)。%對信號進行 FFT 變換 magY=abs(Y)。%求得 FFT 變換后的相位 f=n*fs/N。 plot(f,magY)。f39。 ylabel(39。)。N=51239。 grid。 plot(f,angY)。f39。 ylabel(39。) title(39。)。 圖 48 濾波信號的相位，幅度 利用 FFT提 取離散信號中的最強正弦分量 這里最強是指在信號 [ ], 0 ,1, , 1x n n N??中某個正弦分量的振幅遠大于其它正弦分量的振幅。如果信號 ()xt 是連續(xù)時間形式的，首先還要對其進行抽樣，得到離散時間形式的信號 [ ] ( ) | ( )t n Tx n x t x n T???，根據(jù) Nyquist定理，抽樣間隔 T 應(yīng)滿 max/T ??? ，其中 max? 是 ()xt 中的最大頻率分量。 例 45 太陽耀斑數(shù)據(jù)的分析 太陽耀斑活動的周期是 11 年，這個事實可以 通過提取耀斑數(shù)據(jù)的最強正弦 分量加以證實。網(wǎng)址是： 武漢工程大學(xué)畢業(yè)設(shè)計 (論文 )說明書 25 下載后的數(shù)據(jù)存放在文件“ ”中，里面有四列數(shù)據(jù)，第一年是日期，第三列是太陽耀斑的平均數(shù)，第四列平滑后太陽耀斑的平均數(shù)，可以得到從 1749 年到當前年月（ 2020 年 4 月）的耀斑數(shù)據(jù)。為此先把相關(guān)數(shù)據(jù)復(fù)制到 Excel 表格的第一列中，然后保存到 Matlab 所在目錄下，并命名為“ ”。 spd=csvread(39。,0,0,[0 0 299 0])。 grid。月數(shù) 39。ylabel(39。)?？梢酝ㄟ^數(shù)第一個峰值和第二個峰值之間的月份來估計周期的值。 武漢工程大學(xué)畢業(yè)設(shè)計 (論文 )說明書