非平穩(wěn)信號實時譜分析算法及其FPGA實現(xiàn)

周 圍，梁 琦

(1.重慶郵電大學(xué) 移動通信技術(shù)重慶市重點實驗室，重慶 400065; 2.重慶郵電大學(xué) 光電工程學(xué)院，重慶 400065)

0 引 言

在工業(yè)監(jiān)控領(lǐng)域，通過對旋轉(zhuǎn)設(shè)備中與軸旋轉(zhuǎn)頻率相關(guān)的振動信號與噪聲信號功率譜特征進行分析來判斷旋轉(zhuǎn)機電設(shè)備的運行情況并進行實時故障檢測是十分有效的方法[1]。在電力系統(tǒng)中，通過對電源信號進行譜分析能夠?qū)崟r測量由于大負(fù)載(如：電動機、焊接機)接入或開機運行產(chǎn)生的電壓失真、抖動和不平衡[2]。在頻譜監(jiān)控領(lǐng)域，如感知無線電中通過動態(tài)頻譜接入和快速頻譜感測技術(shù)實時檢測用戶是否存在[3]。因此，譜分析是大多數(shù)學(xué)科中的一個重要任務(wù)，且其主要目的是功率譜測量。傳統(tǒng)的譜分析方法一般采用基于快速傅里葉變換(fast Fourier transform,FFT)及其改進算法和頻譜細(xì)化分析法。FFT算法得到的是整個頻帶的粗略全景譜，頻率細(xì)化算法可以得到頻譜的局部詳細(xì)特征。但由于上述方法只能分析統(tǒng)計平均結(jié)果，因此均適用于平穩(wěn)信號的譜分析，而無法處理非平穩(wěn)信號。實際觀察到的信號(如振動信號)大都為非平穩(wěn)信號，對其非平穩(wěn)性的分析可以實現(xiàn)故障的實時診斷、排除和頻譜實時監(jiān)控。非平穩(wěn)信號典型的處理方法有短時傅里葉變換、Wigner-Ville分布、Gabor變換、小波分析和分?jǐn)?shù)階傅里葉變換?？紤]硬件實現(xiàn)復(fù)雜性和實時性，工程上一般采用短時傅里葉變換。為克服短時傅里葉變換的頻率分辨率完全受限于滑動窗口長度和其無法分析局部頻譜的缺陷。本文擬采用窗口長度逼近平穩(wěn)信號的觀察窗采樣非平穩(wěn)信號，然后利用多相濾波器得到分析頻帶重構(gòu)信號，之后通過加窗線性調(diào)頻Z變化(chirp Z transform,CZT)實現(xiàn)與時間分辨率無關(guān)的頻譜細(xì)化，最后利用累加平均與譜分析解決旁瓣誤差和噪聲誤差，并利用FPGA(field-programmable gate array)硬件實現(xiàn)。

1 輸入輸出描述

圖1為本設(shè)計的全局視圖。圖1中，fmin和fmax分別為譜分析起、止頻率；Ns為分析帶寬內(nèi)最大信號數(shù)目；Δf為頻率分辨率和IDR(instantaneous dynamic range)(單位：dB)為瞬時動態(tài)范圍。假設(shè)輸入平穩(wěn)信號由R個信號組成，采樣窗口內(nèi)非平穩(wěn)信號由R個能量集中的窄帶信號組成，在此假設(shè)的基礎(chǔ)上每個窄帶信號可以作為一個復(fù)數(shù)基。則

(1)

圖1 頻譜分析儀全局視圖Fig.1 Global view of spectrum analyzer

對信號進行譜分析并將各時刻檢測到的信號數(shù)目、各信號中心頻率和功率封裝為信號描述字(signal description word,SDW)。通過SDW的運用將分析帶寬內(nèi)包含的譜信息簡化為一個描述字方便用戶應(yīng)用。

2 關(guān)鍵模塊理論分析

2.1 CZT算法原理

