声振论坛

 找回密码
 我要加入

QQ登录

只需一步,快速开始

查看: 4807|回复: 25

[HHT] HHT和小波的对比范例-----duffing波

  [复制链接]
发表于 2013-3-21 22:06 | 显示全部楼层 |阅读模式

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

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

x
本帖最后由 yghit08 于 2013-3-25 21:45 编辑

程序中,有一小部分是借鉴别人的现成语句。
clear;
N=1000;
Fs=1000;
t=(1:N)/Fs;
x=cos(2*pi*10*t+0.3*sin(2*pi*20*t));
%----------------------------%
% HHT of Duffing signal%
%----------------------------%
imf=emd(x);
figure(1);
a=size(imf);
subplot(a(1)+1,1,1),plot(x);ylabel('x(t)');
for i=1:a(1)
    subplot(a(1)+1,1,i+1);
    plot(imf(i,:));
    ylabel(['imf',int2str(i)]);
end
[A,f,t]=hhspectrum(imf(1:end, :))
[E,t,Cenf]=toimage(A,f,t,length(t));   
figure;
set(gcf,'Color','w');
mesh(t/N,Cenf*Fs,E);
set(gca,'YDir','normal');
axis([0 1 0 50]);
colormap('jet')
colorbar;
xlabel('时间 t/s');
ylabel('频率 f/Hz');
title('Hilbert-Huang Spectrum') ;

N=length(Cenf);
for k=1:size(E,1)
    bjp(k)=sum(E(k,:))*1/Fs;
end
figure;
plot(Cenf(1,:)*Fs,bjp);
xlabel('频率 f/Hz');
ylabel('幅值');
title('边际谱');
%----------------------------%
% tfrscalo(小波) of Duffing signal%
%----------------------------%
x=x';
x=Hilbert(x);
[tfr,t,f]=tfrscalo(x,(1:N),sqrt(N),0.0001,0.1);
figure;
mesh(t/Fs,f*1000,abs(tfr));
xlabel('时间 t/s');
ylabel('频率 f/Hz');
title('tfrscalo') ;
axis([0 1  0 50]);
Duffing信号的频率应围绕10Hz波动,并出现的波内调频现象,即在信号一个完整
周期内,信号频率发生变化,波内调频频率为20Hz。

HHT检测出了20Hz,而小波没有。

04.png
03.png
02.png
01.png
回复
分享到:

使用道具 举报

 楼主| 发表于 2013-3-21 22:07 | 显示全部楼层
图的顺序是反的, 看的时候按所标序号看。
发表于 2013-3-21 22:12 | 显示全部楼层

你至少加个噪声上去啊!你这是考验HHT的分解能力呢还是考验时频工具箱中的瞬时频率的求解呢??
 楼主| 发表于 2013-3-21 22:18 | 显示全部楼层
yghit08 发表于 2013-3-21 22:12
你至少加个噪声上去啊!你这是考验HHT的分解能力呢还是考验时频工具箱中的瞬时频率的求解呢??

我只是随便弄一个,就比较一下,二者不一样。

我前几天看小波去噪方面的资料,你知道SNR为-20左右的信号,怎么去噪吗?或者可以提供点思路吗?
发表于 2013-3-22 08:52 | 显示全部楼层
sunyuxinhe 发表于 2013-3-21 22:18
我只是随便弄一个,就比较一下,二者不一样。

我前几天看小波去噪方面的资料,你知道SNR为-20左右的信 ...

不知道。但是看你的结果和原始信号,好像HHT做的更好:基频10Hz,在这上有个波动!而小波就没能做出来这个效果!
发表于 2013-3-25 21:24 | 显示全部楼层
本帖最后由 bgpz2007 于 2013-3-25 21:27 编辑

我怎么就这一步错了呢?
clear;
N=1000;
Fs=1000;
t=(1:N)/Fs;
x=cos(2*pi*10*t+0.3*sin(2*pi*20*t));
>> imf=emd(x);
figure(1);
>> a=size(imf);
>> subplot(a(1)+1,1,1),plot(x);ylabel('x(t)');
>> for i=1:a(1)
    subplot(a(1)+1,1,i+1);
    plot(imf(i,:));
    ylabel(['imf',int2str(i)]);
