求出HHT频率为何不匹配?

shenlan888

我的程序如下:
clear
fs=200;
t=0:1/fs:1;
s=cos(10*pi*t)+2*cos(40*pi*t);
figure(1)
plot(t,s)
c=emd(s);%进行emd分解
figure(2)
b=size(c,1);
for j=1:1:b
    subplot(b+1,1,j)
    plot(t,c(j,:))
end
subplot(b+1,1,b+1)
plot(t,sum(c))
% c(1,:)=cos(10*pi*t);
% c(2,:)=2*cos(40*pi*t);
l=length(s);
t=t(2:end-1);
figure(3)
% b=1;
for m=1:b
    H(m,:)=hilbert(c(m,:));
%     x(m,:)=real(H(m,:));
%     y(m,:)=imag(H(m,:));
%     fi(m,:)=atan(y(m,:)./x(m,:));%瞬时相位
%     g(m,:)=szwf(fi(m,:),t);%求微分,求瞬时角频率
%     f(m,:)=g(m,:)./(2*pi);%瞬时频率
    f(m,:)=instfreq(H(m,:)');
    f1= f(m,:)*fs;
    subplot(b,1,m);
    plot(t,f1);
    xlabel('t(s)');
    ylabel('freq(hz)')
    title('时频图');
end
求出的频率为什么和采样频率有关,而且和原始信号不匹配?
但是下面这个程序求的信号就是正确的,希望得到大家的帮助
clear
T=10;
N=1000;
n=0:1:N-1;
dt=T/N;%采样周期
t=n*dt;
x=2*cos(40*pi*t);
y=hilbert(x');
[fnor,t]=instfreq(y);
fnor1=fnor*N/T;
fnor1(1,1)
plot(t,fnor1);
这个频率是正确的

wuqiong_cea

黄给出了一个求解瞬时频率的程序，而他本人也关于瞬时频率问题写了专门的文章。下面给出他的程序代码：他的求解方法与instfreq不太相同
% this is a function to calculate instantaneous based on EMD method
%
%   function omega = ifndq(vimf, dt)
%
% INPUT:   
%          vimf:        an IMF;
%          dt:          time interval of the imputted data
% OUTPUT:
%          omega:       instantanesous frequency, which is 2*PI/T, where T
%                       is the period of an ascillation
%
% References can be found in the "Reference" section.
%
% The code is prepared by Zhaohua Wu. For questions, please read the "Q&A" section or
% contact
%   zhwu@cola.iges.org
%

function [omega,amplitude,error] = ifndq(vimf, dt)
Nnormal=5;
rangetop=0.90;

vlength = max( size(vimf) );
vlength_1 = vlength -1;

for i=1:vlength,
    abs_vimf(i)=vimf(i);
    if abs_vimf(i) < 0
        abs_vimf(i)=-vimf(i);
    end
end

for jj=1:Nnormal,%iterate Normal steps,typicall is 2 or 3,here choose 5
    [spmax, spmin, flag]=extrema(abs_vimf);
    dd=1:1:vlength;
    upper= spline(spmax(:,1),spmax(:,2),dd);%e1(t)used to noumalized the data

    for i=1:vlength,
        abs_vimf(i)=abs_vimf(i)/upper(i);%normalized
    end
end
%the abs_vimf is the normalized signal which is called envelop
for i=1:vlength,
    nvimf(i)=abs_vimf(i);
    if vimf(i) < 0;
        nvimf(i)=-abs_vimf(i);
    end
end
amplitude = vimf./nvimf';
error = (abs(hilbert(nvimf))-1).^2;
for i=1:vlength,
    dq(i)=sqrt(1-nvimf(i)*nvimf(i));%the direct quadrature method is used to find the phase function
end
for i=2:vlength_1,
    devi(i)=nvimf(i+1)-nvimf(i-1);%
    if devi(i)>0 & nvimf(i)<1
        dq(i)=-dq(i);
    end
end

rangebot=-rangetop;     
for i=2:(vlength-1),
    if nvimf(i)>rangebot & nvimf(i) < rangetop
        omgcos(i)=abs(nvimf(i+1)-nvimf(i-1))*0.5/sqrt(1-nvimf(i)*nvimf(i));
    else
        omgcos(i)=-9999;
    end
end
omgcos(1)=-9999;
omgcos(vlength)=-9999;

jj=1;
for i=1:vlength,
    if omgcos(i)>-1000
        ddd(jj)=i;
        temp(jj)=omgcos(i);
        jj=jj+1;
    end
end
temp2=spline(ddd,temp,dd);
omgcos=temp2;


for i=1:vlength,
    omega(i)=omgcos(i);
end

pi2=pi*2;
omega=omega/dt;

dailiangren

楼主最好可以把图贴出来，那样有经验的人可以一下子就明白你的问题关键所在了。

一图胜于千言啊:lol

dailiangren

在网上搜索了一下这个程序，找到了这个地方。

  HHT MATLAB program

The current MATLAB program of HHT includes the following modules:



1)      eemd.m. The module that performs EMD/EEMD;



2)      ifndq.m. The module that calculates an instantaneous frequency for a given IMF; and



