声振论坛

 找回密码
 我要加入

QQ登录

只需一步,快速开始

查看: 3820|回复: 12

[HHT] 求出HHT频率为何不匹配?

[复制链接]
发表于 2008-4-14 21:46 | 显示全部楼层 |阅读模式

马上注册,结交更多好友,享用更多功能,让你轻松玩转社区。

您需要 登录 才可以下载或查看,没有账号?我要加入

x
我的程序如下:
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);
这个频率是正确的
回复
分享到:

使用道具 举报

发表于 2008-4-17 15:54 | 显示全部楼层
黄给出了一个求解瞬时频率的程序,而他本人也关于瞬时频率问题写了专门的文章。下面给出他的程序代码:他的求解方法与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;
发表于 2008-4-18 10:08 | 显示全部楼层
楼主最好可以把图贴出来,那样有经验的人可以一下子就明白你的问题关键所在了。

一图胜于千言啊:lol
发表于 2008-4-18 10:18 | 显示全部楼层

回复 2楼 的帖子

在网上搜索了一下这个程序,找到了这个地方。
  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
发表于 2008-4-20 16:18 | 显示全部楼层
nstfreq.m中dt代表什么意思啊,是采样频率分之一吗?好像出来频率不太对啊。
发表于 2008-12-4 11:19 | 显示全部楼层



我替楼主贴出来

1.fig

9.08 KB, 下载次数: 90

发表于 2008-12-4 11:54 | 显示全部楼层
原帖由 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就行了
频谱图见下面

3.fig

12.05 KB, 下载次数: 84

发表于 2009-1-10 10:14 | 显示全部楼层
我现在还在学习中,向大家学习哈
发表于 2009-1-21 10:26 | 显示全部楼层

回复 7楼 baobao1982 的帖子

为什么我直接用你的程序运行结果和你的结果不一样呢?得到的图只在0.4HZ有频率,只出现了一个频率。另一个找不到。请指教。谢谢。
发表于 2010-3-10 10:08 | 显示全部楼层
针对以上程序, omgcos(i)=abs(nvimf(i+1)-nvimf(i-1))*0.5/sqrt(1-nvimf(i)*nvimf(i));这句是什么意思呢,分母中为什么要除以sqrt(1-nvimf(i)*nvimf(i));呢?
发表于 2010-5-8 07:52 | 显示全部楼层
原帖由 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
发表于 2010-5-8 15:13 | 显示全部楼层
imagesc(t/fs,[0,0.5*fs],im,[inf,0]);%我自己改的可以变实际频率的
这个是对的
发表于 2011-5-13 08:04 | 显示全部楼层
回复 2 # wuqiong_cea 的帖子

请问error = (abs(hilbert(nvimf))-1).^2 一步求的是什么?
您需要登录后才可以回帖 登录 | 我要加入

本版积分规则

QQ|小黑屋|Archiver|手机版|联系我们|声振论坛

GMT+8, 2024-11-25 20:37 , Processed in 0.073139 second(s), 22 queries , Gzip On.

Powered by Discuz! X3.4

Copyright © 2001-2021, Tencent Cloud.

快速回复 返回顶部 返回列表