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

sunyuxinhe

程序中，有一小部分是借鉴别人的现成语句。
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，而小波没有。

sunyuxinhe

图的顺序是反的， 看的时候按所标序号看。

yghit08

sunyuxinhe 发表于 2013-3-21 22:07

                               
登录/注册后可看大图

图的顺序是反的， 看的时候按所标序号看。

你至少加个噪声上去啊！你这是考验HHT的分解能力呢还是考验时频工具箱中的瞬时频率的求解呢？？

sunyuxinhe

yghit08 发表于 2013-3-21 22:12

                               
登录/注册后可看大图

你至少加个噪声上去啊！你这是考验HHT的分解能力呢还是考验时频工具箱中的瞬时频率的求解呢？？

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

我前几天看小波去噪方面的资料，你知道SNR为-20左右的信号，怎么去噪吗？或者可以提供点思路吗？

yghit08

sunyuxinhe 发表于 2013-3-21 22:18

                               
登录/注册后可看大图

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

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

不知道。但是看你的结果和原始信号，好像HHT做的更好：基频10Hz，在这上有个波动！而小波就没能做出来这个效果！

bgpz2007

我怎么就这一步错了呢？
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

yghit08

bgpz2007 发表于 2013-3-25 21:24

                               
登录/注册后可看大图

我怎么就这一步错了呢？
clear;
N=1000;

t在你上面的过程中已经使用，换一个量试试？？

bgpz2007

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, :));
而本信号中却出错。
？？？？？？？怎么回事啊？

yghit08

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‘);

bgpz2007

yghit08 发表于 2013-3-25 22:05

                               
登录/注册后可看大图

[A,f,tt]=hhspectrum(imf);

对啊，我就是照你那样该的，不对哦！！！还是出现那个问题！！

yghit08

bgpz2007 发表于 2013-3-25 22:08

                               
登录/注册后可看大图

对啊，我就是照你那样该的，不对哦！！！还是出现那个问题！！

那就不知道了，自己查查看！

yghit08

bgpz2007 发表于 2013-3-25 22:08

                               
登录/注册后可看大图

对啊，我就是照你那样该的，不对哦！！！还是出现那个问题！！

再不行，你就测试一下hhspectrum这个函数，用一正弦信号。接着用两个正弦信号的叠加测试emd这个程序，都没问题的话 ，就是你的问题了！

bgpz2007

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('幅值');

sunyuxinhe

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

bgpz2007

sunyuxinhe 发表于 2013-3-26 21:34

                               
登录/注册后可看大图

我的hhspectrum函数和你的 不完全一样，有时间你可以看看。
%HHSPECTRUM  compute Hilbert-Huang spectr ...

谢谢啊！！！！这个hhs函数到底应该用哪个啊？