設(shè)通過ADC采樣信號為x(n)(n=0,…,N-1)。根據(jù)文獻[4]可知，其Z變換可表示為

(2)

(2)式中：f(n)=x(n)A-nWn2/2；h(n)=W-n2/2，其中，A=A0eiθ0,W=W0e-iφ0。A0為z平面上起始點幅值，θ0為起始點相位(對應(yīng)起始點頻率)，W0為幅值變化梯度倒數(shù)，φ0為相位增量(對應(yīng)頻率分辨率)。A0=W0=1，使CZT變換路徑為單位圓上的一段圓弧，令M=N(M為譜線數(shù))簡化CZT計算復(fù)雜度，通過θ0調(diào)整起始頻率，通過φ0確定選頻范圍內(nèi)細(xì)化倍數(shù)。

(2)式的計算可通過圖2所示步驟來實現(xiàn)。

圖2 CZT計算框圖Fig.2 Block diagram of CZT calculation

圖2中，f(n)與h(n)的離散卷積可通過其循環(huán)卷積實現(xiàn)。

2.2 窗函數(shù)選擇原理

窗函數(shù)的選擇是確保該算法成功應(yīng)用的最關(guān)鍵部分。系統(tǒng)歸一化的頻率分辨率可表示為

(3)

(3)式中，D為頻率細(xì)化倍數(shù)。則窗函數(shù)必須滿足下述條件。

1)窗函數(shù)旁瓣電平(side lobe level,SLL)需滿足SLL>IDR。否則將無法滿足用戶輸入IDR的要求。

2)窗函數(shù)6 dB帶寬需滿足：BW-6 dB<Δfb。在滿足Δf要求下，為正確區(qū)分2個同功率信號，采用文獻[5]中的典型標(biāo)準(zhǔn)BW-6 dB<Δfb。

3)窗口幅值下降IDR(dB)帶寬需滿足：BW-IDR<2Δfb。為保證Δf適用于更大的IDR，IDR對應(yīng)的窗口寬度必須比2Δfb更窄。

為準(zhǔn)確檢測Ns個信號，在此考慮最壞旁瓣條件(即輸入信號為Ns-1個同功率信號和一個滿足IDR的最小功率信號)，則此時窗函數(shù)應(yīng)滿足(4)式條件。

SLL>IDR+20lg(Ns-1)+SL

(4)

(4)式中：Ns>1；SL為窗口扇形衰落即當(dāng)信號頻率為2個CZT頻率點中間時的功率損耗[5]。

根據(jù)IDR的要求，歸一化頻率分辨率可表示為

(5)

(5)式中，BWnull為窗函數(shù)主瓣寬度。

對于非平穩(wěn)或瞬態(tài)信號的分析，窗函數(shù)除滿足上述頻域要求外還應(yīng)根據(jù)時域分辨率和最小檢測信號長度的要求選擇適合的采樣窗長度及重疊幀數(shù)。

2.3 譜分析算法原理

為實現(xiàn)降低測量誤差，得到更精確的測量結(jié)果，本文設(shè)計了3個閾值參數(shù)。該部分對各閾值參數(shù)設(shè)置理論進行了詳細(xì)闡述，具體參數(shù)示意如圖3所示。

圖3 閾值參數(shù)示意圖Fig.3 Schematic diagram of threshold parameter

1)窗口主瓣寬度nb：為了簡化硬件實現(xiàn)將加窗后信號主瓣假設(shè)為由奇數(shù)點CZT分辨率點組成，nb=2|BWnull/2|+1。避免了若nb為偶數(shù)時為降低錯誤檢測概率消除信號兩側(cè)nb/2點還是(nb-1)/2點。

2)噪聲閾值T(單位：dB)：為了避免由噪聲引起的測量誤差，在此引進噪聲閾值T(單位：dB)。根據(jù)文獻[6]定義

(6)

(7)

當(dāng)M?1時，將(7)式代入(6)式得

PFA|M?1≈e-γ

(8)

此時，PFA檢測性能趨于(6)式理論值。

