周 圍，梁 琦

(1.重慶郵電大學(xué) 移動(dòng)通信技術(shù)重慶市重點(diǎn)實(shí)驗(yàn)室，重慶 400065; 2.重慶郵電大學(xué) 光電工程學(xué)院，重慶 400065)

0 引 言

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

1 輸入輸出描述

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

(1)

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

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

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

2.1 CZT算法原理

設(shè)通過(guò)ADC采樣信號(hào)為x(n)(n=0,…,N-1)。根據(jù)文獻(xiàn)[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平面上起始點(diǎn)幅值，θ0為起始點(diǎn)相位(對(duì)應(yīng)起始點(diǎn)頻率)，W0為幅值變化梯度倒數(shù)，φ0為相位增量(對(duì)應(yīng)頻率分辨率)。A0=W0=1，使CZT變換路徑為單位圓上的一段圓弧，令M=N(M為譜線數(shù))簡(jiǎn)化CZT計(jì)算復(fù)雜度，通過(guò)θ0調(diào)整起始頻率，通過(guò)φ0確定選頻范圍內(nèi)細(xì)化倍數(shù)。

(2)式的計(jì)算可通過(guò)圖2所示步驟來(lái)實(shí)現(xiàn)。

圖2 CZT計(jì)算框圖Fig.2 Block diagram of CZT calculation

圖2中，f(n)與h(n)的離散卷積可通過(guò)其循環(huán)卷積實(shí)現(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。否則將無(wú)法滿足用戶輸入IDR的要求。

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

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

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

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

(4)

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

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

(5)

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

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

2.3 譜分析算法原理

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

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

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

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

(6)

(7)

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

PFA|M?1≈e-γ

(8)

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

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

(9)

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

γ≈-ln(PFA,g/MCZT)

(10)

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

2.4 多相濾波原理

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

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

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

(11)

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

(12)

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

(13)

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

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

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

1)根據(jù)用戶輸入IDR=30 dB，Ns=4，由(4)式可知，SLL≥40。在此，窗函數(shù)選擇Hamming窗(根據(jù)文獻(xiàn)[5]可知，Hamming窗SSL=-43 dB，SL=1.78 dB和主瓣點(diǎn)數(shù)=4)以滿足IDR要求，且nb=5，S=33。根據(jù)起止頻率與系統(tǒng)采樣率得D=16，根據(jù)歸一化分辨率公式、頻率分辨率及(5)式得NCZT=M=512。由起始頻率與初始相位對(duì)應(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個(gè)部分組成。硬件實(shí)現(xiàn)的關(guān)鍵問(wèn)題是模塊間的跨時(shí)域傳輸和模塊內(nèi)的流水線處理，因此，本文通過(guò)系統(tǒng)延時(shí)分析確定各時(shí)鐘域時(shí)鐘頻率。下面將對(duì)譜分析儀各主要模塊的硬件實(shí)現(xiàn)進(jìn)行詳細(xì)敘述。

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

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

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

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

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

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

圖8 CZT算法硬件實(shí)現(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 譜分析實(shí)現(xiàn)

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

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

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

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

(14)

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

(15)

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

3.4 系統(tǒng)時(shí)延

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

(16)

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

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

(17)

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

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

(18)

通過(guò)上述各模塊延時(shí)分析可知lT≤3·lFIR。

4 系統(tǒng)驗(yàn)證

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

4.1 平穩(wěn)信號(hào)性能驗(yàn)證

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

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

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

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

(19)

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

由圖11所得檢測(cè)信號(hào)頻率值和功率值。輸入信號(hào)與檢測(cè)信號(hào)頻率和歸一化功率值如表2所示。

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

表2 實(shí)例1信號(hào)測(cè)量結(jié)果對(duì)比

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

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

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

(20)

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

由圖12所得檢測(cè)信號(hào)頻率值和功率值。輸入信號(hào)與測(cè)量信號(hào)頻率和歸一化功率值如表4所示。

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

表4 實(shí)例2信號(hào)測(cè)量結(jié)果對(duì)比

4.2 非平穩(wěn)信號(hào)性能驗(yàn)證

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

(21)

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

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

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

其Modelsim時(shí)序仿真所得波形圖如圖13所示。

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

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

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

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

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

5 與其他實(shí)時(shí)譜分析儀比較

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

表6 非平穩(wěn)信號(hào)測(cè)量結(jié)果對(duì)比

6 結(jié) 論

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

表7 不同實(shí)時(shí)譜分析儀對(duì)比