求助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分解得不够完美么？
请解答。谢谢
请教搂主，你的时频分布图是用什么程序弄出来的？自己写的？还是用的是instfreq函数？我看你的时频图的频率跟实际很吻合，不像instfreq函数是0-0.5，你是怎么计算的？请指教，谢谢
[ 本帖最后由 zhlong 于 2007-10-25 10:22 编辑 ] rc-hw-0002 ，想请教一下。你是用什么画这个HHtP谱的？我用toimage函数，再用disp_hhs画出的结果不对啊 。很急啊，谢谢了:loveliness:
[ 本帖最后由 zhlong 于 2007-10-25 10:23 编辑 ] 怎么做希尔伯特变换哈？？
我做论文需要把数据进行希尔波变换呢。
%对滤波后的数据进行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. %对滤波后的数据进行hilbert变换
load d:\0920\01-06\f06500.txt
size(f06500)
y=abs(hilbert(f06500));
save d:\0920\01-06\06500hil.txt y -ascii
这个为什么允许有错呢？ 主持
回复 #1 rc-hw-0002 的帖子
看样子你的HHT谱效果也好不到哪去，能不能把你做HHT谱还有边际谱的程序发上来，我们大家都分析一下。回复 #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
回复 #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}希尔伯特边际谱');
用了改进的Hilbert变换后效果好一些
回复 #10 rc-hw-0002 的帖子
我看还是差强人意！为什么在10HZ处还是有点不正常？？你用的EMD是新版的吗？回复 #9 rc-hw-0002 的帖子
你的边际谱和我的如出一辙，看来是一个版本！HHT和我的不太一样。回复 #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,:); rc-hw-0002 ,你能把你改进的hilbert变换发过来吗？