本發(fā)明涉及諧波檢測領(lǐng)域，尤其涉及一種基于CEEMD算法和希爾伯特變換的諧波分析方法。

背景技術(shù)：

現(xiàn)在電力電子器件使用的越來越多，導(dǎo)致電力系統(tǒng)中諧波也越來越多，危害系統(tǒng)和各類電器的正常穩(wěn)定運行。

目前的諧波分析方法一般采用傳統(tǒng)經(jīng)驗?zāi)B(tài)分解(Empirical Mode Decomposition，EMD)方法，該方法是由美國國家宇航局于20世紀(jì)末提出，被認(rèn)為是對以傅里葉變換為基礎(chǔ)的線性、穩(wěn)定頻譜分析的一個重大突破。EMD利用信號的局部特征時間尺度，從原信號中提取出若干階固有模態(tài)函數(shù)(Intrinsic Mode Function，IMF)和一個殘余量，分解出的各階IMF分量突出了數(shù)據(jù)的局部特征，殘余分量體現(xiàn)了信號中的緩慢變化量。對各個IMF進行分析，可以更準(zhǔn)確有效地把握原數(shù)據(jù)的特征信息。每階IMF必須滿足以下2個條件：

(1)整個IMF信號中過零點與極值點個數(shù)相等或至多相差1個；

(2)由信號上局部極大值點確定的包絡(luò)線和由局部極小值點確定的包絡(luò)線的均值均為零，即信號關(guān)于時間軸局部對稱。

但該方法的分解效果一般，不能很好地滿足對諧波分析的需求。

技術(shù)實現(xiàn)要素：

針對上述現(xiàn)有技術(shù)中的不足，本發(fā)明提供一種基于CEEMD算法(互補集合經(jīng)驗?zāi)B(tài)分解算法)和希爾伯特變換的諧波分析方法，具有分析效果佳、處理效果好的優(yōu)點。

為了實現(xiàn)上述目的，本發(fā)明提供一種基于CEEMD算法和希爾伯特變換的諧波分析方法，包括步驟：

S1：對一原始信號S(t)進行互補集合經(jīng)驗?zāi)B(tài)分解，獲得一本征模態(tài)函數(shù)和一余量R(t)；其中m為所述固有模態(tài)函數(shù)分量的個數(shù)，C_r(t)為第r個本征模態(tài)函數(shù)的數(shù)列，r為大于等于1的自然數(shù)，t為時間；

S2：對所述本征模態(tài)函數(shù)進行希爾伯特變換，獲得所述本征模態(tài)函數(shù)的希爾伯特變換表達式；

S3：根據(jù)所述本征模態(tài)函數(shù)和所述本征模態(tài)函數(shù)的希爾伯特變換表達式計算獲得所述原始信號的瞬時幅值和瞬時相位；

S4：根據(jù)所述瞬時相位計算獲得所述原始信號的瞬時頻率和瞬時角速度。

本發(fā)明的進一步改進在于，所述本征模態(tài)函數(shù)的希爾伯特變換表達式為：

其中，P為柯西主值，τ為時間窗平移因子。

本發(fā)明的進一步改進在于，所述原始信號的瞬時幅值的表達式為：

本發(fā)明的進一步改進在于，所述原始信號的瞬時相位的表達式為：

本發(fā)明的進一步改進在于，所述原始信號的瞬時頻率的表達式為：

本發(fā)明的進一步改進在于，所述原始信號的瞬時角速度的表達式為：

本發(fā)明由于采用了以上技術(shù)方案，使其具有以下有益效果：

通過互補集合經(jīng)驗?zāi)B(tài)分解的采用，提高了信號的分析處理效果。通過對本征模態(tài)函數(shù)進行希爾伯特變換，實現(xiàn)可得到該信號的希爾伯特譜、時頻能量譜等，達到對信號進行分析的目的。

附圖說明

圖1為本發(fā)明實施例的基于CEEMD算法和希爾伯特變換的諧波分析方法的流程圖；

圖2為本發(fā)明實施例的一工頻標(biāo)準(zhǔn)信號的波形圖；

圖3為本發(fā)明實施例的第一諧波信號的波形圖；

圖4為本發(fā)明實施例的第二諧波信號的波形圖；

圖5為本發(fā)明實施例的目標(biāo)信號的波形圖；

