包络线和包络谱计算程序

%一个求包络线和包络谱的程序%现代希尔伯特变换解调分析：%带通滤波；希尔伯特变换获得信号时域的包络线；用fft变换获得包络谱%如何获得包络线？%信号经希尔伯特变换不能直接得到包络，设信号x的希尔伯特变换为y，则平方和%x.^2+y.^2(或者再开根号，直接取平方和的效果为好)才是信号x的包络。%构造实验数据clearall;closeall;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);holdonplot(t,s);plot(t,y1,'r');title('包络线');holdoff%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)holdonplot((0:N/2-1)/N*fs,p(1:N/2),'g');grid;%包络谱plot((0:N/2-1)/N*fs,p1(1:N/2),'r');%傅里叶变换功率谱holdoff