問題是這樣被發現和引起的：我打算求出各信號各IMF分量的瞬時頻率的一些特征參數（尤其是平均頻率），以便比較這些分量是不是同一物理量引起的。先求TW486012各IMF分量瞬時頻率f_imf_TW486012的平均頻率。根據時頻分析理論，信號譜的平均頻率等于信號瞬時頻率的時間平均mtf_imf_TW486012。

?

mtf_imf_TW486012=sum(f_imf_TW486012.*A_imf_TW486012.^2)'./sum(A_imf_TW486012.^2)';

plot(mtf_imf_TW486012)?

?

????見圖18-1藍線部分。橫軸是IMF分量序號。?


????圖18-1 imf_TW486012瞬時頻率的平均頻率
?

????（附帶說一下：對于這種平均頻率的計算方式，為什么一定要按照信號能量的分布密度來進行加權計算，我不是太明白。因為瞬時頻率與瞬時幅值是兩個獨立的變量，求一個變量的時間平均卻一定要與另一個變量扯上關系，邏輯上好象說不通。雖然理論上可以證明，此時間平均頻率等于信號總譜的平均頻率。當只能得到信號總譜時，平均頻率用信號能量的分布密度來表征是可以理解的。但已經得到了瞬時頻率，就沒有必要考慮與能量分布的關系了。）

????另外再計算一組純粹由瞬時頻率函數本身確定的平均頻率。為了以示區別，就把它叫做瞬時頻率的算術平均吧。程序如下：

?

mf_imf_TW486012=mean(f_imf_TW486012)';

hold on

plot(mf_imf_TW486012,'r-')

?

????曲線見圖18-1紅線部分。?

?

????它們不是單調下降的。這很奇怪哦。什么原因使得低頻分量的平均頻率反而上升了呢？回過頭去再稍微看一下 圖14-2 f_imf_TW486012，很容易注意到圖中象“譜線”一樣的跳變的“高頻”成分。放大曲線，可以看到并不全是“譜線”，有些高頻成分還存在著一個區間。對，就是這些跳變的“高頻”成分拉高了它們的平均頻率！它們是不符合工程實際情況的。實際的瞬時頻率，變化應該是比較連續、緩慢的。

????必須得找出原因，解決這個問題。網上搜索，沒看到相關的文獻、帖子；把它發到論壇里去討論、請教，一直都沒有人回復。不知道他們是不屑于我這個問題呢，還是看不懂我這個問題。沒辦法，只好自己慢慢研究、琢磨。?

????IMF分量瞬時頻率的計算，主要用到了兩個工具箱里的兩個函數：EMD工具箱里的hhspectrum函數、時頻工具箱里的instfreq函數。hhspectrum函數的主要作用就是求出IMF分量的解析信號，然后根據此解析信號調用instfreq函數，求出其瞬時頻率。

?

????先用sptool放大看一看f_imf_TW486012頻率跳變點附近的情況。第5個分量的一個頻率跳變點：

????圖18-2 f_imf_TW486012第5個分量頻率跳變點位置

?

????再看看解析信號對應于頻率跳變點附近的情況：

?