若譜線數(shù)為MCZT點，且其服從獨立同分布，則

(9)

根據(jù)(8)，(9)式可得

γ≈-ln(PFA,g/MCZT)

(10)

3)動態(tài)閾值S(dB)：為保證IDR的要求，避免由于窗函數(shù)旁瓣造成的檢測誤差，設(shè)定動態(tài)閾值S=IDR+SL，其中，SL為窗函數(shù)的扇形衰落。

2.4 多相濾波原理

圖4為子帶設(shè)計示意圖。設(shè)系統(tǒng)采樣率為fs，將其劃分為N個等間隔子帶，則第i個子帶的中心頻率為ωi=2πfi/fs。設(shè)計任何相鄰子帶的阻帶截止頻率為對方的中心頻率，且子帶間中心頻率間隔f1滿足f1=mΔf(m為正整數(shù))[7]。

圖4 子帶劃分示意圖Fig.4 Diagram of subband division

原型濾波器傳輸函數(shù)為

(11)

(11)式中M為濾波器長度，設(shè)子帶i系數(shù)函數(shù)為

(12)

(12)式中WN=ej2π/N，則子帶i傳輸函數(shù)可表示為

(13)

3 FPGA設(shè)計實現(xiàn)

本文設(shè)計譜分析儀期望在算法與硬件實現(xiàn)上高度可配置，因此，設(shè)計時必須考慮多參數(shù)組合。在此以表1中2個實例介紹本算法內(nèi)部參數(shù)計算過程。

表1 設(shè)計實例用戶輸入?yún)?shù)

1)根據(jù)用戶輸入IDR=30 dB，Ns=4，由(4)式可知，SLL≥40。在此，窗函數(shù)選擇Hamming窗(根據(jù)文獻[5]可知，Hamming窗SSL=-43 dB，SL=1.78 dB和主瓣點數(shù)=4)以滿足IDR要求，且nb=5，S=33。根據(jù)起止頻率與系統(tǒng)采樣率得D=16，根據(jù)歸一化分辨率公式、頻率分辨率及(5)式得NCZT=M=512。由起始頻率與初始相位對應(yīng)關(guān)系和NCZT=512得θ0=2π×fmin/fs=2π·35/512。

2)同理可知，在此需選用Blackman-Harris窗，且nb=9，S=78，D=4，NCZT=M=512和θ0=2π·77/512。

圖5為譜分析儀系統(tǒng)結(jié)構(gòu)圖，由3個部分組成。硬件實現(xiàn)的關(guān)鍵問題是模塊間的跨時域傳輸和模塊內(nèi)的流水線處理，因此，本文通過系統(tǒng)延時分析確定各時鐘域時鐘頻率。下面將對譜分析儀各主要模塊的硬件實現(xiàn)進行詳細(xì)敘述。

3.1 濾波器組實現(xiàn)

為使濾波器滿足實時性和固定數(shù)目子帶劃分要求，原型濾波器(抗混疊濾波器)需選用運算量小且截止頻率與采樣率相關(guān)的濾波器。因此，本設(shè)計采用整系數(shù)升余弦濾波器實現(xiàn)。

首先利用Matlab設(shè)計實現(xiàn)升余弦濾波器并導(dǎo)出其系數(shù)，根據(jù)(12)式計算出各子帶濾波系數(shù)。然后將上述濾波系數(shù)存入相應(yīng)ROM(read only memory)中，并通過與采樣數(shù)據(jù)乘累加實現(xiàn)濾波功能，對濾波結(jié)果進行移位實現(xiàn)歸一化處理。為達到測量精度與實時性能的平衡，平穩(wěn)信號與非平穩(wěn)信號原型濾波器帶寬分別為fS/20和fS/50。濾波器組實現(xiàn)結(jié)構(gòu)如圖6所示。

圖5 系統(tǒng)硬件框圖Fig.5 Block diagram of system hardware

圖6 多相濾波器結(jié)構(gòu)圖Fig.6 Structure diagram of multiphase filter

