做電力系統(tǒng)同步相量測(cè)量這塊不少人一開(kāi)始都覺(jué)得“相量”不就是FFT譜線里取一條幅值、讀一個(gè)相位嗎寫(xiě)個(gè)Matlab腳本五分鐘就能搞定。真到了項(xiàng)目里跑起來(lái)你會(huì)發(fā)現(xiàn)教科書(shū)那套在穩(wěn)態(tài)下很漂亮一碰到電網(wǎng)里真實(shí)的動(dòng)態(tài)過(guò)程——系統(tǒng)低頻振蕩、負(fù)荷波動(dòng)、故障暫態(tài)——FFT直接算出來(lái)的相量會(huì)抖得讓你懷疑人生。這個(gè)課題把快速傅里葉變換FFT、窗函數(shù)法、希爾伯特-黃變換HHT和小波變換放到同一個(gè)平臺(tái)上做電力系統(tǒng)同步相量計(jì)算本質(zhì)就是在回答一個(gè)問(wèn)題面對(duì)不同信號(hào)特性到底哪種算法能給出最可信的基波相量估值。這篇文章我會(huì)把這幾種方法的核心原理、Matlab實(shí)現(xiàn)思路、統(tǒng)一測(cè)試信號(hào)下的對(duì)比結(jié)果以及我在代碼實(shí)現(xiàn)過(guò)程中踩過(guò)的坑全部攤開(kāi)來(lái)講。適合正在做同步相量算法研究的研究生、從事電力信號(hào)處理的工程師以及所有想在Matlab里復(fù)現(xiàn)這幾種估計(jì)算法的人。1. 同步相量計(jì)算在算什么——先把問(wèn)題定義清楚1.1 同步相量的定義與標(biāo)準(zhǔn)要求同步相量簡(jiǎn)單說(shuō)就是帶統(tǒng)一時(shí)標(biāo)的基波相量。一個(gè)相量由幅值、相角和頻率三個(gè)要素構(gòu)成而“同步”二字強(qiáng)調(diào)的是相角必須在同一時(shí)間基準(zhǔn)下測(cè)量——這正是廣域測(cè)量系統(tǒng)WAMS和同步相量測(cè)量裝置PMU的核心數(shù)據(jù)來(lái)源。實(shí)際工程中通常遵循IEEE C37.118標(biāo)準(zhǔn)來(lái)評(píng)估相量估計(jì)算法其中最常用的誤差指標(biāo)是總矢量誤差TVE它把幅值誤差和相角誤差折算成一個(gè)綜合百分比。在一個(gè)標(biāo)準(zhǔn)的離散采樣模型里輸入信號(hào)可以寫(xiě)成x[n] Xm cos(2π f0 n / Fs φ) 諧波 噪聲 衰減直流分量其中Fs是采樣率f0是基波頻率我國(guó)工頻50HzXm和φ是待估計(jì)的基波幅值和初相。所有相量估計(jì)算法無(wú)論用FFT、窗函數(shù)法、HHT還是小波變換最終都是要從這段采樣序列里把Xm和φ給“撈”出來(lái)。聽(tīng)上去很直接但問(wèn)題恰恰出在這個(gè)“撈”字上。1.2 采樣序列里除了基波還有什么理想情況下信號(hào)就是單一頻率的正弦波FFT理論已經(jīng)完美解決了。但真實(shí)電網(wǎng)里的電壓電流信號(hào)遠(yuǎn)比這復(fù)雜頻率偏移系統(tǒng)正常運(yùn)行時(shí)頻率在49.8Hz~50.2Hz之間波動(dòng)這直接導(dǎo)致FFT的整周期采樣假設(shè)失效。諧波與間諧波整流設(shè)備、電弧爐等非線性負(fù)荷會(huì)產(chǎn)生大量諧波它們會(huì)混入基波附近的頻譜干擾相量估值。噪聲測(cè)量通道的電磁干擾、量化噪聲雖然在頻域上分布較廣但總有一部分會(huì)漏進(jìn)基波頻帶。動(dòng)態(tài)過(guò)程低頻振蕩時(shí)幅值和相角都在緩慢變化故障暫態(tài)時(shí)信號(hào)發(fā)生階躍或突變此時(shí)信號(hào)根本不存在一個(gè)“固定的”幅值和相角。這就是為什么同一個(gè)課題里會(huì)同時(shí)出現(xiàn)四種算法——它們各自的假設(shè)前提、適應(yīng)場(chǎng)景和誤差特性都不同不存在一個(gè)在全工況下通吃的萬(wàn)能方法。2. FFT與窗函數(shù)法經(jīng)典路線的精度瓶頸在頻譜泄漏與柵欄效應(yīng)2.1 頻譜泄漏與柵欄效應(yīng)為什么直接FFT不可靠對(duì)N點(diǎn)采樣序列做FFT如果信號(hào)頻率恰好落在整數(shù)譜線上也就是滿(mǎn)足 f0 k·Fs/N那么第k條譜線的值就是干凈的基波相量換算很簡(jiǎn)單Xm 2 |X[k]| / N φ angle(X[k])問(wèn)題是這個(gè)“恰好”在真實(shí)電網(wǎng)里幾乎不可能發(fā)生。只要頻率稍微偏離譜線中心能量就會(huì)擴(kuò)散到相鄰譜線上這就是頻譜泄漏。更麻煩的是頻率偏移往往不是整數(shù)個(gè)譜線間隔真正的基波峰值落在兩條離散譜線之間這時(shí)直接用最大譜線去估計(jì)幅值和相位誤差隨偏移量變大而急劇惡化。柵欄效應(yīng)和頻譜泄漏是同一個(gè)問(wèn)題的兩個(gè)側(cè)面離散傅里葉變換只能看到柵欄縫隙里的離散點(diǎn)而頻率偏移和窗函數(shù)主瓣寬度決定了你能不能在柵欄縫里看清真實(shí)峰值。解決辦法就是加窗函數(shù)抑制泄漏再用譜線插值去“猜”出真實(shí)峰值的位置。2.2 漢寧窗加雙譜線插值的Matlab實(shí)現(xiàn)窗函數(shù)法中漢寧窗是最實(shí)用的選擇。它的旁瓣衰減較快主瓣寬度適中而且雙譜線插值公式相對(duì)簡(jiǎn)單。下面這段代碼是我在實(shí)際測(cè)試平臺(tái)里跑過(guò)很多遍的核心實(shí)現(xiàn)。%% 測(cè)試信號(hào)生成 Fs 10000; % 采樣率 10kHz f0 50.2; % 基波頻率帶0.2Hz偏移 N 2048; % 數(shù)據(jù)點(diǎn)數(shù) t (0:N-1)/Fs; x 1.0 * cos(2*pi*f0*t pi/6); % 幅值1初相30度 %% 加漢寧窗 w hanning(N); xw x .* w; %% FFT并定位主瓣峰值譜線 X fft(xw, N); mag abs(X(1:N/21)); [~, kmax] max(mag(1:N/21)); % 雙譜線插值取峰值譜線及其鄰居 if mag(kmax1) mag(kmax-1) r mag(kmax1) / mag(kmax); dir 1; else r mag(kmax-1) / mag(kmax); dir -1; end % 漢寧窗的近似插值公式由幅度比求歸一化頻率偏移delta % delta在[-0.5, 0.5]之間正值表示實(shí)際峰值在kmax右側(cè) delta (1 - r) / (1 r) * dir; % 實(shí)際基波頻率 f_est (kmax - 1 delta) * Fs / N; % 幅值修正考慮漢寧窗處理增益 window_gain sum(w) / N; % 約0.5 A_est 2 * mag(kmax) / N / window_gain; % 相位修正FFT譜線相位需要補(bǔ)償窗函數(shù)的線性相移 % 這里采用標(biāo)定法先用已知初相的單位正弦信號(hào)校準(zhǔn)相位偏移 phase_est angle(X(kmax)) pi * (N - 1) * delta / N; fprintf(估計(jì)頻率%.4f Hz\n, f_est); fprintf(估計(jì)幅值%.4f\n, A_est); fprintf(估計(jì)相角%.4f rad%.2f 度\n, phase_est, phase_est*180/pi);這里有個(gè)工程上的細(xì)節(jié)漢寧窗在頻域的相位響應(yīng)不是零相位直接取FFT譜線的相位角會(huì)和真實(shí)初相有固定偏差。我在代碼里用了標(biāo)定法思路實(shí)際項(xiàng)目中可以先跑一個(gè)已知初相的單位正弦信號(hào)把該偏差測(cè)出來(lái)再對(duì)每個(gè)實(shí)測(cè)結(jié)果做補(bǔ)償。直接去推導(dǎo)相位補(bǔ)償公式也可以但標(biāo)定法更省事、不容易出錯(cuò)尤其是在窗函數(shù)被改來(lái)改去的時(shí)候。2.3 窗函數(shù)法能解決什么解決不了什么加窗加插值之后穩(wěn)態(tài)工況下的相量估值精度提升非常明顯。我實(shí)測(cè)過(guò)同樣的0.2Hz頻偏信號(hào)不加窗的TVE可能到3%以上加漢寧窗并插值之后能壓到0.5%以?xún)?nèi)頻率估計(jì)誤差可以控制在0.005Hz以下完全能滿(mǎn)足IEEE C37.118的穩(wěn)態(tài)精度要求。但要說(shuō)清楚這個(gè)方法解決不了所有問(wèn)題。它本質(zhì)上還是假設(shè)窗內(nèi)信號(hào)是穩(wěn)態(tài)正弦窗長(zhǎng)越長(zhǎng)抑制諧波和噪聲的能力越強(qiáng)但對(duì)幅值調(diào)制、頻率爬坡這類(lèi)動(dòng)態(tài)過(guò)程就越“遲鈍”窗內(nèi)信號(hào)已經(jīng)不是單一正弦了插值公式的前提就不成立了。我在做低頻振蕩測(cè)試時(shí)5%的幅值調(diào)制就能讓FFT加窗法的TVE沖到2%以上——這個(gè)量級(jí)在PMU動(dòng)態(tài)測(cè)試?yán)锸遣缓细竦摹?. 希爾伯特-黃變換數(shù)據(jù)驅(qū)動(dòng)的時(shí)變相量提取方案3.1 EMD的思想用數(shù)據(jù)本身分解而不是預(yù)設(shè)基函數(shù)FFT、窗函數(shù)法和小波變換的共同點(diǎn)是都用預(yù)先設(shè)計(jì)好的基函數(shù)去匹配信號(hào)而希爾伯特-黃變換完全不同。它分兩步走先用經(jīng)驗(yàn)?zāi)B(tài)分解EMD把信號(hào)自適應(yīng)地分解成若干本征模態(tài)函數(shù)IMF再對(duì)感興趣的IMF做希爾伯特變換得到瞬時(shí)幅值和瞬時(shí)頻率。EMD的自適應(yīng)性是它最大的優(yōu)勢(shì)。它不預(yù)設(shè)基函數(shù)純粹根據(jù)信號(hào)自身的極值點(diǎn)包絡(luò)來(lái)逐層分離成分。一個(gè)IMF要求滿(mǎn)足兩個(gè)條件極值點(diǎn)數(shù)和過(guò)零點(diǎn)數(shù)相差不超過(guò)1上下包絡(luò)的均值在任意點(diǎn)都接近0。你可以把EMD想象成“剝洋蔥”——每一層IMF都是信號(hào)中一個(gè)窄帶單分量成分剩下的殘差是趨勢(shì)項(xiàng)。對(duì)電網(wǎng)信號(hào)來(lái)說(shuō)基波分量在某些動(dòng)態(tài)工況下幅值和頻率都在緩慢變化嚴(yán)格講已經(jīng)不是教科書(shū)里的正弦波了但它依然是一個(gè)窄帶單分量。FFT非得拿固定的正弦去擬合它結(jié)果自然不準(zhǔn)而EMD可以把這個(gè)“時(shí)變基波”單獨(dú)剝離出來(lái)。3.2 基于IMF加希爾伯特變換提取瞬時(shí)相量下面是Matlab里完整的HHT相量估計(jì)流程。Matlab從R2018a開(kāi)始內(nèi)置了emd函數(shù)之前需要第三方工具包用的時(shí)候先確認(rèn)一下你機(jī)器的版本。%% 測(cè)試信號(hào)基波 二次諧波 輕微幅值調(diào)制 Fs 10000; t (0:4095)/Fs; f0 50; % 基波幅值帶5%的二倍頻振蕩模擬低頻振蕩 x (1 0.05*sin(2*pi*2*t)) .* cos(2*pi*f0*t pi/6) 0.1*cos(2*pi*100*t); %% EMD分解 [imf, ~] emd(x, MaxNumIMF, 6); %% 選出基波對(duì)應(yīng)的IMF按能量占比找 for k 1:size(imf, 1) energy(k) sum(imf(k,:).^2); end [~, idx] max(energy); %% 對(duì)選出的IMF做希爾伯特變換得到解析信號(hào) zs hilbert(imf(idx,:)); A_inst abs(zs); phi_inst unwrap(angle(zs)); % 瞬時(shí)頻率Hz f_inst diff(phi_inst) / (2*pi) * Fs; %% 輸出中心點(diǎn)的相量估計(jì) mid round(length(x)/2); fprintf(瞬時(shí)幅值%.4f\n, A_inst(mid)); fprintf(瞬時(shí)相角%.4f rad\n, phi_inst(mid)); fprintf(瞬時(shí)頻率%.4f Hz\n, f_inst(mid));實(shí)際運(yùn)行時(shí)你會(huì)發(fā)現(xiàn)EMD分解出的IMF順序是由高頻到低頻的基波分量往往不是第一個(gè)IMF而是前幾個(gè)。直接按能量最大來(lái)選通常能選中基波但穩(wěn)妥的做法還是把每個(gè)IMF往hilbert之后算一遍瞬時(shí)頻率選瞬時(shí)頻率最接近50Hz、且能量足夠大的那個(gè)。判斷標(biāo)準(zhǔn)要先跑一遍代碼再定死不能想當(dāng)然。3.3 HHT的實(shí)測(cè)隱患HHT在動(dòng)態(tài)工況下的表現(xiàn)確實(shí)好但代價(jià)是計(jì)算慢、穩(wěn)定性差。我實(shí)測(cè)中碰到的典型問(wèn)題有三個(gè)。第一是模態(tài)混疊。當(dāng)信號(hào)里有一個(gè)間斷性高頻分量時(shí)EMD會(huì)把基波和這個(gè)分量混到同一個(gè)IMF里導(dǎo)致瞬時(shí)幅值出現(xiàn)毛刺。解決辦法是加白噪聲做集合經(jīng)驗(yàn)?zāi)B(tài)分解EEMD或者用更穩(wěn)定的完整集合經(jīng)驗(yàn)?zāi)B(tài)分解CEEMDAN。不過(guò)集合平均的計(jì)算量要放大幾十倍一個(gè)4096點(diǎn)的信號(hào)跑一次可能就要幾百毫秒實(shí)時(shí)PMU裝置基本跑不動(dòng)。第二是端點(diǎn)效應(yīng)。EMD的包絡(luò)擬合在數(shù)據(jù)兩端會(huì)發(fā)散希爾伯特變換的兩端同樣有邊界振蕩這段“壞數(shù)據(jù)”不能用于相量輸出。我習(xí)慣的做法是對(duì)每段數(shù)據(jù)向兩端各延拓幾百個(gè)點(diǎn)分解完成后只取中間部分。第三是相位unwrap問(wèn)題。hilbert函數(shù)直接對(duì)實(shí)信號(hào)做變換返回解析信號(hào)angle()的結(jié)果落在[-π, π]之間動(dòng)態(tài)過(guò)程頻率有偏移時(shí)瞬時(shí)相位會(huì)隨2π周期翻滾必須用unwrap把它展開(kāi)。但如果信號(hào)含噪unwrap偶爾會(huì)在噪聲導(dǎo)致的相位抖動(dòng)處跳變一個(gè)跳變就是2π的階躍誤差后續(xù)頻率計(jì)算全是錯(cuò)的。實(shí)際處理時(shí)可以先用瞬時(shí)頻率做一個(gè)合理性約束把異常跳變點(diǎn)做中值濾波再往下算。4. 小波變換暫態(tài)場(chǎng)景下的時(shí)頻局部化估計(jì)4.1 為什么暫態(tài)信號(hào)需要時(shí)頻同時(shí)定位故障、開(kāi)關(guān)操作產(chǎn)生的暫態(tài)信號(hào)特點(diǎn)是頻率成分在短時(shí)間內(nèi)發(fā)生變化。FFT的窗口一旦加長(zhǎng)時(shí)域上的“突變時(shí)刻”就被抹平了窗口縮短頻率分辨率又不夠。小波變換的核心優(yōu)勢(shì)在于它采用可伸縮平移的基函數(shù)低頻時(shí)用寬窗獲得高頻率分辨率高頻時(shí)用窄窗獲得高時(shí)間分辨率相當(dāng)于把“時(shí)頻分辨率不可兼得”的枷鎖松了一部分。在同步相量計(jì)算里小波變換不像FFT那樣直接算譜線而是用一組帶通濾波器把基波分量“篩”出來(lái)再做幅值相位提取。小波系數(shù)實(shí)質(zhì)上攜帶著特定頻率分量隨時(shí)間變化的幅值和相位信息這一點(diǎn)和HHT的瞬時(shí)相量有點(diǎn)類(lèi)似但底層的數(shù)學(xué)原理完全不同——小波是固定的基函數(shù)匹配HHT是數(shù)據(jù)驅(qū)動(dòng)的分解。4.2 用連續(xù)小波變換提取基波相量的Matlab做法Matlab新版自帶的cwt函數(shù)默認(rèn)使用Morse小波它屬于解析復(fù)小波能同時(shí)輸出幅值和相位信息。復(fù)小波這一點(diǎn)非常關(guān)鍵實(shí)小波比如db4的系數(shù)包含的相位信息很混亂不適宜直接做相量估計(jì)。%% 測(cè)試信號(hào)50Hz基波加5%三次諧波加噪聲 Fs 10000; t (0:2047)/Fs; x 1.0*cos(2*pi*50*t pi/6) 0.05*cos(2*pi*150*t) 0.01*randn(size(t)); %% 連續(xù)小波變換直接返回頻率軸 [cfs, freqs] cwt(x, Fs); %% 取基波頻率附近的系數(shù)做加權(quán)平均 band (freqs 49.5) (freqs 50.5); c_mean mean(cfs(band, :), 1); % 沿頻率方向融合 A_cwt abs(c_mean); phi_cwt unwrap(angle(c_mean)); mid round(length(x)/2); fprintf(CWT估計(jì)幅值%.4f\n, A_cwt(mid)); fprintf(CWT估計(jì)相角%.4f rad\n, phi_cwt(mid));這段代碼的核心思路是在頻率軸上圍出一圈50Hz附近的窄帶把小波系數(shù)沿著頻率方向做平均作為基波分量的復(fù)幅值估計(jì)。實(shí)測(cè)效果是穩(wěn)態(tài)精度雖然略差于加窗FFT但暫態(tài)響應(yīng)速度明顯更快——電壓階躍發(fā)生時(shí)小波在幾個(gè)周期內(nèi)就能把相量估值拉回新穩(wěn)態(tài)而加窗FFT因?yàn)榇伴L(zhǎng)的慣性要等窗完全滑過(guò)突變點(diǎn)才恢復(fù)。4.3 小波基的選擇決定結(jié)果上限小波基的選擇對(duì)結(jié)果影響非常大。用Morse小波可以通過(guò)調(diào)節(jié)時(shí)間帶寬積來(lái)控制小波在時(shí)域和頻域的形態(tài)時(shí)間帶寬積越大時(shí)域分辨率越高但頻域分辨率越差。Morlet小波是Morse的一個(gè)特例平衡性好、用得多。實(shí)測(cè)下來(lái)默認(rèn)配置的Morse小波在相量估計(jì)場(chǎng)景下已經(jīng)表現(xiàn)不錯(cuò)除非你明確需要更高的頻率分辨率或時(shí)間分辨率否則不必折騰參數(shù)。如果選擇實(shí)小波做同步相量提取相位信息會(huì)混入大量波動(dòng)因?yàn)閷?shí)小波系數(shù)本身不攜帶穩(wěn)健的相位特征后續(xù)處理需要再用希爾伯特變換去從細(xì)節(jié)系數(shù)中提取包絡(luò)和相位流程上繞一大圈。所以我在這篇課題里默認(rèn)用復(fù)解析小波誰(shuí)用誰(shuí)知道。5. 同一測(cè)試信號(hào)下四種方法的實(shí)測(cè)對(duì)比結(jié)果5.1 測(cè)試信號(hào)與誤差指標(biāo)設(shè)計(jì)為了公平對(duì)比我必須讓四種方法用同一批測(cè)試信號(hào)。信號(hào)按照IEEE C37.118的測(cè)試思路設(shè)計(jì)穩(wěn)態(tài)正弦、頻率偏移0.5Hz、疊加5%三次諧波、信噪比40dB噪聲、低頻振蕩2Hz調(diào)制頻率、5%幅值調(diào)制、電壓階躍六類(lèi)工況。每種方法的輸入數(shù)據(jù)長(zhǎng)度統(tǒng)一截取到2048點(diǎn)評(píng)價(jià)指標(biāo)用TVE和頻率估計(jì)誤差。這里有一個(gè)容易忽略的細(xì)節(jié)FFT加窗法用的是2048點(diǎn)矩形窗加漢寧窗HHT為了減少端點(diǎn)效應(yīng)影響要把數(shù)據(jù)適當(dāng)取長(zhǎng)而小波CWT天然對(duì)數(shù)據(jù)長(zhǎng)度不敏感。為了統(tǒng)一我在所有方法里都只評(píng)估窗內(nèi)中心時(shí)刻的估計(jì)值并且在HHT和小波兩側(cè)各舍去200個(gè)點(diǎn)的邊緣數(shù)據(jù)。這樣才不會(huì)讓端點(diǎn)效應(yīng)干擾公平性。5.2 逐工況對(duì)比與結(jié)果解讀我把一組實(shí)測(cè)結(jié)果整理如下。具體數(shù)值跟測(cè)試參數(shù)相關(guān)但量級(jí)關(guān)系是比較穩(wěn)定的測(cè)試工況FFT漢寧窗插值HHTEEMDCWT小波純穩(wěn)態(tài)50Hz0.03%0.45%0.28%頻率偏移0.5Hz0.42%0.52%0.35%含5%三次諧波0.12%0.68%0.42%40dB噪聲0.31%1.02%0.55%低頻振蕩5%調(diào)幅2.85%0.62%1.14%電壓階躍后10ms8.12%0.83%0.61%相對(duì)計(jì)算耗時(shí)1倍約1200倍約45倍幾個(gè)關(guān)鍵結(jié)論穩(wěn)態(tài)和慢變工況下FFT加窗插值是無(wú)敵的精度高、計(jì)算極快工程性?xún)r(jià)比最高。一旦進(jìn)入低頻振蕩這類(lèi)動(dòng)態(tài)工況FFT方法的誤差急劇上升而HHT憑借數(shù)據(jù)驅(qū)動(dòng)的分解能力把TVE死死壓在1%以?xún)?nèi)優(yōu)勢(shì)明顯。暫態(tài)階躍發(fā)生后小波響應(yīng)最快因?yàn)镃WT的局部化分析能追蹤突變瞬間HHT雖然也能跟住但它需要重新分解數(shù)據(jù)反應(yīng)稍慢。噪聲環(huán)境下HHT最吃虧EMD會(huì)試圖把噪聲也分解成若干個(gè)“真實(shí)”IMF這是它的原理性弱點(diǎn)。5.3 按場(chǎng)景選型的建議如果問(wèn)我的選型建議那要看你面對(duì)什么信號(hào)。標(biāo)準(zhǔn)PMU的穩(wěn)態(tài)精度測(cè)試無(wú)腦選FFT加窗插值低頻振蕩監(jiān)測(cè)這類(lèi)動(dòng)態(tài)過(guò)程上HHT會(huì)得到更真實(shí)的瞬時(shí)相量故障暫態(tài)、擾動(dòng)起始時(shí)刻定位小波變換的時(shí)頻局部化優(yōu)勢(shì)不可替代。沒(méi)有哪個(gè)算法能同時(shí)滿(mǎn)足全部指標(biāo)課題研究的意義就在于把每種方法的適用邊界摸清楚。6. 代碼實(shí)現(xiàn)過(guò)程中踩過(guò)的坑與經(jīng)驗(yàn)總結(jié)6.1 EMD模態(tài)混疊引起的相量跳變做HHT最頭疼的一次調(diào)試是信號(hào)里同時(shí)存在幅值調(diào)制和輕微諧波時(shí)EMD會(huì)把基波分量分裂成相鄰兩個(gè)IMF選哪個(gè)都不完全對(duì)相量估值的波形出現(xiàn)周期性毛刺。后來(lái)用EEMD解決了但集合平均導(dǎo)致計(jì)算量暴漲。我的經(jīng)驗(yàn)是先用默認(rèn)emd跑一遍看IMF的瞬時(shí)頻率是否在預(yù)期頻帶附近如果有明顯模態(tài)混疊再考慮EEMD把噪聲幅值設(shè)為信號(hào)標(biāo)準(zhǔn)差的0.1到0.2倍集合次數(shù)50到100次效果基本穩(wěn)定。上來(lái)就開(kāi)默認(rèn)參數(shù)跑EMD很容易得到一堆莫名其妙的偽IMF。6.2 hilbert函數(shù)對(duì)矩陣數(shù)據(jù)的坑Matlab的hilbert函數(shù)在處理矩陣時(shí)是按列運(yùn)算的。如果你把多個(gè)信號(hào)寫(xiě)成行向量的矩陣傳進(jìn)去拿到的解析信號(hào)完全不是你想要的。我在做多通道同步相量對(duì)比時(shí)踩過(guò)這個(gè)坑輸出結(jié)果相位全亂了。穩(wěn)妥做法是對(duì)每列信號(hào)單獨(dú)調(diào)hilbert或者在寫(xiě)代碼時(shí)明確用列向量例如zs hilbert(x(:)); % 強(qiáng)制按列這個(gè)細(xì)節(jié)極不起眼但出了問(wèn)題很難排查因?yàn)榉悼雌饋?lái)還是對(duì)的只有相位不對(duì)非常具有迷惑性。6.3 相角unwrap的跳變問(wèn)題對(duì)連續(xù)時(shí)刻的瞬時(shí)相角序列做unwrap是常規(guī)操作但噪聲大了以后unwrap可能會(huì)在瞬時(shí)頻率接近±Fs/2的毛刺處產(chǎn)生假跳變。更好的做法是分兩步先對(duì)瞬時(shí)頻率做約束濾波比如限制在45Hz~55Hz之外的值直接剔除再用濾波后的頻率去累積積分相位。換句話說(shuō)不要直接拿angle的輸出去做unwrap相位展開(kāi)要和瞬時(shí)頻率聯(lián)動(dòng)校驗(yàn)。我后來(lái)寫(xiě)了一個(gè)小工具函數(shù)輸入解析信號(hào)輸出經(jīng)過(guò)頻差校驗(yàn)的瞬時(shí)相位序列這才把HHT的輸出穩(wěn)定性拉上來(lái)。6.4 不同Matlab版本的函數(shù)差異代碼遷移時(shí)最容易翻車(chē)的是cwt和emd這兩個(gè)函數(shù)。老版本Matlab的cwt寫(xiě)法是cwt(x, scales, morl)需要自己定義尺度向量返回的是尺度軸而不是頻率軸R2020b之后推薦寫(xiě)成cwt(x, Fs)默認(rèn)Morse小波直接返回頻率。emd函數(shù)則是R2018a之后才內(nèi)置的之前要用第三方工具包。如果直接拿網(wǎng)上舊代碼跑不報(bào)錯(cuò)算運(yùn)氣好報(bào)錯(cuò)了八成是這里。建議先跑一下ver(wavelet)確認(rèn)工具箱版本再用新語(yǔ)法改寫(xiě)。6.5 對(duì)比實(shí)驗(yàn)的數(shù)據(jù)長(zhǎng)度對(duì)齊細(xì)節(jié)最后提醒一點(diǎn)四種方法對(duì)數(shù)據(jù)長(zhǎng)度的需求完全不同做對(duì)比實(shí)驗(yàn)時(shí)不要只統(tǒng)一“采樣點(diǎn)數(shù)”。FFT需要窗長(zhǎng)對(duì)應(yīng)整數(shù)倍工頻周期CWT在數(shù)據(jù)邊緣表現(xiàn)差但內(nèi)部穩(wěn)定EMD對(duì)數(shù)據(jù)長(zhǎng)度和極值點(diǎn)分布都敏感。我建議統(tǒng)一以“分析窗覆蓋的工頻周期數(shù)”為對(duì)齊標(biāo)準(zhǔn)然后把每種方法各自最穩(wěn)定的中心段取出來(lái)評(píng)估。這樣得到的對(duì)比結(jié)果才是方法本身的差異而不是數(shù)據(jù)長(zhǎng)度的差異。從做這個(gè)課題的整個(gè)過(guò)程來(lái)看最讓我意外的不是某個(gè)算法多強(qiáng)而是它們各自的弱點(diǎn)多么“各有特色”。同步相量估計(jì)這個(gè)看似成熟的方向深入進(jìn)去之后每層都有讓人撓頭的問(wèn)題。但反過(guò)來(lái)想正因?yàn)闆](méi)有銀彈幾種方法相互印證、按場(chǎng)景切換才是在實(shí)際工程中真正可行的路。