rc-hw-0002 发表于 2007-4-20 20:02

求助EMD分解误差和Hilbert变换

t = linspace(1,2,1000);
x1 =cos(2*pi*5*t);
x2 =0.5*cos(2*pi*10*t);
x3 =0.7*cos(2*pi*25*t);
x = x1+x2+x3;
x=x+t;
%在后面添加一个趋势向可以把三个分量完全分解出来

得时频分布:

但是感觉时频分布端点处偏差比较严重,是因为Hibert变换的端点效应么?
还有10HZ处波动比较严重,是因为IMF分解得不够完美么?
请解答。谢谢

fay1014 发表于 2007-4-21 08:55

请教搂主,你的时频分布图是用什么程序弄出来的?自己写的?还是用的是instfreq函数?我看你的时频图的频率跟实际很吻合,不像instfreq函数是0-0.5,你是怎么计算的?请指教,谢谢

[ 本帖最后由 zhlong 于 2007-10-25 10:22 编辑 ]

linqin1201 发表于 2007-4-21 10:20

rc-hw-0002    ,想请教一下。你是用什么画这个HHtP谱的?我用toimage函数,再用disp_hhs画出的结果不对啊 。很急啊,谢谢了:loveliness:

[ 本帖最后由 zhlong 于 2007-10-25 10:23 编辑 ]

suj568 发表于 2007-4-21 20:25

怎么做希尔伯特变换哈??
我做论文需要把数据进行希尔波变换呢。

%对滤波后的数据进行hilbert变换
A = load('d:\shuju\hgt.2000stdf.txt');
    if ~isstruct(A)
      error('The mat file should contain a structure array!');
    end
    cFieldNames = fieldnames(A);
    B = A.(cFieldNames{1});
size(hgt.2000stdf)
y=abs(imag(hilbert(hgt.2000stdf)));
save d:\shuju\hgt.2000stdhil.txt y -ascii

之前弄的总是得到这样的错误提示:
??? Error: File: h.m Line: 8 Column: 9
Unexpected MATLAB expression.

suj568 发表于 2007-4-21 20:26

%对滤波后的数据进行hilbert变换
load d:\0920\01-06\f06500.txt
size(f06500)
y=abs(hilbert(f06500));
save d:\0920\01-06\06500hil.txt y -ascii

这个为什么允许有错呢?

liliang 发表于 2007-4-21 20:44

主持

zhangnan3509 发表于 2007-4-21 22:15

回复 #1 rc-hw-0002 的帖子

看样子你的HHT谱效果也好不到哪去,能不能把你做HHT谱还有边际谱的程序发上来,我们大家都分析一下。

rc-hw-0002 发表于 2007-4-21 23:01

回复 #7 zhangnan3509 的帖子

在别人的基础上改了一下:
function = HHTspe(imf,N,fs);
%%N--------分辨率
%fs-------采样频率
%s--------时间-频率-幅值矩阵
if(nargin<3)
    fs=1;
end
L = size(imf,1);
S = [];
   
clear x z m p freq
   
x = imf';
z = hilbtmf(x); % 用改进的Hilbert变换,效果好一些。
m = abs(z);
p = angle(z);

for i = 1:L
      freq(:,i) = instfreq(z(:,i))*fs; %乘以了采样频率
                  ceilfreq(:,i) = ceil(freq(:,i)*N);
            for j = 1:length(x)-2
            S(ceilfreq(j,i),j+1) = m(j+1,i);
         end
   end
%plot S
figure;
t=(1:length(x))/fs;
f=(1:size(S,1))/N; %
imagesc(t,f,S);
colorbar;
set(gca,'YDir','normal');
xlabel('Time t/s');
ylabel(' Frequency/HZ');

return

rc-hw-0002 发表于 2007-4-21 23:05

回复 #7 zhangnan3509 的帖子