3.2 CZT實現(xiàn)

根據(jù)圖2算法流程可知，該算法核心為f(n)和h(n)的卷積。由傅里葉變換卷積定理可知循環(huán)卷積可通過FFT快速計算[8]，計算流程如圖7所示。

根據(jù)CZT原理可知，產(chǎn)生信號與用戶輸入起始頻率、終止頻率、細(xì)化倍數(shù)和輸出譜線數(shù)有關(guān)，因此信號產(chǎn)生可通過2個DDS(direct digital synthesizer)模塊和一個數(shù)據(jù)緩存RAM實現(xiàn)。為提高系統(tǒng)運行速率將2個DDS模塊與2個FFT均采用并行架構(gòu)實現(xiàn)，為節(jié)約資源提高芯片利用率將復(fù)數(shù)乘法器復(fù)用并將FFT1復(fù)用至FFT3。此時FPGA實現(xiàn)CZT模塊數(shù)據(jù)流如圖8所示。

圖7 FFT實現(xiàn)循環(huán)卷積CZT流程Fig.7 FFT implementation of the convolution CZT process

圖8 CZT算法硬件實現(xiàn)數(shù)據(jù)流圖Fig.8 CZT algorithm hardware implementation data flow diagram

圖9 頻譜分析結(jié)構(gòu)框圖Fig.9 Structure diagram of spectrum analysis

3.3 譜分析實現(xiàn)

通過對CZT模塊輸出數(shù)據(jù)進行累加平均降低瞬態(tài)誤差，將累加平均結(jié)果通過頻譜分析得到SDW。具體實現(xiàn)步驟如圖9所示。

累加模塊在FPGA中通過2個移位寄存器與加法器實現(xiàn)多個功率譜對應(yīng)位置數(shù)據(jù)的累加。累加譜數(shù)量由用戶輸入?yún)?shù)決定，為提高系統(tǒng)運行速率將其設(shè)為2的整數(shù)次冪，則平均可通過移位實現(xiàn)。

創(chuàng)建2個長度為Ns的存儲器，分別存儲功率值和對應(yīng)數(shù)據(jù)位置。則譜分析可通過下述步驟實現(xiàn)。

2)噪聲估計。噪聲功率為不考慮與上述Ns個信號相關(guān)的輸出功率平均值：

(14)

(14)式中，Ms=M-nbNs，b[k]為空白窗,即

(15)

4)SDW產(chǎn)生。將分析帶寬內(nèi)信號數(shù)ND、信號檢測頻率和功率通過輸出端口輸出。

3.4 系統(tǒng)時延

圖10為譜分析儀數(shù)據(jù)流結(jié)構(gòu)。在此定義2個概念：塊延時lB(數(shù)據(jù)幀進入系統(tǒng)的時間間隔)和總延時lT(從第一個數(shù)據(jù)進入到SDW輸出時間間隔)。根據(jù)圖10可知，lT=lFIR+lCZT+lPA，因此,lB

(16)

圖10 譜分析儀時延示意圖Fig.10 Schematic diagram of spectrum analyzer time-delay

多相濾波模塊以采樣率fs連續(xù)處理數(shù)據(jù)并將結(jié)果存入緩存RAM，因此，該模塊的延時與塊延時相等即lFIR=lB。當(dāng)?shù)?個數(shù)據(jù)塊通過多相濾波模塊時，CZT模塊將開始處理數(shù)據(jù)，因此，CZT模塊必須在下一個數(shù)據(jù)塊完成多相濾波之前處理完該數(shù)據(jù)塊，即lCZT≤lFIR。CZT模塊的延時由10個周期地址運算時延、3個周期復(fù)數(shù)乘法器、2·NCZT個周期緩存時延和2·NCZT個周期IFFT時延組成，lCZT可表示為

(17)

(17)式中，fCZT/2由FFT的結(jié)構(gòu)決定。

