HHT和小波的对比范例-----duffing波
本帖最后由 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
=hhspectrum(imf(1:end, :))
=toimage(A,f,t,length(t));
figure;
set(gcf,'Color','w');
mesh(t/N,Cenf*Fs,E);
set(gca,'YDir','normal');
axis();
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);
=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();
Duffing信号的频率应围绕10Hz波动，并出现的波内调频现象，即在信号一个完整
周期内，信号频率发生变化，波内调频频率为20Hz。
HHT检测出了20Hz，而小波没有。
图的顺序是反的， 看的时候按所标序号看。 sunyuxinhe 发表于 2013-3-21 22:07 static/image/common/back.gif
图的顺序是反的， 看的时候按所标序号看。
你至少加个噪声上去啊！你这是考验HHT的分解能力呢还是考验时频工具箱中的瞬时频率的求解呢？？ yghit08 发表于 2013-3-21 22:12 static/image/common/back.gif
你至少加个噪声上去啊！你这是考验HHT的分解能力呢还是考验时频工具箱中的瞬时频率的求解呢？？
我只是随便弄一个，就比较一下，二者不一样。
我前几天看小波去噪方面的资料，你知道SNR为-20左右的信号，怎么去噪吗？或者可以提供点思路吗？ sunyuxinhe 发表于 2013-3-21 22:18 static/image/common/back.gif
我只是随便弄一个，就比较一下，二者不一样。
我前几天看小波去噪方面的资料，你知道SNR为-20左右的信 ...
不知道。但是看你的结果和原始信号，好像HHT做的更好：基频10Hz，在这上有个波动！而小波就没能做出来这个效果！ 本帖最后由 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
>> =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 = hhspectrum(imf,t,l,aff)
% = 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
bgpz2007 发表于 2013-3-25 21:24 static/image/common/back.gif
我怎么就这一步错了呢？
clear;
N=1000;
t在你上面的过程中已经使用，换一个量试试？？ yghit08 发表于 2013-3-25 21:29 static/image/common/back.gif
t在你上面的过程中已经使用，换一个量试试？？
不对吧，程序的第一行就错了
=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
有些别的信号可以用
=hhspectrum(imf(1:end, :));
而本信号中却出错。
？？？？？？？怎么回事啊？ 本帖最后由 yghit08 于 2013-3-25 22:08 编辑
bgpz2007 发表于 2013-3-25 22:02 static/image/common/back.gif
不对吧，程序的第一行就错了
=hhspectrum(imf);
??? Output argument "A" (and maybe others)...在第一行加上clear,clc再试试
=hhspectrum(imf);不行的话
=hhspectrum(imf‘);
yghit08 发表于 2013-3-25 22:05 static/image/common/back.gif
=hhspectrum(imf);
对啊，我就是照你那样该的，不对哦！！！还是出现那个问题！！ bgpz2007 发表于 2013-3-25 22:08 static/image/common/back.gif
对啊，我就是照你那样该的，不对哦！！！还是出现那个问题！！
那就不知道了，自己查查看！ bgpz2007 发表于 2013-3-25 22:08 static/image/common/back.gif
对啊，我就是照你那样该的，不对哦！！！还是出现那个问题！！
再不行，你就测试一下hhspectrum这个函数，用一正弦信号。接着用两个正弦信号的叠加测试emd这个程序，都没问题的话 ，就是你的问题了！ yghit08 发表于 2013-3-25 22:11 static/image/common/back.gif
再不行，你就测试一下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分解
=hhspectrum(imf); %对IMF分量求取瞬时频率与振幅：A：是每个IMF的振幅向量,f:每个IMF对应的瞬时频率，t：时间序列号
=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('幅值'); bgpz2007 发表于 2013-3-26 08:03 static/image/common/back.gif
余弦函数没问题啊！！！我运行了下！！！为什么楼主的就会出错？ 那我的问题是在哪里？
clear;
N=1000; ...
我的hhspectrum函数和你的 不完全一样，有时间你可以看看。
%HHSPECTRUMcompute Hilbert-Huang spectrum
%
% = 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);
% = hhspectrum(imf(1:end-1,:));
%
%s = randn(10,512);
% = hhspectrum(s,1:512,2,1);
%
%
% See also
%emd, toimage, disp_hhs
%
% G. Rilling, last modification 3.2007
% gabriel.rilling@ens-lyon.fr
function = 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
sunyuxinhe 发表于 2013-3-26 21:34 static/image/common/back.gif
我的hhspectrum函数和你的 不完全一样，有时间你可以看看。
%HHSPECTRUMcompute Hilbert-Huang spectr ...
谢谢啊！！！！这个hhs函数到底应该用哪个啊？
页:
[1]
2