现代测试技术李德胜2.docVIP

现代测试技术第二次作业 李德胜 3112058010 硕2007班 实验目的 通过matlab编程，对短时傅里叶变换，小波包分解和EMD分解详细了解，对信号的处理方法作形象的比较。 二、实验题目 对下列信号进行短时傅里叶变换，小波包分解和EMD分解： 其中 实验结果 波形图程序如下： clc; clear all; close all; Ts=0.0002; N=5000; t=0:0.0002:1; x=zeros(); for i=1:(length(t)) if t(i)0.2 x(i)=8*sin(60*pi*t(i)+pi/12)+4*sin(120*pi*t(i)+pi/2)+sin(300*pi*t(i)+pi/6)+20*sin(2*pi*250*t(i)+pi/3)*exp(-20*(t(i)-0.2)); else x(i)=8*sin(60*pi*t(i)+pi/12)+4*sin(120*pi*t(i)+pi/2)+sin(300*pi*t(i)+pi/6); end end figure(1); plot(t,x); x（t）的波形图如下图。 短时傅里叶变换程序如下： figure(2); spectrogram(x,128,120,1024,1000); 采样频率为1000Hz，点数为1000时，得到的短时傅里叶变换图如下图： 对波形进行小波包分解时，选用db1函数和shannon函数。 程序如下： wpt=wpdec(x,3,db1,shannon); bt=besttree(wpt); plot(bt); plot_hht00(x,Ts,t) 得到的图形如下图。 （0，0）节点的波形图： （1，0）节点的波形图： （1,1）节点的波形图： （2，0）节点的波形图： （2，1）节点的波形图： （2，2）节点的波形图： （2，3）节点的波形图： 利用besttree函数计算可知，（2,3）节点的分解最佳。波形时长1s，在0.2s时有明显的波动，波动信号又随时间逐渐降低。求解图形如下。 EMS分解的程序： % FUNCTIONS function u = ismonotonic(x) %u=0表示x不是单调函数，u=1表示x为单调的 u1 = length(findpeaks(x))*length(findpeaks(-x)); if u1 0, u = 0; else, u = 1; end function u = isimf(x) %u=0表示x不是固有模式函数，u=1表示x是固有模式函数 N = length(x); u1 = sum(x(1:N-1).*x(2:N) 0); u2 = length(findpeaks(x))+length(findpeaks(-x)); if abs(u1-u2) 1, u = 0; else, u = 1; end function s = getspline(x) %三次样条函数拟合成元数据包络线 N = length(x); p = findpeaks(x); s = spline([0 p N+1],[0 x(p) 0],1:N); function n = findpeaks(x) % Find peaks.找到极值 ,n为极值点所在位置 % n = findpeaks(x) n = find(diff(diff(x) 0) 0); u = find(x(n+1) x(n)); n(u) = n(u)+1; function imf = emd(x) x= transpose(x(:));%转置为行矩阵 imf = []; while ~ismonotonic(x) %当x不是单调函数，分解终止条件 x1 = x; sd = Inf;%均值 %直到x1满足IMF条件，得c1 while (sd 0.1) | ~isimf(x1) %当标准偏差系数sd大于0.1或x1不是固有模态函数时，分量终止条件 s1 = getspline(x1);%上包络线 s2 = -getspline(-x1);%下包络线 x2 = x1-(s1+s2)/2;%此处的x2为文章中的h sd = sum((x1-x2).^2)/sum(x1.^2); x1 = x2; end imf{end+1} = x1; x = x-x1; end imf{end+1} = x; functi