问题是这样被发现和引起的：我打算求出各信号各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化”处理之后，也会得到相应的改变。具体情况如何，以后再研究。

 

 

（本文首发于：http://blog.sina.com.cn/s/blog_6ad0d3de0100zxfq.html

首发时间：2011-11-03 13:57:00）