在CZT模塊加載新的一組數(shù)據(jù)之前，譜分析模塊應(yīng)完成上一組數(shù)據(jù)處理，即lPA≤lFIR。譜分析首先需要將CZT模塊輸出數(shù)據(jù)加載到譜分析模塊，再對NCZT點CZT輸出結(jié)果中進行Ns次最大值搜索，一次噪聲估計和280個時鐘周期的其他處理(譜調(diào)整及SDW生成)。因此，譜分析延時lPA可表示為

(18)

通過上述各模塊延時分析可知lT≤3·lFIR。

4 系統(tǒng)驗證

在此部分將對設(shè)計譜分析儀在平穩(wěn)信號和非平穩(wěn)信號實時分析方面性能進行驗證，硬件平臺選用Altera Cyclone IVE:EP4CE115F29I7 FPGA。

4.1 平穩(wěn)信號性能驗證

為驗證本設(shè)計系統(tǒng)在不同用戶輸入?yún)?shù)組合下測量性能，在此仿真表1中2個實例進行系統(tǒng)性能驗證。在實例中選擇各緩沖延時周期為NCZT個樣點周期。

實例1對于高采樣率高頻率分辨局部譜分析。根據(jù)分析帶寬和濾波器子帶帶寬，在此實例中選擇子帶2，3和4，為滿足譜分析性能約束需對CZT模塊與PA模塊的時鐘進行適當(dāng)設(shè)置。

根據(jù)文獻[9]可知，F(xiàn)FT模塊最大頻率fFFT=400 ΜΗz，因此，fCZTmax=400 MHz。由(16)式可得，lFIR=lB=15.36 μs，由(17)式和CZT模塊時序約束可得fCZTmin=268 MHz，由(18)式和譜分析模塊時序約束可得fPAmin=212 MHz。為降低系統(tǒng)功耗并保證信號流實時處理，根據(jù)fCZTmax=400 MHz,可得fPAmax=250 MHz。以fCZT=400 MHz，fPA=250 MHz配置CZT和PA模塊，此時系統(tǒng)延時lT=38.27 μs。

設(shè)輸入平穩(wěn)信號x1(i)表達式為

(19)

(19)式中，f1=8 MHz；f2=9 MHz；f3=10 MHz；f4=12 MHz和fs=100 MHz。其Modelsim時序仿真所得波形圖如圖11。

由圖11所得檢測信號頻率值和功率值。輸入信號與檢測信號頻率和歸一化功率值如表2所示。

圖11 實例1時序仿真波形圖Fig.11 Timing simulation waveform of example 1

表2 實例1信號測量結(jié)果對比

實例2對于高采樣率較低頻率分辨局部譜分析。同理可知，選擇子帶為4，5，6，7，8，9，lFIR=15.36 μs，fCZTmin=268 MHz，fPAmin=174 MHz。滿足要求的時鐘頻率配置如表3所示。在此選用fCZT=400 MHz和fPA=250 MHz配置CZT及PA模塊時鐘，此時系統(tǒng)總延時lT=lFIR+lCZT+lPA=36.225 μs。

表3 實例2滿足連續(xù)實時處理時鐘頻率

設(shè)輸入平穩(wěn)信號x2(i)表達式為

(20)

(20)式中：f1=35 MHz；f2=36 MHz；f3=60 MHz；f4=72 MHz和fs=200 MHz，其Modelsim時序仿真所得波形圖如圖12所示。

由圖12所得檢測信號頻率值和功率值。輸入信號與測量信號頻率和歸一化功率值如表4所示。

圖12 實例2時序仿真波形圖Fig.12 Timing simulation waveform of example 2

表4 實例2信號測量結(jié)果對比

4.2 非平穩(wěn)信號性能驗證

為驗證本文實現(xiàn)系統(tǒng)具有非平穩(wěn)信號譜分析性能，對其進行測試。設(shè)一個非平穩(wěn)信號f(i)，分析帶寬為0～30 MHz，每個子帶最多能量集中點為3，f(i)表達式為

(21)

