求助EMD分解误差和Hilbert变换

rc-hw-0002

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;
%在后面添加一个趋势向可以把三个分量完全分解出来

EMD分解结果

得时频分布：

时频分布

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

边际谱

fay1014

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

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

linqin1201

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

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

suj568

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

%对滤波后的数据进行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

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

这个为什么允许有错呢？

liliang

主持

zhangnan3509

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

rc-hw-0002

在别人的基础上改了一下：
function [S,freq] = 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

画边际谱
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

zhangnan3509

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

zhangnan3509

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

zhangnan3509

这是我的HHT谱
function h1= nspab(data,nyy,minw,maxw,dt)
%----- Get dimensions (number of time points and components)
[npt,knb] = 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
[nx,ny]=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

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

rc-hw-0002

不是新版本，你分解的结果如何？发个图上来看一下吧