先聲明一下這個(gè)項(xiàng)目本身不算新但確實(shí)是電力系統(tǒng)同步相量測(cè)量領(lǐng)域繞不開的經(jīng)典組合。FFT、窗函數(shù)法、希爾伯特-黃變換、小波變換四種方法放在一起做電力系統(tǒng)同步相量計(jì)算聽起來像是拼盤實(shí)際做下來會(huì)發(fā)現(xiàn)它們各有各的脾氣有的適合穩(wěn)態(tài)精確測(cè)量有的擅長動(dòng)態(tài)跟蹤有的純粹是離線分析利器。本文就圍繞這個(gè)研究題目把這四種方法在Matlab里的實(shí)現(xiàn)思路、代碼框架、實(shí)測(cè)對(duì)比和那些文檔里不會(huì)寫的坑完完整整梳理一遍。如果你是為了課程設(shè)計(jì)、畢業(yè)設(shè)計(jì)或者找算法落地方向這篇可以直接當(dāng)參考如果你是剛接觸同步相量這一塊我也盡量把原理用大白話講透保證你能拿著代碼跑起來、看得懂結(jié)果。1. 同步相量計(jì)算要解決什么問題1.1 從PMU說起同步相量測(cè)量到底在測(cè)什么電力系統(tǒng)同步相量Synchrophasor本質(zhì)上是借助全球衛(wèi)星定位系統(tǒng)北斗或GPS這里僅指通用授時(shí)服務(wù)提供的統(tǒng)一時(shí)標(biāo)對(duì)不同變電站的電壓或電流波形進(jìn)行高精度同步采樣再通過算法估算出基波分量的幅值、相位和頻率從而刻畫整個(gè)電網(wǎng)的動(dòng)態(tài)運(yùn)行狀態(tài)。這個(gè)測(cè)量裝置就是同步相量測(cè)量單元PMU它輸出的數(shù)據(jù)廣泛應(yīng)用于動(dòng)態(tài)監(jiān)測(cè)、故障分析和廣域控制。你可能覺得測(cè)量電壓電流的幅值和相位不就是一個(gè)傅里葉變換的事嗎理論上是這樣但工程上有個(gè)核心矛盾PMU要求算法既能適應(yīng)穩(wěn)態(tài)電網(wǎng)又能在系統(tǒng)頻率偏離工頻、幅值波動(dòng)、諧波污染、低頻振蕩等動(dòng)態(tài)條件下保持測(cè)量精度。國家標(biāo)準(zhǔn)和IEEE標(biāo)準(zhǔn)比如C37.118系列對(duì)相量測(cè)量誤差都有硬性指標(biāo)總向量誤差TVE通常要求穩(wěn)態(tài)時(shí)小于1%動(dòng)態(tài)條件下也有嚴(yán)格的容限。這就決定了你不能簡單地套一個(gè)DFT就完事而是要根據(jù)測(cè)量場(chǎng)景選擇合適的算法組合。我在做這個(gè)項(xiàng)目的初期第一個(gè)跳進(jìn)腦子的問題就是無非就是算個(gè)基波相量需要把FFT、窗函數(shù)、HHT、小波變換全用上嗎是不是過度設(shè)計(jì)了等我把穩(wěn)態(tài)和動(dòng)態(tài)的測(cè)試場(chǎng)景擺出來后才發(fā)現(xiàn)每一種方法都有它不可替代的位置。FFT是基準(zhǔn)窗函數(shù)法是工程上最實(shí)用的改進(jìn)HHT和小波變換則是處理非平穩(wěn)信號(hào)、揭示相量動(dòng)態(tài)過程的重要手段。所以這個(gè)題目看著像算法合集實(shí)際是一個(gè)由淺入深、從工程到學(xué)術(shù)的完整鏈路。1.2 算法研究的痛點(diǎn)傳統(tǒng)DFT在動(dòng)態(tài)條件下為什么不夠用先看最經(jīng)典的離散傅里葉變換DFT。教科書告訴我們DFT的前提是信號(hào)是周期性的并且采樣窗長度恰好等于信號(hào)周期的整數(shù)倍?？呻娋W(wǎng)的頻率是會(huì)波動(dòng)的比如正常是50Hz負(fù)荷突變、故障、新能源接入時(shí)可能變成49.8Hz或者50.3Hz。此時(shí)你用固定采樣率、固定窗長的DFT去做計(jì)算窗內(nèi)不再是整數(shù)個(gè)基波周期頻譜發(fā)生了泄漏——本來應(yīng)該集中在50Hz頻點(diǎn)上的能量分散到了相鄰頻點(diǎn)測(cè)出來的幅值偏小、相位偏歪TF值電網(wǎng)頻率也失真。想象一下你拿著一把固定長度的尺子去量一根熱脹冷縮的鋼筋鋼筋變長了但尺子沒變測(cè)量結(jié)果必然有誤差。DFT就是這個(gè)“固定尺子”電網(wǎng)頻率一動(dòng)它就露餡。所以傳統(tǒng)DFT在穩(wěn)態(tài)條件下確實(shí)夠用誤差可以做到非常小但一旦進(jìn)入動(dòng)態(tài)場(chǎng)景頻率偏移、諧波疊加、幅值調(diào)制就必須引入窗函數(shù)來抑制頻譜泄漏或者改用其他時(shí)頻分析手段。這個(gè)項(xiàng)目里實(shí)際用到的四類方法本質(zhì)上就是圍繞“如何在不同工況下把基波相量算準(zhǔn)”這個(gè)問題展開的。2. 四種算法的原理與選型邏輯2.1 FFT與窗函數(shù)法工程精度最優(yōu)解FFT快速傅里葉變換只是DFT的高效算法實(shí)現(xiàn)核心思路沒變把時(shí)域信號(hào)分解成一系列正弦分量的疊加。在同步相量計(jì)算中我們關(guān)心的是基波50Hz附近對(duì)應(yīng)的那個(gè)復(fù)數(shù)譜線幅值對(duì)應(yīng)相量幅值角度對(duì)應(yīng)相量相位。Matlab里的fft函數(shù)直接可用但這只是起點(diǎn)。窗函數(shù)法解決的是FFT的頻譜泄漏問題。常見的窗有漢寧窗Hanning、海明窗Hamming、布萊克曼窗Blackman等選取原則是在主瓣寬度和旁瓣衰減之間做權(quán)衡。漢寧窗因?yàn)榕园晁p快、實(shí)現(xiàn)簡單在電力系統(tǒng)諧波分析里用得最多。加窗的本質(zhì)是把信號(hào)兩端不連續(xù)的地方進(jìn)行平滑抑制截?cái)嘁鸬念l譜虛假擴(kuò)展。我實(shí)測(cè)下來加窗后幅值恢復(fù)要乘一個(gè)修正系數(shù)這個(gè)很多人會(huì)漏。Hanning窗在N點(diǎn)下的相干增益不等于N而是約等于N/2。如果你用sum(win)去做歸一化幅值計(jì)算結(jié)果會(huì)準(zhǔn)確如果直接用N/2近似點(diǎn)數(shù)非整周期時(shí)會(huì)有偏差。這個(gè)細(xì)節(jié)直接影響TVE指標(biāo)后面代碼環(huán)節(jié)再細(xì)說。不過加窗也有代價(jià)主瓣變寬頻率分辨率變差。也就是說你能測(cè)出基波附近有能量但區(qū)分不開49Hz和51Hz兩個(gè)鄰近分量。在同步相量計(jì)算這個(gè)場(chǎng)景里基波和鄰近間諧波的分辨通常不是首要矛盾首要矛盾是基波自身的幅值精度和相位精度所以加窗是個(gè)劃算的買賣。2.2 希爾伯特-黃變換自適應(yīng)的非平穩(wěn)信號(hào)分析利器希爾伯特-黃變換Hilbert-Huang TransformHHT是黃鍔提出的方法核心分為兩步第一步用經(jīng)驗(yàn)?zāi)B(tài)分解EMD把信號(hào)自適應(yīng)地分解成本征模態(tài)函數(shù)IMF第二步對(duì)每個(gè)IMF做Hilbert變換得到瞬時(shí)幅值和瞬時(shí)頻率。為什么要先EMD再Hilbert因?yàn)镠ilbert變換要求信號(hào)是窄帶的單分量信號(hào)否則求出來的瞬時(shí)頻率沒有物理意義。EMD就是把一個(gè)復(fù)雜信號(hào)拆成若干個(gè)“相對(duì)單分量”的IMF這些IMF按頻率從高到低排列每個(gè)IMF滿足兩個(gè)條件極值點(diǎn)數(shù)與過零點(diǎn)數(shù)相等或最多差1上下包絡(luò)關(guān)于時(shí)間軸對(duì)稱。你可以把EMD理解為一種自適應(yīng)的濾波器組它不需要預(yù)先知道信號(hào)的頻率范圍也不需要選擇小波基所有分解都由信號(hào)自身驅(qū)動(dòng)。在同步相量計(jì)算場(chǎng)景里HHT的價(jià)值特別明顯。傳統(tǒng)FFT假定信號(hào)在一個(gè)時(shí)間窗內(nèi)是平穩(wěn)的可實(shí)際電網(wǎng)中頻率是連續(xù)變化的幅值也可能因?yàn)榈皖l振蕩而波動(dòng)。HHT通過EMD把基波分量抽取出來再用Hilbert變換獲得每個(gè)采樣時(shí)刻的瞬時(shí)頻率和瞬時(shí)相位天然適合描述這類時(shí)變過程。它的缺點(diǎn)也很突出EMD的分解過程沒有嚴(yán)格的數(shù)學(xué)公式支撐端點(diǎn)效應(yīng)、模態(tài)混疊等問題處理起來很考驗(yàn)經(jīng)驗(yàn)算法計(jì)算量比較大實(shí)時(shí)性不如FFT。所以我在項(xiàng)目里把HHT定位為離線動(dòng)態(tài)分析工具用來觀察相量的瞬時(shí)軌跡而不是作為在線測(cè)量主算法。2.3 小波變換帶有頻率伸縮能力的顯微鏡小波變換Wavelet Transform和FFT、窗函數(shù)法的區(qū)別在于它通過一個(gè)可以伸縮和平移的小波基函數(shù)去匹配信號(hào)在不同時(shí)間位置的局部特征。低頻部分用寬窗口看頻率分辨率高、時(shí)間分辨率低高頻部分用窄窗口看時(shí)間分辨率高、頻率分辨率低。這種多分辨率特性非常適合分析非平穩(wěn)信號(hào)。在電力系統(tǒng)同步相量領(lǐng)域連續(xù)小波變換CWT用得比較多因?yàn)樗茉跁r(shí)頻平面上畫出幅值和相位隨時(shí)間的變化。比如你選擇復(fù)Morlet小波它本質(zhì)上是復(fù)指數(shù)乘高斯窗幅頻特性近似帶通濾波器尺度參數(shù)a對(duì)應(yīng)中心頻率f fc/(a*fs)這里fc是小波的中心頻率fs是采樣頻率。通過遍歷不同尺度就能在不同頻率點(diǎn)上看信號(hào)成分的強(qiáng)度變化找到基波頻率的位置后直接提取該位置的小波系數(shù)得到基波相量的瞬時(shí)估計(jì)。Matlab里CWT的實(shí)現(xiàn)經(jīng)歷了幾個(gè)版本的變化。早期用cwt(x, scales, morl)后來R2016b引入了新接口cwt(x, fs)默認(rèn)使用Morse小波R2021a之后默認(rèn)的解析小波變成了Morletcwt(x, amor, fs)或cwt(x, morse, fs)。這些接口變化是很多初學(xué)者運(yùn)行舊代碼報(bào)錯(cuò)的重災(zāi)區(qū)我在后面會(huì)專門列出來。2.4 方法對(duì)比什么時(shí)候用什么為什么把四種方法放在同一張表里對(duì)比方便你根據(jù)項(xiàng)目需求做選型方法原理核心優(yōu)勢(shì)劣勢(shì)適用場(chǎng)景FFT頻域分解計(jì)算快、實(shí)現(xiàn)簡單、穩(wěn)態(tài)精度極高頻譜泄漏、柵欄效應(yīng)穩(wěn)態(tài)相量測(cè)量基準(zhǔn)窗函數(shù)法加窗抑制泄漏工程實(shí)用、在線可用主瓣變寬、動(dòng)態(tài)跟蹤能力弱頻率偏移小、諧波場(chǎng)景希爾伯特-黃變換EMD瞬時(shí)頻率自適應(yīng)、適合非平穩(wěn)信號(hào)實(shí)時(shí)性差、端點(diǎn)效應(yīng)、經(jīng)驗(yàn)性強(qiáng)離線動(dòng)態(tài)分析、振蕩分析小波變換時(shí)頻域多分辨率時(shí)間頻率同時(shí)定位、抗噪能力強(qiáng)尺度選擇影響大、計(jì)算量大動(dòng)態(tài)相量跟蹤、時(shí)頻圖譜從工程落地角度看PMU現(xiàn)有標(biāo)準(zhǔn)算法大多還是DFT/加窗DFT這一模式因?yàn)樗诰€性能好、確定性高、便于實(shí)時(shí)實(shí)現(xiàn)。HHT和小波變換更像是“研究工具”用來分析動(dòng)態(tài)行為、驗(yàn)證理論或者在特殊故障信號(hào)中做后處理。我在項(xiàng)目里把四種方法都跑了基準(zhǔn)測(cè)試這樣既能說明各方法的適用邊界也能給后續(xù)改進(jìn)提供對(duì)照。3. 基于Matlab的仿真系統(tǒng)搭建3.1 仿真信號(hào)模型與測(cè)試場(chǎng)景設(shè)計(jì)做算法研究的第一步不是寫代碼是先把測(cè)試信號(hào)定義清楚。我按IEEE標(biāo)準(zhǔn)和電力系統(tǒng)動(dòng)態(tài)特性的常見情況設(shè)計(jì)了以下幾種仿真場(chǎng)景穩(wěn)態(tài)場(chǎng)景信號(hào)為額定頻率50Hz、幅值100V或標(biāo)幺值1.0、初相角30度的標(biāo)準(zhǔn)正弦波。理論上此時(shí)任何算法都能測(cè)準(zhǔn)主要用來驗(yàn)證代碼正確性和基線誤差。頻率偏移場(chǎng)景頻率分別設(shè)為49.5Hz、50.5Hz、51Hz模擬電網(wǎng)偏離工頻。此時(shí)固定窗長的DFT會(huì)產(chǎn)生頻譜泄漏是最考驗(yàn)窗函數(shù)法和插值算法的場(chǎng)景。諧波污染場(chǎng)景在基波上疊加3、5、7次諧波各次諧波幅值按基波的10%、5%、3%取值模擬實(shí)際電網(wǎng)中非線性負(fù)載導(dǎo)致的波形畸變。動(dòng)態(tài)調(diào)制場(chǎng)景幅值按1Hz正弦調(diào)制調(diào)制深度0.1模擬低頻振蕩。這種場(chǎng)景下幅值和相位是時(shí)變的適合用HHT和小波變換展示瞬時(shí)相量軌跡。采樣參數(shù)我統(tǒng)一設(shè)為采樣率fs 1600Hz即每個(gè)工頻周期采32個(gè)點(diǎn)符合同步相量測(cè)量的常見配置分析數(shù)據(jù)窗長度取N 128約一個(gè)多工頻周期FFT點(diǎn)數(shù)也用128窗函數(shù)選Hanning窗。為什么選1600Hz而不是更高的采樣率一方面是奈奎斯特采樣定理的一般要求分析到幾十次諧波綽綽有余另一方面是數(shù)據(jù)量不至于太大方便批量仿真對(duì)比算法。如果你要分析更高次諧波或間諧波采樣率需要相應(yīng)提高但基礎(chǔ)框架不變。3.2 Matlab代碼實(shí)現(xiàn)基線與加窗DFT信號(hào)生成的Matlab代碼如下這段代碼放在整個(gè)仿真程序的開頭作為信號(hào)源% 仿真信號(hào)生成函數(shù) % fs: 采樣率, f0: 基波頻率, A: 基波幅值, phi0: 初相角 % harm_amps: 諧波幅值向量, harm_freqs: 諧波頻率向量 function x gen_signal(t, f0, A, phi0, harm_amps, harm_freqs) x A * cos(2*pi*f0*t phi0); for k 1:length(harm_amps) x x harm_amps(k) * cos(2*pi*harm_freqs(k)*t phi0*randn); end end % 主程序示例 fs 1600; % 采樣率 duration 0.2; % 仿真時(shí)長, 0.2秒 t (0:round(duration*fs)-1)/fs; % 穩(wěn)態(tài): 50Hz, 幅值100, 初相30度 x_steady gen_signal(t, 50, 100, 30*pi/180, [], []); % 頻率偏移: 49.5Hz x_off gen_signal(t, 49.5, 100, 30*pi/180, [], []); % 諧波: 3次10%, 5次5%, 7次3% x_harm gen_signal(t, 50, 100, 30*pi/180, [10 5 3], [150 250 350]);接下來是基波相量提取函數(shù)。我用一個(gè)子函數(shù)封裝DFT/加窗DFT的相量計(jì)算方便后面批量調(diào)用function [amp, pha] phasor_dft(x, fs, f0, win_type) N length(x); % 構(gòu)造基波對(duì)應(yīng)的FFT頻點(diǎn)索引 k0 round(f0 / fs * N) 1; % Matlab索引從1開始 % 加窗 switch win_type case rect w ones(N,1); case hann w hann(N, periodic); case hamming w hamming(N, periodic); end % 加窗FFT X fft(x(:) .* w(:)); Xk X(k0); % 幅值恢復(fù)除以窗的相干增益 cg sum(w); % 相干增益, 矩形窗為N, hann窗約為N/2 amp 2 * abs(Xk) / cg; pha angle(Xk); % 注意對(duì)于全周期采樣角度直接就對(duì)應(yīng)相位非全周期時(shí)要修正 end這一段有幾個(gè)關(guān)鍵點(diǎn)值得解釋。第一hann(N, periodic)和hann(N, symmetric)的區(qū)別要搞清楚periodic窗的首尾值都是0更適合頻譜分析它保證窗長為N時(shí)N-1個(gè)點(diǎn)的周期延拓是連續(xù)的symmetric則用于FIR濾波器設(shè)計(jì)。同步相量計(jì)算中建議用periodic。第二幅值恢復(fù)公式amp 2*abs(Xk)/cg是全周期采樣情況下嚴(yán)格執(zhí)行的公式。矩形窗時(shí)cgN加窗后幅值衰減了一半所以要乘2恢復(fù)。Hanning窗時(shí)cg≈N/2所以2*N/(N/2)4的修正系數(shù)很多初學(xué)者直接用2/N恢復(fù)結(jié)果幅值偏大或偏小這其實(shí)是對(duì)“2”的來歷理解不透導(dǎo)致的。第三頻點(diǎn)索引k0 round(f0/fs*N)1。當(dāng)f0不是fs/N的整數(shù)倍時(shí)比如fs1600N128頻率分辨率12.5Hz50Hz對(duì)應(yīng)的k549.5Hz對(duì)應(yīng)的k4.96取整到5FFT峰值譜線的位置是固定的但因?yàn)閷?shí)際頻率不在整數(shù)頻點(diǎn)上就存在柵欄效應(yīng)。這種情況下加窗也只能抑制泄漏并不能消除頻率偏差帶來的幅值和相位誤差。要想真正解決需要插值算法我后面專門寫一節(jié)。3.3 HHT算法實(shí)現(xiàn)EMD分解Hilbert瞬時(shí)參數(shù)Matlab在Signal Processing Toolbox里提供了內(nèi)置的emd函數(shù)這個(gè)函數(shù)從R2018a開始可以用。用法如下% 對(duì)含諧波和調(diào)制的信號(hào)做EMD [imf, residual, info] emd(x_harm); % 畫出IMF觀察哪些分量對(duì)應(yīng)基波 figure; for i 1:4 subplot(4,1,i); plot(t, imf(i,:)); endEMD的運(yùn)行結(jié)果通常把最高頻分量放在第一階IMF依次往低頻排。對(duì)典型的電網(wǎng)信號(hào)基波50Hz通常出現(xiàn)在第一或第二階IMF。你可能會(huì)遇到諧波和基波被分到同一階IMF的情況這通常是模態(tài)混疊需要做預(yù)處理。提取出基波分量后用Hilbert變換求瞬時(shí)幅值和瞬時(shí)頻率% 假設(shè) imf1 是基波分量 z hilbert(imf1); % 解析信號(hào) inst_amp abs(z); % 瞬時(shí)幅值 inst_phase unwrap(angle(z)); % 瞬時(shí)相位, 注意unwrap inst_freq diff(inst_phase)/(2*pi) * fs; % 瞬時(shí)頻率 % 瞬時(shí)頻率序列比原信號(hào)短1點(diǎn), 可做線性平滑或中值濾波 inst_freq medfilt1(inst_freq, 5);這里有一個(gè)經(jīng)驗(yàn)直接對(duì)原始信號(hào)做Hilbert變換求瞬時(shí)頻率結(jié)果會(huì)非常抖因?yàn)樵夹盘?hào)不是單分量信號(hào)。必須先EMD分解成IMF再做Hilbert變換這恰恰是HHT的完整含義。另一個(gè)坑是邊緣效應(yīng)。EMD分解的端點(diǎn)會(huì)出現(xiàn)包絡(luò)發(fā)散導(dǎo)致第一階IMF在首尾幾十個(gè)采樣點(diǎn)明顯失真。Hilbert變換本身在端點(diǎn)也有Gibbs效應(yīng)所以你看到的瞬時(shí)幅值曲線兩端會(huì)有抖動(dòng)。我通常把首尾各40個(gè)點(diǎn)直接舍棄只保留中間段的瞬時(shí)參數(shù)用于分析。如果你需要整段時(shí)間的準(zhǔn)確相量建議給兩端各延長半周期數(shù)據(jù)分解完再裁掉。實(shí)際做動(dòng)態(tài)調(diào)制場(chǎng)景時(shí)HHT的瞬時(shí)幅值曲線能清晰地看到1Hz的調(diào)制頻率瞬時(shí)相位經(jīng)過unwrap后是一條斜線斜率除以2π乘以fs就是實(shí)測(cè)頻率。這個(gè)能力是FFT加窗無法直接提供的FFT給的是一個(gè)窗內(nèi)平均的“靜態(tài)相量”HHT給的是逐點(diǎn)的“動(dòng)態(tài)相量”。3.4 CWT連續(xù)小波變換實(shí)現(xiàn)頻率脊提取Matlab中CWT的使用在2016年以后發(fā)生了明顯變化。新接口推薦用% 推薦做法直接指定解析小波和采樣率 [wt, freqs] cwt(x, amor, fs); % amor表示復(fù)Morlet小波 % 或者用默認(rèn)的Morse小波 [wt, freqs] cwt(x, fs);wt是復(fù)數(shù)矩陣行對(duì)應(yīng)頻率列對(duì)應(yīng)時(shí)間。找出每個(gè)時(shí)刻幅度最大對(duì)應(yīng)的那一行頻率就是基波頻率的估計(jì)軌跡這叫做小波脊wavelet ridge提取[~, idx] max(abs(wt), [], 1); f_est freqs(idx); % 頻率估計(jì)軌跡 % 提取基波頻率處的小波系數(shù) base_idx find(freqs 49.5 freqs 50.5); % 基波頻率附近 cw wt(base_idx, :); % 取每個(gè)時(shí)刻小波系數(shù) cw_phase angle(cw); cw_amp abs(cw);如果要在CWT中指定尺度而非頻率需要知道尺度與頻率的換算關(guān)系。對(duì)于Morlet小波中心頻率fc6默認(rèn)頻率f與尺度a的關(guān)系是a fc / (2*pi*f/fs)更精確地說在Matlab里用scal2frq(scales, wavelet, deltat)來換算。新接口cwt已經(jīng)直接返回freqs省去了手工換算的麻煩所以強(qiáng)烈建議用新接口。CWT的效果有一個(gè)直觀特征在時(shí)頻譜圖上基波的頻率脊會(huì)像一條隨著頻率波動(dòng)蜿蜒的線如果用FFT你只能看到一條固定的譜線看不出頻率隨時(shí)間的變化。這就是多分辨率時(shí)頻分析的視覺優(yōu)勢(shì)。實(shí)際提取動(dòng)態(tài)相量時(shí)CWT的瞬時(shí)幅值在首尾也有邊緣效應(yīng)小波支持域的長度限制處理思路和HHT類似。另外CWT用Morlet小波時(shí)小波基在頻域有一定的帶寬如果信號(hào)頻率偏移比較大而基波附近仍有能量提取頻帶就不能選太窄否則部分信號(hào)能量落在窗外。我在代碼里選擇了freqs 49.5 freqs 50.5這個(gè)范圍實(shí)際使用可根據(jù)你關(guān)心的頻偏范圍靈活調(diào)整。3.5 其他輔助實(shí)現(xiàn)代碼還包含了幾項(xiàng)輔助功能頻率粗測(cè)、真值計(jì)算、誤差指標(biāo)計(jì)算。頻率粗測(cè)先對(duì)信號(hào)去直流求過零間隔得到頻率的粗略估計(jì)再用這個(gè)估計(jì)值修正FFT頻點(diǎn)索引。function f_est rough_freq(x, fs) % 簡單過零檢測(cè) signs sign(x - mean(x)); zero_idx find(diff(signs) ~ 0); if length(zero_idx) 2 f_est 50; % 默認(rèn) return; end T_est (zero_idx(end) - zero_idx(1)) / fs; cycles length(zero_idx) - 1; f_est cycles / T_est; end真值計(jì)算仿真信號(hào)是我們自己生成的所以真值已知。對(duì)于調(diào)幅信號(hào)基波相位是phi(t) 2*pi*f0*t phi0頻率是f0本身實(shí)際調(diào)制會(huì)引入額外頻率成分相位是時(shí)間的一次函數(shù)這些信息直接由參數(shù)生成時(shí)保留即可。計(jì)算TVE的公式為TVE mean(abs((phaseor_est - phaseor_true) ./ phaseor_true)) * 100; % 百分比因?yàn)橄嗔渴菑?fù)數(shù)這里TVE定義為復(fù)數(shù)向量差與真值模的比值是IEEE標(biāo)準(zhǔn)中的總向量誤差。更多府上把幅值誤差和相位誤差分開看amp_err abs(amp_est - amp_true) / amp_true * 100; phase_err abs(wrapToPi(pha_est - pha_true)) * 180 / pi; % 度4. 實(shí)驗(yàn)結(jié)果對(duì)比與誤差分析4.1 穩(wěn)態(tài)與頻率偏移場(chǎng)景誰更穩(wěn)下面這個(gè)表是我在測(cè)試仿真中記錄的典型結(jié)果供參考場(chǎng)景方法幅值誤差(%)相位誤差(度)備注50Hz穩(wěn)態(tài)矩形窗FFT0.0002以下0.001以下全周期采樣理想情況50Hz穩(wěn)態(tài)Hanning窗FFT0.0005以下0.002以下窗功能量泄漏極小49.5Hz矩形窗FFT約1.5~2約5~10頻譜泄漏明顯49.5HzHanning窗FFT約0.2~0.5約1~2泄漏抑制有效49.5Hz加窗插值DFT小于0.1小于0.5顯著優(yōu)于單純加窗頻率偏移場(chǎng)景是最能體現(xiàn)窗函數(shù)價(jià)值的場(chǎng)景。矩形窗在50.5Hz時(shí)TVE可能超過2%而Hanning窗加插值修正后TVE基本能控制在1%以下。這里為什么插值這么關(guān)鍵因?yàn)镕FT只能給出離散頻點(diǎn)上的譜線當(dāng)信號(hào)頻率不在頻點(diǎn)上時(shí)峰值譜線兩側(cè)的譜線也能用來插值還原真實(shí)譜峰位置和幅度。最常見的是利用對(duì)數(shù)幅度譜在峰值點(diǎn)附近的二次插值或者用兩點(diǎn)比值法。以最簡單的三點(diǎn)插值Rife-Jane方法為例設(shè)k1是峰值譜線索引則頻率修正量Δ ((X_{k1}| - |X_{k-1}|) / (2|X_k| - |X_{k-1}| - |X_{k1}|))然后真實(shí)頻率就是(k1Δ)*fs/N。加窗后用這個(gè)修正量可以大幅提升頻率和幅值精度。篇幅所限我不把完整推導(dǎo)寫出來具體公式可以參考《FFT插值算法在PMU相量測(cè)量中的應(yīng)用》類文獻(xiàn)。4.2 諧波與動(dòng)態(tài)調(diào)制場(chǎng)景時(shí)頻方法顯身手諧波場(chǎng)景下的現(xiàn)象也很有意思。對(duì)50Hz基波疊加3、5、7次諧波矩形窗FFT在基波頻點(diǎn)也會(huì)受諧波旁瓣的污染TVE明顯變大加Hanning窗后因?yàn)榇暗母哳l旁瓣衰減很快諧波影響被有效抑制TVE能回落不少。動(dòng)態(tài)調(diào)制場(chǎng)景下FFT和窗函數(shù)法只能給出一個(gè)統(tǒng)計(jì)平均的相量值無法反映瞬時(shí)幅值變化。而HHT和CWT能直接輸出瞬時(shí)幅值軌跡兩者的波形都能跟著1Hz調(diào)制頻率同步起伏時(shí)間上的對(duì)位關(guān)系完全一致。我實(shí)際測(cè)試中HHT提取的基波瞬時(shí)幅值與理論調(diào)制曲線的一致性在信號(hào)中間段非常好兩端因?yàn)槎它c(diǎn)效應(yīng)誤差偏大CWT用小波脊提取的瞬時(shí)幅值在首尾各損失約一個(gè)支持長度但抗噪能力比HHT更好。下面這張表是動(dòng)態(tài)調(diào)制場(chǎng)景的對(duì)比方法能輸出瞬時(shí)相量能輸出瞬時(shí)頻率抗噪能力計(jì)算耗時(shí)(0.2s數(shù)據(jù))矩形窗FFT否否一般1msHanning窗DFT否否較好1msHHT是是一般約2sCWT是是較好約0.5sHHT的耗時(shí)主要花在EMD篩選過程的迭代上如果數(shù)據(jù)長度增加到幾秒鐘耗時(shí)增長明顯。CWT耗時(shí)可接受但想要更高的頻率分辨率需要增加尺度數(shù)耗時(shí)也會(huì)線性增長。這個(gè)項(xiàng)目本身是研究性質(zhì)的離線分析對(duì)實(shí)時(shí)性要求不高所以都可以接受。5. 常見問題與排查技巧實(shí)錄5.1 問題速查表把我在實(shí)際跑代碼過程中遇到的經(jīng)典問題整理成一張表問題現(xiàn)象可能原因解決方案幅值總是偏小或偏大窗函數(shù)的相干增益歸一化不對(duì)用sum(win)做歸一化不要直接用N/2近似瞬時(shí)頻率曲線在兩端劇烈抖動(dòng)Hilbert變換端點(diǎn)效應(yīng)、EMD端點(diǎn)發(fā)散舍棄首尾50個(gè)采樣點(diǎn)或用延拓預(yù)處理EMD把諧波和基波混在一起模態(tài)混疊先帶通濾波預(yù)處理或增加篩迭代次數(shù)cwt函數(shù)報(bào)錯(cuò)或提示未定義變量舊代碼用了cwtft/cwt(..., morl)舊語法改用新接口cwt(x, amor, fs)頻率偏移大時(shí)加窗FFT仍誤差較大柵欄效應(yīng)相位誤差結(jié)合插值算法如三點(diǎn)插值相位結(jié)果要加/減某個(gè)固定偏差沒有修正窗函數(shù)的相位延遲對(duì)加窗DFT相位做群延遲補(bǔ)償每條里面的問題我都踩過下面撿幾個(gè)展開說。5.2 FFT頻譜泄漏與柵欄效應(yīng)的連鎖反應(yīng)頻譜泄漏來自非整周期截?cái)鄸艡谛?yīng)來自離散頻點(diǎn)插值兩者經(jīng)常同時(shí)出現(xiàn)一起導(dǎo)致誤差。你單獨(dú)測(cè)49.5Hz信號(hào)時(shí)FFT的峰值頻點(diǎn)實(shí)際上是第5根譜線對(duì)應(yīng)50Hz而真實(shí)的49.5Hz落在了第4.96根譜線上這0.04的偏差就是柵欄效應(yīng)加窗只能減輕泄漏不能消除柵欄效應(yīng)。所以行內(nèi)經(jīng)驗(yàn)是把兩者聯(lián)合解決加窗抑制泄漏插值修正柵欄。常用的是雙譜線插值算法——在峰值譜和相鄰譜之間通過比值或?qū)?shù)比值求修正量計(jì)算量增加不多精度提升明顯。把這段代碼封裝好后在49.5Hz場(chǎng)景下TVE可以從2%級(jí)別降到0.5%級(jí)別效果立竿見影。5.3 Hilbert變換的相位解纏問題一定不能忘求瞬時(shí)相位時(shí)angle(z)得到的是主值相位在(-π, π]區(qū)間內(nèi)跳變。如果頻率穩(wěn)定瞬時(shí)相位是一條直線但主值相位會(huì)有鋸齒跳變所以必須用unwrap進(jìn)行解纏否則后續(xù)求導(dǎo)得到的瞬時(shí)頻率全是巨大尖峰。另外hilbert函數(shù)是對(duì)整段數(shù)據(jù)做FFT再變換的它的輸出長度和輸入長度相同邊界處有Gibbs振鈴所以提取低頻瞬時(shí)幅值的時(shí)候兩端也會(huì)出現(xiàn)波動(dòng)。處理辦法是拿信號(hào)中間段評(píng)估精度或者前后加反射延拓。5.4 小波變換尺度選擇與頻率分辨率的權(quán)衡CWT的尺度數(shù)決定了頻率分辨率。如果尺度數(shù)太少頻率脊可能忽上忽下不平滑尺度數(shù)太多計(jì)算耗時(shí)大且低頻段時(shí)間分辨率又不足。對(duì)同步相量場(chǎng)景我建議把頻率范圍設(shè)置在10Hz~200Hz之間尺度數(shù)以幾百到幾千為量級(jí)。新接口cwt默認(rèn)會(huì)基于信號(hào)長度自動(dòng)選頻率你只需要指定寬度即可。如果想要更高的頻率分辨率可以用cwtfreqbounds函數(shù)查看可用頻率范圍再手動(dòng)指定voices參數(shù)控制頻率細(xì)分。5.5 EMD的模態(tài)混疊與預(yù)處理手段模態(tài)混疊是EMD處理電網(wǎng)信號(hào)最常見的坑。如果諧波頻率和基波頻率不是整數(shù)倍比如存在間諧波49HzEMD容易出現(xiàn)一次篩選把兩個(gè)成分混在一起的情況導(dǎo)致IMF失去單分量特性。經(jīng)驗(yàn)做法是對(duì)信號(hào)先做帶通濾波比如30~70Hz濾掉諧波和噪聲再送入EMD。雖然這樣等于先入為主地選擇了頻帶但電力系統(tǒng)同步相量場(chǎng)景下我們明確關(guān)心基波所以帶通預(yù)處理是合理的。另一種方法是增大EMD的篩迭代次數(shù)讓分解更充分代價(jià)是耗時(shí)增加。具體迭代次數(shù)可在emd函數(shù)的名稱-值對(duì)參數(shù)中設(shè)置比如MaxNumIMF、SiftMaxIterations。6. 工程落地與擴(kuò)展思考做這類研究的時(shí)間久了我個(gè)人的體會(huì)是不管算法多花哨工程上衡量一個(gè)同步相量算法的標(biāo)準(zhǔn)永遠(yuǎn)是四個(gè)字——穩(wěn)、準(zhǔn)、快、簡單。在這個(gè)項(xiàng)目里加窗插值DFT是那個(gè)最“穩(wěn)”的方案它既有理論依據(jù)又能在線實(shí)時(shí)實(shí)現(xiàn)HHT和CWT雖然精度和動(dòng)態(tài)能力更強(qiáng)但它們更像是分析工具適合在離線場(chǎng)景或者非平穩(wěn)分析中發(fā)揮作用。如果后續(xù)要把這套方法推向?qū)嶋H系統(tǒng)我建議從兩個(gè)方向擴(kuò)展。第一個(gè)方向是自適應(yīng)參數(shù)優(yōu)化比如根據(jù)頻率粗測(cè)結(jié)果動(dòng)態(tài)調(diào)整窗長或者采樣點(diǎn)數(shù)盡量保持整周期采樣這樣一來加窗DFT的誤差能進(jìn)一步壓縮。第二個(gè)方向是結(jié)合卡爾曼濾波或自適應(yīng)陷波器用狀態(tài)空間模型把基波相量建模為復(fù)指數(shù)通過遞推估計(jì)逐點(diǎn)更新這樣既能得到比FFT更平滑的動(dòng)態(tài)軌跡又比HHT實(shí)時(shí)性好很多。實(shí)際上現(xiàn)在很多動(dòng)態(tài)相量測(cè)量研究都在往這個(gè)方向走我和同行交流時(shí)也常聽到類似的思路。最后分享一個(gè)實(shí)用小技巧我在項(xiàng)目測(cè)試時(shí)會(huì)把所有算法的估計(jì)結(jié)果保存成統(tǒng)一格式的數(shù)據(jù)結(jié)構(gòu)含估計(jì)幅值、相位、頻率、TVE寫一個(gè)通用的批量測(cè)試腳本遍歷所有場(chǎng)景和算法。這樣一來改一個(gè)參數(shù)、重跑一次仿真、對(duì)比一張表整個(gè)流程非常順手。這個(gè)習(xí)慣也推薦給每一個(gè)做算法研究的人。