(21)式中，信號頻率為f1=1 MHz，f2=9 MHz，f3=20 MHz和fs=250 MHz。為保證頻率分析精度和實時性要求，在此采用8倍細(xì)化512點CZT，每CZT周期窗口時移樣點數(shù)為52，各緩沖器延時為NCZT個樣點周期。

根據(jù)分析帶寬和濾波器子帶帶寬，在此選擇子帶1，2，3，4，5和6。lFIR=12.288 μs，fCZTmin=335 MHz和fPAmin=169 MHz。滿足要求的時鐘配置如表5所示。配置CZT和PA模塊時鐘頻率分別為400 MHz和250 MHz，此時系統(tǒng)延時lT=lFIR+lCZT+lPA=31.098 μs。

表5 非平穩(wěn)分析滿足連續(xù)實時處理時鐘頻率

其Modelsim時序仿真所得波形圖如圖13所示。

時頻分析結(jié)果如圖14所示。

通過圖13和圖14可得各信號頻率、功率和出現(xiàn)時刻(二者所含信息相同)。具體檢測參數(shù)如表6所示。

本設(shè)計檢測時間由計算時延和重疊幀時延組成，其中，計算時延可通過減少重疊點數(shù)和提高系統(tǒng)頻率與采樣率比值降低。通過表6可知，本設(shè)計對非平穩(wěn)信號檢測具有很好的實時性能和準(zhǔn)確度。

圖13 非平穩(wěn)信號時序仿真波形圖Fig.13 Timing simulation waveform of non-stationary signal

圖14 Modelsim檢測結(jié)果時頻圖Fig.14 Modelsim test results time-frequency diagram

5 與其他實時譜分析儀比較

表7對比了不同實時譜分析儀在核心算法、頻點數(shù)、實時帶寬(最大分析帶寬)、最小頻率分辨率(最大實時帶寬下)、近似成本和是否可分析非平穩(wěn)信號方面的異同。通過分析可知,本設(shè)計系統(tǒng)在實時帶寬和最小頻率分辨率方面具有與安捷倫和泰克實時譜分析儀相近的性能，且都可以對非平穩(wěn)信號進行分析，因此在分析信號基本譜信息應(yīng)用場合具有更高的性價比。同另外2種基于FPGA的實時譜分析儀相比，本設(shè)計的譜算法具有更優(yōu)越的性能，且成本相近。

表6 非平穩(wěn)信號測量結(jié)果對比

6 結(jié) 論

本文在結(jié)合短時傅里葉變換對非平穩(wěn)信號的分析方法和CZT算法優(yōu)點的基礎(chǔ)上，針對平穩(wěn)信號及非平穩(wěn)信號功率譜的實時測量，提出了一種基于多相濾波原理、CZT變換和譜分析相結(jié)合的便于硬件實現(xiàn)的實時譜分析方法，并利用FPGA平臺對該方法硬件實現(xiàn)。該方法利用多相濾波重構(gòu)用戶關(guān)注譜，利用CZT頻譜細(xì)化算法解決了窗口長度對頻譜分辨率的限定，譜分析算法提供了精確的譜分析結(jié)果。仿真結(jié)果及誤差分析表明,該系統(tǒng)可準(zhǔn)確檢測平穩(wěn)信號并對其進行譜分析，檢測頻率誤差小于0.6%、功率誤差小于4.5%和各實例系統(tǒng)最小延時均小于37 μs。而對于長度為32.768 μs非平穩(wěn)信號該系統(tǒng)最大頻率檢測誤差為94 kHz(頻譜分辨率為61 kHz)，最大時間檢測誤差為0.836 μs,最大時間延時為69 μs (重疊樣點為460點)，最小時間分辨率和頻率分辨率可分別為1.66 μs和40.69 kHz。仿真分析表明本設(shè)計對于平穩(wěn)信號及非平穩(wěn)信號譜分析都具有很好的檢測準(zhǔn)確度、時頻精度和實時性。

表7 不同實時譜分析儀對比