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

背景技術(shù)：

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

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

(1)整個(gè)IMF信號(hào)中過(guò)零點(diǎn)與極值點(diǎn)個(gè)數(shù)相等或至多相差1個(gè)；

(2)由信號(hào)上局部極大值點(diǎn)確定的包絡(luò)線(xiàn)和由局部極小值點(diǎn)確定的包絡(luò)線(xiàn)的均值均為零，即信號(hào)關(guān)于時(shí)間軸局部對(duì)稱(chēng)。

但該方法的分解效果一般，不能很好地滿(mǎn)足對(duì)諧波分析的需求。

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

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

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

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

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

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

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

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

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

本發(fā)明的進(jìn)一步改進(jìn)在于，所述原始信號(hào)的瞬時(shí)幅值的表達(dá)式為：

本發(fā)明的進(jìn)一步改進(jìn)在于，所述原始信號(hào)的瞬時(shí)相位的表達(dá)式為：

本發(fā)明的進(jìn)一步改進(jìn)在于，所述原始信號(hào)的瞬時(shí)頻率的表達(dá)式為：

本發(fā)明的進(jìn)一步改進(jìn)在于，所述原始信號(hào)的瞬時(shí)角速度的表達(dá)式為：

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

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

附圖說(shuō)明

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

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

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

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

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

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

圖7為圖6中第一分量的時(shí)頻特征曲線(xiàn)圖；

圖8為圖6中第一分量的時(shí)幅特征曲線(xiàn)圖；

圖9為圖6中第二分量的時(shí)頻特征曲線(xiàn)圖；

圖10為圖6中第二分量的時(shí)幅特征曲線(xiàn)圖；

圖11為圖6中第三分量的時(shí)頻特征曲線(xiàn)圖；

圖12為圖6中第三分量的時(shí)幅特征曲線(xiàn)圖；

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

具體實(shí)施方式

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

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

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

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

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

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

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

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

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

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

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

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

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

其中，原始信號(hào)的瞬時(shí)幅值的表達(dá)式為：

原始信號(hào)的瞬時(shí)相位的表達(dá)式為：

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

其中，原始信號(hào)的瞬時(shí)頻率的表達(dá)式為：

原始信號(hào)的瞬時(shí)角速度的表達(dá)式為：

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

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

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