hil_imf_TW486012=hilbert(imf_TW486012');

?

????復信號裝入sptool之后，可以直接查看它的實部、虛部、幅值與相位函數。

?

????圖18-3 imf_TW486012第5個分量（亦即解析信號的實部）對應于其頻率跳變點位置

?

????圖18-4 imf_TW486012第5個分量的解析信號的虛部對應于其頻率跳變點位置

???

????咋一看，這虛部信號上下包絡線對稱，局部均值似等于0，好象是一基本模式函數IMF。但再看看圖中標尺所示的兩個局部極大值點，其值都小于0，說明虛部信號并不是真正的基本模式函數IMF。

圖18-5 hil_imf_TW486012第5個分量第45100_45143點包括頻率跳變點的一段數據

?

????通過數據可以更加清楚地看出解析信號在頻率跳變點處的變化情況。那么，虛部信號局部極大值小于0是否必定會發生頻率跳變呢？

?

????圖18-6 imf_TW486012第4個分量解析信號虛部某小于0的局部極大值?

?

????圖18-7 上述小于0的局部極大值對應的瞬時頻率位置沒有發生頻率跳變

?

?

?

????圖18-8 imf_TW486012第5個分量的解析信號的相位函數對應于其頻率跳變點位置

?

????再看imf_TW486012第5個分量解析信號的相位函數。可以看出，在頻率跳變點處，相位函數發生了遞減的情況。由于瞬時頻率等于相位函數的導數，這說明此處出現了負頻率。如果瞬時頻率確定是正頻率的話，那么除了在＋pi變為－pi處之外，相位函數應該是單增函數才對。

????圖18-9 f_imf_TW486012第4個分量頻率跳變點位置

?

????虛部局部極大（小）值小（大）于0處，瞬時頻率不一定會發生跳變，那么瞬時頻率發生跳變的地方，是不是一定會發生虛部局部極大（小）值小（大）于0的情況呢？圖圖18-9為第4個分量一個頻率跳變點。

????圖18-10 imf_TW486012第4個分量（亦即解析信號的實部）對應于其頻率跳變點位置

?

????圖18-11 imf_TW486012第4個分量的解析信號的虛部對應于其頻率跳變點位置

?

????可以看出，并沒有發生虛部信號局部極大（小）值小（大）于0的情況。

圖18-12 hil_imf_TW486012第4個分量的解析信號第26476_26500點包括頻率跳變點的一段數據

?

????從具體的數據可以更清楚地了解瞬時頻率跳變點處的情況。

????圖18-13 imf_TW486012第4個分量的解析信號的相位對應于其頻率跳變點位置

?

????相位函數仍然存在遞減區間。?

????圖18-14 f_imf_TW486012第14個分量頻率跳變點位置

?

????再看f_imf_TW486012第14個分量頻率跳變點位置。?

????圖18-15 imf_TW486012第14個分量的解析信號實部對應于其頻率跳變點位置

?

????實部。

?

????圖18-16 imf_TW486012第14個分量的解析信號虛部對應于其頻率跳變點位置

?

???虛部。

?

????圖18-17 imf_TW486012第14個分量的解析信號的相位對應于其頻率跳變點位置

?

????相位。

?

????通過上面的分析，除了相位遞減與頻率跳變存在必然關系外，似乎很難再找出其它的必然關系。

????上面幾個頻率跳變點，還有一個共同點，那就是其實部、虛部數據都極度接近0。那么是不是幅值靠近0的地方就必定會發生頻率跳變呢？?將信號幅值按升序排列，同時讓瞬時頻率跟著幅值重新排列，就可以看出上面問題的答案。 ?

??

Af=[A_imf_TW486012(4,:)',f_imf_TW486012(4,:)'];

srAf=sortrows(Af);

subplot(1,2,1)

plot(srAf(:,1))

subplot(1,2,2)

plot(srAf(:,2))

?

????運行，得：

????圖18-18 imf_TW486012第4個分量解析信號幅值按升序排列

?

????咋一看，頻率跳變點好象的確集中在幅值靠近0的地方。但將上面右圖的左側部分橫軸放大，就可以看出，幅值靠近0的地方并不一定就是頻率跳變點。見圖18-19。?

?

????圖18-19 imf_TW486012第4個分量瞬時頻率跟隨升序幅值重新排列

?

????招數都用盡了，還是找不出一個扼要的方法解決問題。最后一招，就是頭痛醫頭腳痛醫腳的笨辦法，根據頻率跳變處的波形特點，編長程序消除此頻率跳變。

?

????程序如下：

?

g=0.071893186;
d=0.33;
h=0.35;
nf_imf_TW486012=f_imf_TW486012;
M=cell(size(nf_imf_TW486012,1)-1,2);
df_imf_TW486012=diff(f_imf_TW486012)';
for i=1:size(nf_imf_TW486012,1)-1
a=1;b=1;
??for j=2:size(nf_imf_TW486012,2)-1
????if df_imf_TW486012(i,j-1)*df_imf_TW486012(i,j)<0,
????????if nf_imf_TW486012(i,j)<g,
????????????if df_imf_TW486012(i,j)>d,??
???????????????M{i,1}(a)=j;
???????????????a=a+1;
????????????elseif df_imf_TW486012(i,j-1)<-d,
???????????????M{i,2}(b)=j;
???????????????b=b+1;
????????????end
????????end
????end
??end
??for k=1:length(M{i,1})
????for l=1:length(M{i,2})
?????????mnf=mean(nf_imf_TW486012(i,M{i,1}(k)+1:M{i,2}(l)-1));
??????if mnf>h
?????????nf_imf_TW486012(i,M{i,1}(k)+1:M{i,2}(l)-1)=nf_imf_TW486012(i,M{i,1}(k)+1:M{i,2}(l)-1)-0.5;
??????end
????end
??end
end

%上面只解決了數據中段（從第3個（含）至倒數第3個（含）數據之間的

%頻率跳變問題。對于從第1、第2個數據就已經處于高頻位，以及在倒數

%第2、倒數第1個數據都未回歸低頻位的跳頻問題，需下面程序部分解決。

h1_end=0.48;
d1_end=0.45;
nf_imf_TW486012=nf_imf_TW486012';%行變列
dnf_imf_TW486012=diff(nf_imf_TW486012);
[Cmin,Imin]=min(dnf_imf_TW486012);
[Cmax,Imax]=max(dnf_imf_TW486012);
for i=1:size(nf_imf_TW486012,2)-1
if Cmin(i)<-d1_end
????if mean(nf_imf_TW486012(1:Imin(i),i))>h1_end
???????nf_imf_TW486012(1:Imin(i),i)=nf_imf_TW486012(1:Imin(i),i)-0.5;
????end
end??
if Cmax(i)>d1_end
????if mean(nf_imf_TW486012(Imax(i)+1:end,i))>h1_end
???????nf_imf_TW486012(Imax(i)+1:end,i)=nf_imf_TW486012(Imax(i)+1:end,i)-0.5;
????end
end???
end

for k=1:14
subplot(5,3,k)
plot(nf_imf_TW486012(:,k))
end

?

????運行，得：?

????圖18-20 f_imf_TW486012消除頻率跳變現象之后的瞬時頻率nf_imf_TW486012

?

????可以看到都沒有頻率跳變了。再看看圖18-2與圖18-9兩處頻率跳變點局部情況：

????圖18-21 圖18-2消除跳變之后的局部波形

?

????圖18-22 圖18-9消除跳變之后的局部波形

?

????曲線變得連續、光滑了。

?

????問題雖然解決了，但是用這種笨辦法解決，心中總覺不甘。再研究研究。

????先把消除頻率跳變現象之后的瞬時頻率nf_imf_TW486012的頻率跳變點數據序號都找出來。將nf_imf_TW486012的局部極小值都找出來，其中最小的那幾個就是頻率跳變點。

?

nf4_imf_TW486012=nf_imf_TW486012(:,4);

nf4=-nf4_imf_TW486012;

[nf4c,nf4cn]=findpeaks(nf4,'sortstr','ascend');

nf4c=fliplr(nf4c);%跳頻點放到序列最前面

nf4cn=fliplr(nf4cn);

?

?%找出解析信號虛部負極大值與正極小值序號

ihil_imf_TW486012=imag(hil_imf_TW486012);
ihtw4=ihil_imf_TW486012(:,4);
[ihtw4p,ihtw4pn]=findpeaks(ihtw4,'sortstr','ascend');

ihtw4=-ihtw4;
[ihtw4c,ihtw4cn]=findpeaks(ihtw4,'sortstr','ascend');

ihtw4c=-ihtw4c;

nf4c=nf4c';nf4cn=nf4cn';ihtw4p=ihtw4p';ihtw4pn=ihtw4pn';ihtw4c=ihtw4c';ihtw4cn=ihtw4cn';%行變列

?

????將上面6列數據放在一起觀察，見圖18-23。

????圖18-23

?

????圖中第1、2列前5個點是頻率跳變點；第3、4列前11個點是虛部負極大值點；第5、6列前10個點是虛部正極小值點。可以看出，在5個頻率跳變點中，有3個是信號虛部負極大值點，且最靠近0處；有1個是信號虛部正極小值點，且與0值只隔了一個數據。

?

????初步結論：虛部存在負極大值點與正極小值點，且它們極度靠近0，是發生頻率跳變的主要原因。

????但上面的數據，有兩處例外。先看看最右邊一列粉色方框標示的、序號為48568的數據點是怎么回事。它是虛部正極小值點，且極度靠近0。但它沒有被列入頻率跳變點。用sptool放大觀察原信號局部，見圖18-24。

?

????圖18-24

?

????原來這也是一個堪稱頻率跳變的點。因為它的波形，不符合我程序中的搜索條件

if df_imf_TW486012(i,j-1)*df_imf_TW486012(i,j)<0，所以沒有被我的程序列為頻率跳變點。但它的高度也超過了0.5的一半0.25。從圖18-23數據的“同類相聚”性來說，它無疑也是一個頻率跳變點。需要修改的是我的搜索條件。

?

???再看看圖18-23中的第1個頻率跳變點，nf_imf_TW486012中序號為26493的點。它在虛部極大值點中緊挨著0，但是是正的，不符合我的“初步結論”的條件。可以看出，這個數據點是個罕有的例子（緊挨著0的正極大值），碰巧被我在圖18-11中引用到。所以不能援引它來說明“在負（正）極大（小）值點之外也有很多地方發生頻率跳變”。

?

????現在把瞬時頻率計算過程中，包含這個數據點的一段數據，調出來看一看。?可以分步了解跳變頻率是怎樣產生的。instfreq函數計算瞬時頻率的主要語句是fnormhat=0.5*(angle(-x(t+1).*conj(x(t-1)))+pi)/(2*pi)。所以要看這項計算中具體數據的變化情況。

?

hiltw4=hil_imf_TW486012(:,4);

tt=2:size(hil_imf_TW486012,1)-1;

x=-hiltw4(tt+1).*conj(hiltw4(tt-1));

angx=angle(x);

f_angx=0.5*angx/(2*pi);

%f_imf_TW486012=0.25+f_angx?

?

????運行，打開hiltw4、x、f_angx變量，將包含序號26493數據點的一段數據截圖。見下。

?????圖18-25

?

????x、f_angx的長度比hiltw4要短2個點，所以x、f_angx的26494、就是hiltw4的26495。注意這是原解析信號hiltw4實部負極小值點，而它的虛部正極大值點是26492。

?

????還可以將這一小段數據作圖，可以更直觀。

?

plot(hiltw4(26490:26496),'-bo')

grid

hold on

plot(x(26489:26495),'-g*')

plot(10.*x(26489:26495),'-r*')%x(26489:26495)變化范圍太小，看不清楚。擴大10倍看。

?

?????圖18-26

?

????上圖橫軸是實部，縱軸是虛部。可以看出，這一小段數據，解析信號（右邊藍線）首先在第4象限（26490點），然后跳到第1象限，再回到第4象限，再跳到第3象限，最后再回到第4象限。將此軌跡看成一個閉合曲線的話，原點正好被排除在此閉合曲線之外。這就是發生頻率跳變的根本原因。至于x=-hiltw4(tt+1).*conj(hiltw4(tt-1))，它只是解析信號的一個映射。當解析信號上面的特征具備后，x產生跳變就是必然的了。圖中左邊紅（綠）線首先在第3象限（26489點），然后跳入第1象限（26492、26493點），再跳入第2象限（26494點），并在此形成最大相位角（注意此曲線各點相位角就相當于各點瞬時頻率），最后回到第3象限。整個軌跡避開了第4象限。

?

????通過上面的分析，現在基本上可以得出如下的結論：?當利用信號IMF分量的解析信號對分量進行瞬時頻率估計時，需要對解析信號的虛部進行“IMF化”處理，使其也成為一個IMF基本模式函數，并與原分量構成一新的復信號。它們的過0點需要按下面的順序排列：實部降0點－>虛部降0點－>實部升0點－>虛部升0點－>實部降0點……。當復信號模值較大時，此條件基本上可自然滿足；但當模值極靠近0時，需要對信號特別處理才可以滿足。此處“降0點”“升0點”表示信號從局部極大值降到局部極小值、局部極小值升到局部極大值時所經過的0點。

?

????對信號作這樣的處理，是基于信號處理的結果，要盡量符合現實情況的要求。瞬時頻率應該是連續的，光滑的。這種思想與最初建立EMD分解、HHT變換的思想是一致的。經典的Hilbert變換，是從信號的全局處理角度提出的，因此其構建的目標函數是復信號的頻譜沒有負頻率成分。現在進入到時頻分析時代，這一目標函數也應該改變了。

?

????上面的結論主要來自于對信號數據中段（從第3個（含）至倒數第3個（含）數據之間）的

頻率跳變點的分析。至于數據兩端的頻率跳變，我相信信號虛部經過“IMF化”處理之后，也會得到相應的改變。具體情況如何，以后再研究。

與50位技術專家面對面20年技術見證，附贈技術全景圖

總結

以上是生活随笔為你收集整理的imf瞬时频率跳变问题的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。