end
>> [A,f,t]=hhspectrum(imf(1:end, :))
??? Output argument "A" (and maybe others) not assigned during call to "D:\matlab 2007\work\pack_emd\package_emd\utils\hhspectrum.m (hhspectrum)".

Error in ==> hhspectrum at 20
if nargin < 2

这是我的hhspectrum.m文件
function [A,f,tt] = hhspectrum(imf,t,l,aff)

% [A,f,tt] = HHSPECTRUM(imf,t,l,aff) computes the Hilbert-Huang spectrum
%
% inputs:
%         - imf : matrix with one IMF per row
%   - t   : time instants
%   - l   : estimation parameter for instfreq
%   - aff : if 1, displays the computation evolution
%
% outputs:
%   - A   : amplitudes of IMF rows
%   - f   : instantaneous frequencies
%   - tt  : truncated time instants
%
% calls:
%   - hilbert  : computes the analytic signal
%   - instfreq : computes the instantaneous frequency

if nargin < 2

  t=1:size(imf,2);

end

if nargin < 3

  l=1;

end

if nargin < 4

  aff = 0;

end

lt=length(t);

tt=t((l+1):(lt-l));

for i=1:(size(imf,1)-1) % This is the original statement
%     for i=1:(size(imf,1)) % I modify this only to veryfy if the HHT spectrum is caculated from the imfs exclucding the residual

  an(i,:)=hilbert(imf(i,:)')';
  f(i,:)=instfreq(an(i,:)',tt,l)';
  A=abs(an(:,l+1:end-l));

  if aff

    disp(['mode ',int2str(i),' trait?'])

  end

end

发表于 2013-3-25 21:29 | 显示全部楼层
bgpz2007 发表于 2013-3-25 21:24
我怎么就这一步错了呢?
clear;
N=1000;

t在你上面的过程中已经使用,换一个量试试??
发表于 2013-3-25 22:02 | 显示全部楼层
yghit08 发表于 2013-3-25 21:29
t在你上面的过程中已经使用,换一个量试试??

不对吧,程序的第一行就错了
[A,f,tt]=hhspectrum(imf);
??? Output argument "A" (and maybe others) not assigned during call to "D:\matlab 2007\work\pack_emd\package_emd\utils\hhspectrum.m (hhspectrum)".

Error in ==> hhspectrum at 20
if nargin < 2

有些别的信号可以用
[A,f,tt]=hhspectrum(imf(1:end, :));
而本信号中却出错。
???????怎么回事啊?
发表于 2013-3-25 22:05 | 显示全部楼层
本帖最后由 yghit08 于 2013-3-25 22:08 编辑
bgpz2007 发表于 2013-3-25 22:02
不对吧,程序的第一行就错了
[A,f,tt]=hhspectrum(imf);
??? Output argument "A" (and maybe others)  ...
在第一行加上clear,clc再试试
[A,f,tt]=hhspectrum(imf);不行的话
[A,f,tt]=hhspectrum(imf‘);
发表于 2013-3-25 22:08 | 显示全部楼层
yghit08 发表于 2013-3-25 22:05
[A,f,tt]=hhspectrum(imf);

对啊,我就是照你那样该的,不对哦!!!还是出现那个问题!!
发表于 2013-3-25 22:10 | 显示全部楼层
bgpz2007 发表于 2013-3-25 22:08
对啊,我就是照你那样该的,不对哦!!!还是出现那个问题!!

那就不知道了,自己查查看!
发表于 2013-3-25 22:11 | 显示全部楼层
bgpz2007 发表于 2013-3-25 22:08
对啊,我就是照你那样该的,不对哦!!!还是出现那个问题!!

再不行,你就测试一下hhspectrum这个函数,用一正弦信号。接着用两个正弦信号的叠加测试emd这个程序,都没问题的话 ,就是你的问题了!
发表于 2013-3-26 08:03 | 显示全部楼层
yghit08 发表于 2013-3-25 22:11
再不行,你就测试一下hhspectrum这个函数,用一正弦信号。接着用两个正弦信号的叠加测试emd这个程序,都没 ...

余弦函数没问题啊!!!我运行了下!!!为什么楼主的就会出错? 那我的问题是在哪里?
clear;
N=1000;
n=1:N;
fs=1000;
t=n/fs;
fx=10;
fy=50;
x=cos(2*pi*fx*t);
y=10*cos(2*pi*fy*t);
z=x+y;
data=z;
imf=emd(data);                        %对输入信号进行EMD分解   
[A,f,t]=hhspectrum(imf);            %对IMF分量求取瞬时频率与振幅:A:是每个IMF的振幅向量,f:每个IMF对应的瞬时频率,t:时间序列号
[E,t,Cenf]=toimage(A,f);            %将每个IMF信号合成求取Hilbert谱,E:对应的振幅值,Cenf:每个网格对应的中心频率  这里横轴为时间,纵轴为频率        
                                                   %即时频图(用颜色表示第三维值的大小)和三维图(三维坐标系:时间,中心频率,振幅)         
cemd_visu(data,1:length(data),imf);   %显示每个IMF分量及残余信号--------------------------------------------
disp_hhs(E);                          %希尔伯特谱----------------------------------------------------------
%画出边际谱
%N=length(Cenf);%设置频率点数   %完全从理论公式出发。网格化后中心频率很重要,大家从连续数据变为离散的角度去思考,相信应该很容易理解
for k=1:size(E,1)
    bjp(k)=sum(E(k,:))*1/fs;
end
figure(3);
plot(Cenf(1,:)*fs,bjp);  % 作边际谱图   进行求取Hilbert谱时频率已经被抽样成具有一定窗长的离散频率,所以此时的频率轴已经是中心频率
xlabel('频率 / Hz');
ylabel('幅值');
 楼主| 发表于 2013-3-26 21:34 | 显示全部楼层
bgpz2007 发表于 2013-3-26 08:03
余弦函数没问题啊!!!我运行了下!!!为什么楼主的就会出错? 那我的问题是在哪里?
clear;
N=1000; ...

我的hhspectrum函数和你的 不完全一样,有时间你可以看看。
%HHSPECTRUM  compute Hilbert-Huang spectrum
%
% [A,f,tt] = HHSPECTRUM(x,t,l,aff) computes the Hilbert-Huang spectrum
%
% inputs:
%   - x   : matrix with one signal per row
%   - t   : time instants
%   - l   : estimation parameter for instfreq (integer >=1 (1:default))
%   - aff : if 1, displays the computation evolution
%
% outputs:
%   - A   : instantaneous amplitudes
%   - f   : instantaneous frequencies
%   - tt  : truncated time instants
%
% calls:
%   - hilbert  : computes the analytic signal
%   - instfreq : computes the instantaneous frequency
%   - disprog : displays the computation evolution
%
%Examples:
%
%s = randn(1,512);
%imf = emd(s);
%[A,f,tt] = hhspectrum(imf(1:end-1,:));
%
%s = randn(10,512);
%[A,f,tt] = hhspectrum(s,1:512,2,1);
%
%
% See also
%  emd, toimage, disp_hhs
%
% G. Rilling, last modification 3.2007
% gabriel.rilling@ens-lyon.fr

function [A,f,tt] = hhspectrum(x,t,l,aff)

error(nargchk(1,4,nargin));

if nargin < 2

  t=1:size(x,2);

end

if nargin < 3

  l=1;

end

if nargin < 4

  aff = 0;

end

if min(size(x)) == 1
        if size(x,2) == 1
                x = x';
                if nargin < 2
                        t = 1:size(x,2);
                end
        end
        Nmodes = 1;
else
        Nmodes = size(x,1);
end

lt=length(t);

tt=t((l+1):(lt-l));% this reduce the length of t from 1000 to 998

for i=1:Nmodes

  an(i,:)=hilbert(x(i,:)')';
  f(i,:)=instfreq(an(i,:)',tt,l)';
  A=abs(an(:,l+1:end-l));

  if aff
        disprog(i,Nmodes,max(Nmodes,100))
  end

end
发表于 2013-3-27 10:04 | 显示全部楼层
sunyuxinhe 发表于 2013-3-26 21:34
我的hhspectrum函数和你的 不完全一样,有时间你可以看看。
%HHSPECTRUM  compute Hilbert-Huang spectr ...

谢谢啊!!!!这个hhs函数到底应该用哪个啊?
您需要登录后才可以回帖 登录 | 我要加入

本版积分规则

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

GMT+8, 2024-11-25 03:07 , Processed in 0.070510 second(s), 21 queries , Gzip On.

Powered by Discuz! X3.4

Copyright © 2001-2021, Tencent Cloud.

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