3)      extrema.m. A supporting module that is used by the previous two modules.

　

4)    significance.m This is used to obtain the "percenta" line based on Wu and Huang

      (2004).

　

5)    dist_value.m  This is a utility program being called by "significance.m".

　

The modules of the program are suggested to be in a folder (for example folder “EEMD”) under the MATLAB “toolbox” folder (C:\MATLAB6p5\toolbox\EEMD). By adding the MATLAB path to that directory, the programs are ready to be used. Some interface related information is given in each program. The “help” command can be used to get these information.



The above programs are the preliminary version of an HHT software. In future, more modules that provide more features such statistical test, calculation of marginal spectrum, and visualization in time-frequency-energy domain, will be added.

http://rcada.ncu.edu.tw/research1_clip_program.htm

jjzhuzhu

nstfreq.m中dt代表什么意思啊，是采样频率分之一吗？好像出来频率不太对啊。

baobao1982

原帖由 dailiangren 于 2008-4-18 10:08 发表

                               
登录/注册后可看大图

楼主最好可以把图贴出来，那样有经验的人可以一下子就明白你的问题关键所在了。

一图胜于千言啊:lol

我替楼主贴出来

baobao1982

原帖由 shenlan888 于 2008-4-14 21:46 发表

                               
登录/注册后可看大图

我的程序如下:
clear
fs=200;
t=0:1/fs:1;
s=cos(10*pi*t)+2*cos(40*pi*t);
figure(1)
plot(t,s)
c=emd(s);%进行emd分解
figure(2)
b=size(c,1);
for j=1:1:b
    subplot(b+1,1,j)
    plot(t,c(j,:))
...

估计是你的程序有问题
你可以试试这个
clear
fs=200;
t=0:1/fs:1;
s=cos(10*pi*t)+2*cos(40*pi*t);
figure(1)
plot(t,s)
c=emd(s);%进行emd分解
figure(2)
b=size(c,1);
for j=1:1:b
    subplot(b+1,1,j)
    plot(t,c(j,:))
end
subplot(b+1,1,b+1)
plot(t,sum(c))
% c(1,:)=cos(10*pi*t);
% c(2,:)=2*cos(40*pi*t);
l=length(s);
t=t(2:end-1);
[A,f,tt]=hhspectrum(c);%生成希尔伯特谱
[im,tt]=toimage(A,f);%画出希尔伯特谱
disp_hhs(im);
colormap(flipud(gray)) %设置背景颜色

求得的频率是归一化频率，在乘以200就行了
频谱图见下面

jhhxl

我现在还在学习中，向大家学习哈

zhuxiaoxun

为什么我直接用你的程序运行结果和你的结果不一样呢？得到的图只在0.4HZ有频率，只出现了一个频率。另一个找不到。请指教。谢谢。

cboboc

针对以上程序， omgcos(i)=abs(nvimf(i+1)-nvimf(i-1))*0.5/sqrt(1-nvimf(i)*nvimf(i));这句是什么意思呢，分母中为什么要除以sqrt(1-nvimf(i)*nvimf(i));呢？

devil7663579

原帖由 baobao1982 于 2008-12-4 11:54 发表

                               
登录/注册后可看大图

求得的频率是归一化频率，在乘以200就行了
频谱图见下面

请问：想得到真实信号的时间--频率图，您说乘以200，是在dis_hhts.m文件中乘以200吗？能不能附上部分程序呢？
以下是dis_hhts.m文件中部分程序。您说的是不是在imagesc(t,[0,0.5],im,[inf,0])中乘以fs变成imagesc(t,[0,0.5*fs],im,[inf,0]);就i行了呢？

if fs == 0
  imagesc(t,[0,0.5],im,[inf,0]);
  ylabel('normalized frequency')
else
  imagesc(t,[0,0.5*fs],im,[inf,0]);
% imagesc(t/fs,[0,0.5*fs],im,[inf,0]);%我自己改的可以变实际频率的
  ylabel('frequency')
end

dzkt

imagesc(t/fs,[0,0.5*fs],im,[inf,0]);%我自己改的可以变实际频率的
这个是对的

yrj01

回复 2 # wuqiong_cea 的帖子

请问error = (abs(hilbert(nvimf))-1).^2 一步求的是什么？