Hilbert变换及其边际谱求解代码

%%查找最大值
for i=2:point*N
if (x(i)>x(i-1)) &&(x(i)>x(i+1))
j=j+1;
max(j)=x(i);
max_num(j)=i/point/f1;
point=32; %每周期采样点数
fmax=110;
CN=3; %模态数量
cn=0;
M=1; %求每个模态分量时的迭代次数
m=0;
f0=50; %信号1的频率
f2=1; %信号1的频率
g0=0;
xlabel('t/s');
ylabel('f/Hz');
zlabel('a/V');
if i>1
diffANG(i)=ANG(i)-ANG(i-1);
f(cn,i)=(ANG(i)-ANG(i-1))/(1/f0/point)/(2*pi);
if(f(cn,i)>fmax)|(f(cn,i)<(0-fmax))
min_num(j)=i/point/f1;
end
end
%最小值三次样条插值
t=0: 1/f0/point: 1/f0*N;
min_spline=spline(min_num,min,t);
t=0: 1/f0/point: 1/f0*N;
for i=1:length(t);
temp(i)= a(1,i);
end
figure(3)
plot(t,temp)
title('谐波分量的幅值-时间图')
%%%%%%%%%%%%%%%%%%%%%%%%
figure(2)
t=0: 1/f1/point: 1/f1*N;
for i=1:length(t);
for j=1:length(fview);
aview(j,i)=0;
aview2(j)=0;
end
end
for i=1:length(t);
for j=1:length(fview);
for k=1:CN;
title('第1个IMF的相角');
hold on
plot(t,diffANG,'g')
hold off
end
end
t=0: 1/f0/point: 1/f0*N;
fview=0:1:fmax;
for j=1:CN
for i=1:length(t);
temp(i)=C(j,i);
end
subplot(CN,1,j);
plot(t,temp)
end
figure(21)
[xx,yy]=meshgrid(t,fview);
surf(xx,yy,aview);
t=0: 1/f0/point: 1/f0*N;
plot(t,ADdata)
for cn=1:CN
for m=1:M
j=0;
clear max;
clear max_num;
clear min;
clear min_num;
%%求插值后最大值、最小值包络线的平均值
i=0;
for t=0: 1/f0/point: 1/f0*N;
i=i+1;
m1_array(i)=(max_spline(i)+min_spline(i))/2;
%查找最小值
j=0;
for i=2:point*N
if (x(i)<x(i-1)) &&(x(i)<x(i+1))
j=j+1;
min(j)=x(i);
t=0: 1/f0/point: 1/f0*N;
for i=1:length(t);
temp(i)=f(1,i);
end
figure(4)
plot(t,temp)
title('频率-时间图')
%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%省去了判断：h1_array连续几次的极值点个数与过零点数相同或者相差1
i=0;
for t=0: 1/f0/point: 1/f0*N;
i=i+1;
x(i)=h1_array(i); %%H1的查找多做几次
f(cn,i)=f(cn,i-1);
end
end
end
if cn==1;
t=0: 1/f0/point: 1/f0*N;
figure(24)
plot(t,ANG,'r')
i=i+1;
a(cn,i)=abs(hil(i));
end
%%%计算相角
i=0;
f(1)=0;
for t=0: 1/f0/point: 1/f0*N;
i=i+1;
ANG(i)=atan(imag(hil(i))/real(hil(i)));
end
end
%最大值三次样条插值
t=0: 1/f0/point: 1/f0*N;
max_spline=spline(max_num,max,t); %？？？？两边的两点不能加入运算，导致两边的两点插值误差非常大
ADdata(i)=x(i);
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
hil=hilbert(h1_array);
%%%%计算幅值
i=0;
for t=0: 1/f0/point: 1/f0*N;
if (abs(fview(j)-f(k,i))) <
end
end
end
end
%========边缘谱==========
for j=1:length(fview);
for i=2*point:length(t)-2*point;
aview2(j)=aview2(j)+aview(j,i);
end
aview2(j)=aview2(j)/(length(t)-4*point);
end
%%%%%%%%%%%%%%%%%%%%%%%%
end
%%求H1
i=0;
for t=0: 1/f0/point: 1/f0*N;
i=i+1;
h1_array(i)=x(i)-m1_array(i);
end
%===============================
%======Hilbert边缘谱============
%调试记录：
%时变系统的HHT:f时变的HHT，含50Hz,51Hz
%===============================
clc;
clear;
N=100; %采样周期
end
ADdata(i)=220*sin(2*pi*f1*t+g0); %+sin(2*pi*f2*t); %+0.1*rand();
ADdata0(i)=ADdata(i);
x(i)=ADdata(i);
end
figure(1)
end
end
i=0;
for t=0: 1/f0/point: 1/f0*N;
i=i+1;
C(cn,i)=h1_array(i);
x(i)=ADdata(i)-h1_array(i); %%H1的查找多做几次
t=0: 1/f0/point: 1/f0*N;
%AD采样数据
i=0;
f1=50;
g0=0;
for t=0: 1/f0/point: 1/f0*N
i=i+1;
if t > 1/f0*N/2;
f1=51;
g0=0;%34.56/180*pi;