圖6為本發(fā)明實施例的對目標(biāo)信號進行總體平均經(jīng)驗?zāi)B(tài)分解的分量波形對比圖；

圖7為圖6中第一分量的時頻特征曲線圖；

圖8為圖6中第一分量的時幅特征曲線圖；

圖9為圖6中第二分量的時頻特征曲線圖；

圖10為圖6中第二分量的時幅特征曲線圖；

圖11為圖6中第三分量的時頻特征曲線圖；

圖12為圖6中第三分量的時幅特征曲線圖；

圖13為本發(fā)明實施例的希爾伯特譜圖。

具體實施方式

下面根據(jù)附圖1-13，給出本發(fā)明的較佳實施例，并予以詳細描述，使能更好地理解本發(fā)明的功能、特點。

請參閱圖1，本發(fā)明實施例的一種基于CEEMD算法和希爾伯特變換的諧波分析方法，包括步驟：

S1：對一原始信號S(t)進行互補集合經(jīng)驗?zāi)B(tài)分解，獲得一本征模態(tài)函數(shù)和一余量R(t)；其中m為固有模態(tài)函數(shù)分量的個數(shù)，C_r(t)為第r個本征模態(tài)函數(shù)的數(shù)列，r為大于等于1的自然數(shù)，t為時間。

互補集合經(jīng)驗?zāi)B(tài)分解可歸納為以下幾個步驟:

(1)在原始信號S(t)中加入模值相等的正負(fù)兩組白噪聲信號w(t)與-w(t),獲得加入了正、負(fù)白噪聲后的信號S1(t)、S2(t)：

S1(t)＝S(t)+w(t)；

S2(t)＝S(t)-w(t)；

(2)對信號S1(t)，S2(t)分別進行傳統(tǒng)經(jīng)驗?zāi)B(tài)分解(EMD),獲得分解后的本征模態(tài)函數(shù)的數(shù)列C_i+(t)與C_i-(t)；

(3)添加不同的白噪聲重復(fù)步驟(1)、(2)N次，并進行集總平均,

則原始信號(r＝1，2，……，m)。

S2：對本征模態(tài)函數(shù)進行希爾伯特變換，獲得本征模態(tài)函數(shù)的希爾伯特變換表達式

其中，P為柯西主值，τ為時間窗平移因子。

S3：根據(jù)本征模態(tài)函數(shù)和本征模態(tài)函數(shù)的希爾伯特變換表達式計算獲得原始信號的瞬時幅值和瞬時相位；

其中，原始信號的瞬時幅值的表達式為：

原始信號的瞬時相位的表達式為：

S4：根據(jù)瞬時相位計算獲得原始信號的瞬時頻率和瞬時角速度。

其中，原始信號的瞬時頻率的表達式為：

原始信號的瞬時角速度的表達式為：

為了驗證本方法的效果，通過模擬獲得一初始信號f(t)，初始信號f(t)可通過在工頻標(biāo)準(zhǔn)信號的基礎(chǔ)上添加第一諧波信號和第二諧波信號獲得，工頻標(biāo)準(zhǔn)信號、第一諧波信號、第二諧波信號和初始信號的波形請參見附圖2～圖5。

f(t)＝18cos(2π×10×t)+60cos(2π×50×t)+10cos(2π×150×t)；通過本發(fā)明的基于CEEMD算法和希爾伯特變換的諧波分析方法對其做分析處理，可獲得其總體平均經(jīng)驗?zāi)B(tài)分解的分量波形對比圖，請參閱圖6，圖中IMF1為第一分量、IMF2為第二分量、IMF3為第三分量、r為余量。其中，第一分量的時頻特征曲線請參閱圖7，第一分量的時幅特征曲線請參閱圖8，第二分量的時頻特征曲線請參閱圖9，第二分量的時幅特征曲線請參閱圖10，第三分量的時頻特征曲線請參閱圖11，第三分量的時幅特征曲線請參閱圖12，希爾伯特譜請參閱圖13。

以上結(jié)合附圖實施例對本發(fā)明進行了詳細說明，本領(lǐng)域中普通技術(shù)人員可根據(jù)上述說明對本發(fā)明做出種種變化例。因而，實施例中的某些細節(jié)不應(yīng)構(gòu)成對本發(fā)明的限定，本發(fā)明將以所附權(quán)利要求書界定的范圍作為本發(fā)明的保護范圍。