包络线和包络谱计算程序

%一个求包络线和包络谱的程序

%现代希尔伯特变换解调分析：

%带通滤波；希尔伯特变换获得信号时域的包络线；用fft变换获得包络谱

%如何获得包络线？

%信号经希尔伯特变换不能直接得到包络，设信号x的希尔伯特变换为y，则平方和

%x.^2+y.^2(或者再开根号，直接取平方和的效果为好)才是信号x的包络。

%构造实验数据

clear all;close all;

t=0:0.005:1*pi;

fs=10000;

s=4*sin(2*200*pi*t).*(sin(2*4500*pi*t))+25*(sin(2*4500*pi*t));

figure(1);

subplot(211);plot(t,s);title('原始信号');

%运用小波方法滤波

[c,l]=wavedec(s,1,'db10');

d1=wrcoef('d',c,l,'db10');

a1=0;

subplot(212);plot(d1);title('滤波后重构的高频信号');

%希尔伯特变换求包络线

y=hilbert(d1);

y1=abs(d1+y*j); %这是取得包络线的三种方程。看一看哪种效果好。

%y1=abs(y); %或者z=x.^2+y.^2; 有的取得是abs(y),但是不推荐用。%y1=d1.^2+y.^2; %通过分析，该方程在包络谱中的效果最好，即取二者平方和。

figure(2);

subplot(211);

hold on

plot(t,s);

plot(t,y1,'r');title('包络线');

hold off

%FFT求包络谱

N=1024;

p=abs(fft(y1,N));

subplot(212);

plot((0:N/2-1)/N*fs,p(1:N/2));%只需取到半频，即fs/2

%f=(0:N-1)*fs/N; plot(f,p); %横坐标是在fs上，其中以fs/2为轴中心对称。

title('包络谱');xlabel('频率');ylabel('功率谱');

%对比信号直接的傅里叶变换功率谱与包络谱

p1=abs(fft(s,N));

figure(3)

hold on

plot((0:N/2-1)/N*fs,p(1:N/2),'g'); grid;%包络谱

plot((0:N/2-1)/N*fs,p1(1:N/2),'r'); %傅里叶变换功率谱hold off