l利用希尔伯特直接正交发求瞬时频率

一蓑烟雨

%l利用希尔伯特直接正交发求瞬时频率
clear ;clc;
t=0:300;
x0=@(t)sin(0.027*pi*t)-sin(2*(0.027*pi)*t+pi/2)+sin(4*(0.027*pi)*t);%信号x(t)
y0=x0(t);
figure(1);
plot(t,y0);
axis([0 300 -4 4]);
x1=@(t)(((t-120)==0)|((t-240)==0));%脉冲信号v(t)
y1=x1(t);
figure(2);
plot(t,y1);
axis([0 300 -0.1 0.1]);
y2=y0+y1;
imf = emd(y2);%EMD 分解
x3=imf(1,:);%用第一个imf
x4=x3';
figure(3);
subplot(411);plot(t,imf(1,:));
axis([0 300 -1 1]);
subplot(412);plot(t,imf(2,:));
axis([0 300 -1 1]);
subplot(413);plot(t,imf(3,:));
axis([0 300 -1 1]);
subplot(414);plot(t,imf(4,:));
axis([0 300 -1 1]);
[instf,t]=instfreq(x4);%用HT求得瞬时频率
instf = sym('instf');
y=int(instf);%积分
b=cos(y);%求的F(T)
z1=1-b^2;
z2=sqrt(z1);%求得Y(T)
z3=b./z2;
z4=atan(z3);
dt=diff(t);
dx=diff(z4);
f=0.5.*pi*dx./dt;%求瞬时频率
f1=sym(f);
figure(4);
plot(t,f1);
大家帮我看看程序错在哪里了?尤其是最后图显示不出来，

shuihai707

使用直接正交法之前要对IMF进行调幅调频分解，直接对IMF进行DQ法效果很差。

shuihai707

三维图我还没画过，我现在做AM-FM分解，论坛里有一个AM-FM程序，但有点问题，我正修改呢。

lzqfhz

本人新手，没搞得太明白。不过有几点想说的是
1、看你的思路好像是求的瞬时频率，[instf,t]=instfreq(x4）；再瞬时频率积分y=int(instf)求相位；然后求归一化调频信号b=cos(y);之后求正交分量(sin(y))z1=1-b^2;z2=sqrt(z1);再求相位z4=atan(z3);最后求频率dx=diff(z4);这想法很新颖。但是这不是Huang所说的直接正交求分量的方法。因为huang给出的方法就是为了避开hilbert,instfreq函数输入一般要求是解析信号的，而解析信号一般还是通过Hilbert求解的。还有你这里z2=sqrt(z1),结果永远为正。所以相位z4=atan(z3);肯定有问题。
2、instf = sym('instf'); y=int(instf); 运行完后instf就成了一个字符串了，而从你后面的思路来看可能还想让它是个瞬时频率的值构成的向量（不知可否这样理解）。所以后面的表达式里即使有instf也不是瞬时频率了，它仅仅是一个串。
3、Huang的文章里给的思路应该是，拿到一个IMF，直接求包络进行归一化，从而得到相应的调频分量，然后再求正交分量，最后求瞬时相位和瞬时频率。
拙见！