imf瞬时频率跳变问题

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

imf瞬时频率跳变问题相关推荐

气象数据分析之EMD方法介绍及python的实现
文章目录前言一.EMD方法介绍二.在python中的实现 1.引入库 2.生成一个随机的信号 3.做EMD分解,提取IMF和res 4.可视化最后前言经验模态分解(Empirical Mo ...
离散信号的希尔伯特变换的计算公式_希尔伯特变换和瞬时频率问题--连载（二）...
写在开始的一段话: PS:OK,上一期关于希尔伯特变换的文章发出后,有知友在评论区说"看到最后--居然这--",哈哈,其实我也挺愧疚大家的,明明一篇知识分享的文章,却写到结尾都没进 ...
希尔伯特-黄变换（HHT）的前世今生——一个从瞬时频率讲起的故事
一.来自小X的疑问从前有一个国家,叫做实国(一维国度),里边有个叫小X的小人儿. 小X是一根线段,他每天最爱做的事情就是跳舞. 因为小X的舞姿十分稳定,同伴们都说他的头部跳出了频率为1hz的正弦曲线 ...
希尔伯特谱、边际谱、包络谱、瞬时频率/幅值/相位——Hilbert分析衍生方法及MATLAB实现
上一篇文章对希尔伯特-黄变换(HHT)的前世今生进行了介绍. 不过在研究中通常并不是到希尔伯特-黄变换就停止了. 而是要用到诸如希尔伯特谱.包络谱.边际谱.瞬时频率/幅值/相位等方法进一步分析. 这些 ...
hs = hht(imf)的输出介绍
这是MATLAB官方网站对它的介绍 hs = hht(imf) 返回hs由固有模式函数指定的信号的希尔伯特频谱imf.hs对于分析包含频谱含量随时间变化的信号混合的信号很有用.用于hht对信号执行希尔 ...
[HT/NHT/DQ]-三种基于EMD的瞬时频率计算方法的比较
文章目录 0 前言 1 瞬时频率和经验模态分解 1.1 瞬时频率的定义 1.2 经验模态分解 2 瞬时频率的计算方法 2.1 HT方法 2.2 NHT方法 2.3 DQ方法 3 三种瞬时频率计算方法的 ...
负值瞬时频率出现原因及解决方法总结
文章目录 0 英文缩写 1 负值瞬时频率出现的原因 2 三次样条插值 3 IMF与EMD 4 经验AM-FM分解 4.1 直接正交 4.2 归一化希尔伯特变换 5 总结如果觉得乱,可以直接看最后的[ ...
“类EMD”算法分解后要怎样使用（1）——内涵模态分量IMF的方差贡献率、平均周期、相关系数的计算及MATLAB代码实现
之前我们有了十几篇文章讲述了EMD算法的基础理论.IMF的含义.EMD的MATLAB实现方法,EEMD.CEEMD.CEEMDAN.VMD.ICEEMDAN.LMD.EWT的理论及代码实现,还讲到了H ...
Google Instant 瞬时搜索上手指南
Google Instant是Google刚刚发布的一种新的搜索方式,随着你在搜索框里输入文字,Google将同时给出搜索结果,同时在搜索框里还会根据你输入的关键字给出搜索建议,通过上下键即可切换.按 ...

imf瞬时频率跳变问题

imf瞬时频率跳变问题相关推荐

最新文章

热门文章