画边际谱
function ms=HHTms(x,N)
% Input-
%      x      - 2-D matrix x(k,n) of the HHT spectrum
%      N       - 分辨率
% Output-
%      ms      - vector ms(k) that specifies the marginal spectrum      surf(h)   shading interp
if(nargin<2)
    N=1;
end
n=size(x);
k=n(1);
n=n(2);
ms=sum(x')'/n;

w=1:length(ms);
w=w/N;
%画边际谱:
figure
plot(w,ms);
xlabel('\fontsize{18}频率 / Hz');
ylabel('\fontsize{18}幅值');
legend('\fontsize{18}希尔伯特边际谱');

rc-hw-0002 发表于 2007-4-21 23:07

用了改进的Hilbert变换后效果好一些

zhangnan3509 发表于 2007-4-21 23:13

回复 #10 rc-hw-0002 的帖子

我看还是差强人意!为什么在10HZ处还是有点不正常??你用的EMD是新版的吗?

zhangnan3509 发表于 2007-4-21 23:16

回复 #9 rc-hw-0002 的帖子

你的边际谱和我的如出一辙,看来是一个版本!HHT和我的不太一样。

zhangnan3509 发表于 2007-4-21 23:20

回复 #8 rc-hw-0002 的帖子

这是我的HHT谱
function h1= nspab(data,nyy,minw,maxw,dt)
%----- Get dimensions (number of time points and components)
= size(data);

%----- Get time interval


%----- Apply Hilbert Transform
data=hilbert(data);
a=abs(data);
omg=abs(diff(unwrap(angle(data))))/(2*pi*dt);

%----- Smooth amplitude and frequency   
filtr=fir1(8,.1);
for i=1:knb
    a(:,i)=filtfilt(filtr,1,a(:,i));
   omg(:,i)=filtfilt(filtr,1,omg(:,i));
end
for i=1:knb
   a(:,i)=filtfilt(filtr,1,a(:,i));
omg(:,i)=filtfilt(filtr,1,omg(:,i));
end

%----- Limit frequency and amplitude
for i=1:knb
   for i1=1:npt-1
      if omg(i1,i) >=maxw,
         omg(i1,i)=maxw;
       a(i1,i)=0;
      elseif omg(i1,i)<=minw,
         omg(i1,i)=minw;
       a(i1,i)=0;
      else
      end
   end
end
clear filtr data

%va=var(omg(200:1200))

%----- Get local frequency
dw=maxw - minw;
wmx=maxw;
wmn=minw;

%----- Construct the ploting matrix
clear p;
h1=zeros(npt-1,nyy+1);
p=round(nyy*(omg-wmn)/dw)+1;
for j1=1:npt-1
   for i1=1:knb
      ii1=p(j1,i1);       
      h1(j1,ii1)=h1(j1,ii1)+a(j1,i1);
   end
end

%----- Do 3-point to 1-point averaging
=size(h1);
%n1=fix(nx/3);
%h=zeros(n1,ny);
%for i1=1:n1
   %h(i1,:)=(h1(3*i1,:)+h1(3*i1-1,:)+h1(3*i1-2,:));
%end
%clear h1;

%----- Do 3-points smoothing in x-direction
fltr=1./3*ones(3,1);
for j1=1:ny
h1(:,j1)=filtfilt(fltr,1,h1(:,j1));
end
clear fltr;


%----- Define the results
%w=linspace(wmn,wmx,ny-1)';
%xs=linspace(t0,t1,nx)';
h1=flipud(rot90(h1));
h1=h1(1:ny-1,:);

tangaoming 发表于 2007-4-22 09:13

rc-hw-0002 ,你能把你改进的hilbert变换发过来吗?

rc-hw-0002 发表于 2007-4-22 10:43

回复 #11 zhangnan3509 的帖子

不是新版本,你分解的结果如何?发个图上来看一下吧
页: [1] 2 3
查看完整版本: 求助EMD分解误差和Hilbert变换