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