请各位大虾过来瞅瞅哦：HHT变换中的 瞬时能量谱

civil1314

根据破凰的边际谱程序，改了一下，做瞬间能量的图。
但是由于本人信号处理不熟悉，无法判断结果正确与否。
先把程序和图贴上来，大家帮忙看看啊。
%Hilbert瞬时能量谱图
function [N,IE]=ie(x,n,dt,fs)
%x 原始信号
%n采样点数
%dt采样周期
%fs=1/dt;%采样频率
imf=emd(x);
%[S,freq]=hhte(imf,2) %Y.Lzhang
[A,fa,tt]=hhspectrum(imf);
[E,tt1]=toimage(A,fa,tt,length(tt));
for k=1:size(E,2)
    IE(k)=sum(E(k,:))^2*2*pi*fs;  %瞬时能量谱
end
N=(0:n-3)/n

clear all
n=1000;
dt=0.002;
%n=250;
%dt=0.002;
fs=1/dt;%采样频率
for i=1:n
    t(i)=i*dt;
x(i)=exp(-0.5*t(i))*(sin(10*pi*t(i))+t(i))+0.2*cos(2*pi*t(i));
%x(i)=cos(2*pi*t(i)+0.5*sin(2*pi*t(i)))*exp(-0.2*t(i))+0.05*sin(30*i*t(i));
%%任春论文
end
[N,IE]=ie(x,1000,0.002,1/0.002);
%plot(N(1:25),IE(1:25),'k');
plot(N(1:50),IE(1:50),'k');
xlabel('时间/ Hz','Fontsize',10);
set(gca,'XTick',0:0.01:0.05)
ylabel('瞬时能量','Fontsize',10);
title('时间','FontSize',10)

真人

??? Error: File: httshun.m Line: 9 Column: 1
Function definitions are not permitted at the prompt or in scripts.
%Hilbert瞬时能量谱图
function [N,IE]=ie(x,n,dt,fs)

=============
么意思啊？

civil1314

IE是自定义的一个函数啊，我贴的程序里面有这个啊。
我用的emd是法国的那个程序啊，网站里面有下载的地址的。

blusun_2008

想问问楼主：
for k=1:size(E,2)
    IE(k)=sum(E(k,:))^2*2*pi*fs;  %瞬时能量谱

为什么不是size(E,1)？E是两维的，取它的谱值是用size(E,1)还是用size(E,2)？

civil1314

不好意思，这一块我也不是很明白。

haifeng1213

用1是为了求固有模态函数的个数把
新手小见解

jhhxl

我运行这个程序也出现了同样的问题，请问你是怎么解决的呢

qqvirile

期待问题的解决

xiaokang

能量谱就是H的平方然后在对频率积分，那IE(k)=sum(E(k,:))^2*2*pi*fs;  %瞬时能量谱这句怎么理解呢？E(k,:)是先平方再求和还是先求和在平方？2*pi*fs是什么含义啊？？最终求出来横坐标是什么含义呢？N=(0:n-3)/n为什么这么取？

nkdtxf

期待！没人解答呢

starshine1985

回复 civil1314 的帖子

不对吧！
1、破凰边际谱是 bjp(k)=sum(E(k,:))*1/fs; 1//fs为时间分辨率
所以，瞬时能量谱应该为IE(k)=sum(E(k,:).^2)*2*pi*fs/(2*N);fs/(2*N)为频率分辨率
总之，都相当用矩形面积代替原来的积分式！
欢迎各位大牛拍砖！！！

凤鸣荆楚

回复 11 # starshine1985 的帖子

按照变成定积分的定义来求的话，为什么求和不是整体求而是求和后相当于乘以dt和dw啊，而且这里面2pi是怎么出来的啊，

凤鸣荆楚

希望各位大牛解答一下啊，不应该求能量是出现2pi啊

fred110

到底是谁对啊？是楼主还是1985呢

王若尧

for k=1:size(E,2)
IE(k)=sum(E(:,k).^2)*2*pi*fs/(2*N)
end
这样子对吗