现代数字信号处理课程报告

1、现代数字信号处理课程报告 问题：需要利用4种时频分析方法对给定的6个信号文件进行处理，给出每种分析方法的结果图。由于时频分析方法众多，选取的用于MATLAB编程的4种时频分析方法分别为：短时傅立叶变换（STFT）、连续小波变换、Wigner-Ville分布、希尔伯特黄变换（HHT）。 短时傅立叶变换： 傅立叶分析一直被认为是最完美的数学理论和最实用的方法之一。但是用傅立叶分析只能获得信号的整个频谱，而难以获得信号的局部特性，特别是对于突变信号和非平稳信号难以获得希望的结果。 为了克服经典傅立叶分析本身的弱点，1946年Gabor提出了短时傅立叶变换，其基本思想是首先利用窗函数截取信号，假设信号

2、在各个窗内都是平稳的，再对第一个窗内的信号采用傅立叶变换进行分析，以确定此时窗内存在的频率成分，然后沿着信号时间方向移动窗函数，得到频率随时间的变换关系，即所需要的时频分布。 解题思路与流程： 分别读取6个信号，对读取的信号，计算其长度，从而选取窗口长度为100，因此重叠宽度为50，从而计算需要分段的数目，进而计算做fft的点数，此处采取的窗函数为矩形窗，沿着信号时间方向移动窗函数，从而求出每段的fft，拼接后画出时频分析二维与三维图。6个信号的图如图1-6所示。 NACB.BHR信号STFT nf a. NACB.BHR三维b. NACB.BHR二维 图1NACB.BHR nf a. NAC

3、B.BHT三维b. NACB.BHT二维 图 2NACB.BHT NACB.BHZ信号STFT nf a. NACB.BHZ三维b. NACB.BHZ二维 图3NACB.BHZ nf b. YHNB.BHR三维b. YHNB.BHR二维 图4YHNB.BHR YHNB.BHT信号 STFT nf c. YHNB.BHT三维b. YHNB.BHT二维 图5YHNB.BHT nf d. YHNB.BHZ三维b. YHNB.BHZ二维 图 6YHNB.BHZ 连续小波变换： 由于STFT采用的滑动窗函数一经选定就固定不变，故决定了其时频分辨率固定不变，不具备自适应能力，因此短时傅立叶变换没有从根本上

4、解决傅立叶分析的固有问题。小波变换的诞生，正是为了克服经典傅立叶分析本身的不足。 小波变换的基本思想是通过伸缩平移运算对信号逐步进行多尺度细化，最终达到高频处时间细分和低频处频率细分，从而聚焦到信号的任意细节。 解题思路与流程： 由于小波分析方法是一种窗口面积固定但其形状可以改变的方法，因此我们弹性的时频窗，将信号分为了128个尺度，不同的尺度对应不同的频率，频率高的地方，需要时间间隔小，因此高频区会有一个较高的时间分辨率，而低频的地方会有较高的频率分辨率，这样就保证了时间分辨率和频率分辨率在各自需要的范围都达到很高的精度。 在本次实验中选取总尺度为128，选取sym2函数为小波函数，计算出小

5、波基的中心频率，根据中心频率与总尺度计算出各个尺度划分，从而求出小波系数，进而画出如图7-12的小波时频图。 NACB.BHR信号sym2小波时频图 时间 t/s-6频率 f/Hz a. NACB.BHR三维b. NACB.BHR二维 图7NACB.BHR 时间 t/s-6频率 f/Hz a. NACB.BHT三维b. NACB.BHT二维 图 8NACB.BHT NACB.BHZ信号sym2小波时频图 时间 t/s-6频率 f/Hz a. NACB.BHZ三维b. NACB.BHZ二维 图9NACB.BHZ 时间 t/s-5频率 f/Hz a. YHNB.BHR三维b. YHNB.BHR二维

6、 图10YHNB.BHR YHNB.BHT信号sym2 小波时频图 时间 t/s-5频率 f/Hz a. YHNB.BHT三维b. YHNB.BHT二维 图11YHNB.BHT 时间 t/s-5频率 f/Hz a. YHNB.BHZ三维b. YHNB.BHZ二维 图 12YHNB.BHZ Wigner-Ville分布： Wigner-Ville distribution (WVD)方法，是由Wigner1932提出并由Ville于1948年引入信号分析领域。WVD可看作信号能量在联合的时间和频率域中的分布，其功能是分析时间序列中的频率成分及其能量，是分析非平稳和时变信号的重要工具。 WVD基本

7、思想为信号中心协方差函数的傅立叶变换，它具有许多优良的性能，如对称性、时移性、组合性、复共轭关系等，不会损失信号的幅值与相位信息，对瞬时频率和群延时有清晰的概念。其不足是不能保证非负性，尤其是对多分量信号或具有复杂调制规律的信号会产生严重的交叉项干扰，这是二次型时频分布的固有结果，大量的交叉项会淹没或严重干扰信号的自项，模糊信号的原始特征。 解题思路与流程： 对于x(t)为一连续时间信号，则公式（1）称为信号x(t)的自WVD。 ? ?,? = ?(?+?/2)?(?/2)exp?(?)?（1） 式（1）中，?(?)为?(?)的复共轭，?为频率轴划分区间，值为2的n次方，因此对x(t)信号进行

8、处理时，可以先计算信号的?(?+?/2)?(?/2)，然后对结果取fft，从而画出如图13-18所示时频分析图。 图13NACB.BHR图14NACB.BHT 图15NACB.BHZ图16YHNB.BHR 图17 YHNB.BHT图18YHNB.BHZ HHT变换： HHT就是Hilbert-Huang Transform（希尔伯特黄变换），是由美国航天航空局黄鳄教授于1998年提出。HHT方法从信号自身特征出发，用经验模态分解（EMD）方法把信号分解成一系列的本征模态函数（IMF），然后对这些IMF分量进行Hilbert变换，从而得到时频平面上能量分布的Hilbert谱图，打破了测不准原理的

9、限制，可以准确地表达信号在时频面上的各类信息。 HHT的功能是可以分析并提取时间序列中的主要频率（即经验模态）成分。其特点在于可以提取（而不仅仅是分析）信号中的主要成分。HHT从数学理论上很完美，但在应用时常常受到攻击：其所谓经验模态缺乏物理本质或意义。 图 21 NACB.BHZ 图22 YHNB.BHR 图23 YHNB.BHT图24 YHNB.BHZ 结果分析： 根据图1.b-6.b（利用STFT处理6个信号二维时频图）可知，图1.b信号的时间分辨率较高，但是频率分辨率能力差；图2.b信号的频率分辨率较高，而时间分辨率较差，这是因为图1.b在利用短时傅立叶变换时，选取的窗口长度（Nw=2

10、0）过短，导致有较高的时间分辨率，但是频率分辨率能力差；图2.b选取的窗口长度（Nw=1000）过长，导致有较高的频率分辨率，但是时间分辨率能力差。 根据图7.b-12.b（利用连续小波变换处理6个信号二维时频图）可知，连续小波变换能给出比较好的时间精度，在低频区也有很好的频率精度，但是在高频区频率分辨能力较弱。解决了STFT时间分辨率差的问题。 根据图13-18（利用WVD处理6个信号二维时频图）可知，WVD方法在低频区有很好的时间分辨率和频率分辨率，即WVD有很好的时频聚焦性。但是式（1）中，?(?+?/2)?(?/2)导致信号x(t)出现了两次，而且WVD里不含有任何窗函数，因此必然会产

11、生交叉项，而这些交叉项会在信号的各个分量之间产生虚假信号，导致分析不准确。 根据图19-24（利用HHT处理6个信号二维时频图）可知，HHT变换有很高的时间分辨率和频率分辨率，然而EMD方法还不成熟，不知道为什么求出来的二维时频谱只有离散的一些点（只画出了其中一个信号，如下图）。猜想是因为在低频区出现了一些不合理的频率成分，掩盖了低能量的频率成分。 二维hilbert时-频-幅度谱 6-654频率 f/Hz32 1 200400600800 n1000120014001600 程序代码： STFT： clear all； clc; close all; fid=fopen(NACB.BHR,r

12、,ieee-le); A=fread(fid,70,1,float32); B=fread(fid,40,1,int32); C=char(fread(fid,1,192,char); x1=fread(fid,float32); A(A=-12345.0)=NaN; B(B=-12345)=NaN; fclose(fid); figure(1); timefreq(x1,20,rec); title(NACB.BHR信号STFT); fid=fopen(NACB.BHT,r,ieee-le); A=fread(fid,70,1,float32); B=fread(fid,40,1,int32

13、); C=char(fread(fid,1,192,char); x2=fread(fid,float32); A(A=-12345.0)=NaN; B(B=-12345)=NaN; fclose(fid); figure(2); timefreq(x2,1000,rec); title(NACB.BHT信号STFT); fid=fopen(NACB.BHZ,r,ieee-le); A=fread(fid,70,1,float32); B=fread(fid,40,1,int32); C=char(fread(fid,1,192,char); x3=fread(fid,float32); A(

14、A=-12345.0)=NaN; B(B=-12345)=NaN; fclose(fid); figure(3); timefreq(x3,100,rec); title(NACB.BHZ信号STFT); fid=fopen(YHNB.BHR,r,ieee-le); A=fread(fid,70,1,float32); B=fread(fid,40,1,int32); C=char(fread(fid,1,192,char); x4=fread(fid,float32); A(A=-12345.0)=NaN; B(B=-12345)=NaN; fclose(fid); figure(4); t

15、imefreq(x4,100,rec); title(YHNB.BHR信号STFT); fid=fopen(YHNB.BHT,r,ieee-le); A=fread(fid,70,1,float32); B=fread(fid,40,1,int32); C=char(fread(fid,1,192,char); x5=fread(fid,float32); A(A=-12345.0)=NaN; B(B=-12345)=NaN; fclose(fid); figure(5); timefreq(x5,100,rec); title(YHNB.BHT信号STFT); fid=fopen(YHNB.

16、BHZ,r,ieee-le); A=fread(fid,70,1,float32); B=fread(fid,40,1,int32); C=char(fread(fid,1,192,char); x6=fread(fid,float32); A(A=-12345.0)=NaN; B(B=-12345)=NaN; fclose(fid); figure(6); timefreq(x6,100,rec); title(YHNB.BHZ信号STFT); 函数timefreq： function timefreq(x,Nw,window) Lap=Nw/2;%重叠宽度 Tn=floor(length(

17、x)-Lap)/(Nw-Lap);%计算分段数目 nfft=3ceil(log2(Nw);%做fft的点数 TF=zeros(Tn,nfft); %时频矩阵 for i=1:Tn if(strcmp(window,rec) Xw=x(i-1)*10+1:i*10+10);%加窗矩形处理 else return; end temp=fft(Xw,nfft);%求fft temp=fftshift(temp);%调整0频位置 TF(i,:)=temp;%保存分段fft结果 end %绘制时频分析结果 fnew=(1:nfft)-nfft/2)/nfft; tnew=(1:Tn)*Lap; % F,

18、T=meshgrid(fnew,tnew); % mesh(T,F,abs(TF); %画二维图 imagesc(tnew,fnew,abs(TF);set(gca,YDir,normal);%画三维图 xlabel(n);ylabel(f);zlabel(Gf); xiaobo： clear all； clc; close all; fid=fopen(NACB.BHR,r,ieee-le); A=fread(fid,70,1,float32); B=fread(fid,40,1,int32); C=char(fread(fid,1,192,char); data1=fread(fid,fl

19、oat32); A(A=-12345.0)=NaN; B(B=-12345)=NaN; fclose(fid); L=length(data1); %数据长度 t=1:L; totalscal=128; Fc=centfrq(sym2);% 小波基的中心频率 scals=(2*Fc*totalscal)./(1:totalscal);%尺度 f=scal2frq(scals,sym2,1); %将尺度转换为频率 coefs=cwt(data1,scals,sym2);%求小波系数 figure; %画小波时频图 % imagesc(t,f,abs(coefs);set(gca,YDir,nor

20、mal);%二维图 surf(t,f,abs(coefs),edgecolor,none);%三维图 colorbar;xlabel(时间 t/s);ylabel(频率 f/Hz);title(NACB.BHR信号sym2小波时频图); fid=fopen(NACB.BHT,r,ieee-le); A=fread(fid,70,1,float32); B=fread(fid,40,1,int32); C=char(fread(fid,1,192,char); data2=fread(fid,float32); A(A=-12345.0)=NaN; B(B=-12345)=NaN; fclose

21、(fid); L=length(data2); %数据长度 t=1:L; totalscal=128; Fc=centfrq(sym2);% 各种小波基的中心频率 scals=(2*Fc*totalscal)./(1:totalscal);%尺度 f=scal2frq(scals,sym2,1); %将尺度转换为频率 coefs=cwt(data2,scals,sym2);%求小波系数 figure; %画小波时频图 surf(t,f,abs(coefs),edgecolor,none); colorbar;xlabel(时间 t/s);ylabel(频率 f/Hz);title(NACB.B

22、HT信号sym2小波时频图); fid=fopen(NACB.BHZ,r,ieee-le); A=fread(fid,70,1,float32); B=fread(fid,40,1,int32); C=char(fread(fid,1,192,char); data3=fread(fid,float32); A(A=-12345.0)=NaN; B(B=-12345)=NaN; fclose(fid); L=length(data3); %数据长度 t=1:L; totalscal=128; Fc=centfrq(sym2);% 各种小波基的中心频率 scals=(2*Fc*totalscal

23、)./(1:totalscal);%尺度 f=scal2frq(scals,sym2,1); %将尺度转换为频率 coefs=cwt(data3,scals,sym2);%求小波系数 figure; %画小波时频图 surf(t,f,abs(coefs),edgecolor,none); colorbar;xlabel(时间 t/s);ylabel(频率 f/Hz);title(NACB.BHZ信号sym2小波时频图); fid=fopen(YHNB.BHR,r,ieee-le); A=fread(fid,70,1,float32); B=fread(fid,40,1,int32); C=ch

24、ar(fread(fid,1,192,char); data4=fread(fid,float32); A(A=-12345.0)=NaN; B(B=-12345)=NaN; fclose(fid); L=length(data4); %数据长度 t=1:L; totalscal=128; Fc=centfrq(sym2);% 各种小波基的中心频率 scals=(2*Fc*totalscal)./(1:totalscal);%尺度 f=scal2frq(scals,sym2,1); %将尺度转换为频率 coefs=cwt(data4,scals,sym2);%求小波系数 figure; %画小

25、波时频图 surf(t,f,abs(coefs),edgecolor,none); colorbar;xlabel(时间 t/s);ylabel(频率 f/Hz);title(YHNB.BHR信号sym2小波时频图); fid=fopen(YHNB.BHT,r,ieee-le); A=fread(fid,70,1,float32); B=fread(fid,40,1,int32); C=char(fread(fid,1,192,char); data5=fread(fid,float32); A(A=-12345.0)=NaN; B(B=-12345)=NaN; fclose(fid); L=

26、length(data5); %数据长度 t=1:L; totalscal=128; Fc=centfrq(sym2);% 各种小波基的中心频率 scals=(2*Fc*totalscal)./(1:totalscal);%尺度 f=scal2frq(scals,sym2,1); %将尺度转换为频率 coefs=cwt(data5,scals,sym2);%求小波系数 figure; %画小波时频图 surf(t,f,abs(coefs),edgecolor,none); colorbar; xlabel(时间 t/s); ylabel(频率 f/Hz); title(YHNB.BHT信号sy

27、m2小波时频图); fid=fopen(YHNB.BHZ,r,ieee-le); A=fread(fid,70,1,float32); B=fread(fid,40,1,int32); C=char(fread(fid,1,192,char); data6=fread(fid,float32); A(A=-12345.0)=NaN; B(B=-12345)=NaN; fclose(fid); L=length(data6); %数据长度 t=1:L; totalscal=128; Fc=centfrq(sym2);% 各种小波基的中心频率 scals=(2*Fc*totalscal)./(1:

28、totalscal);%尺度 f=scal2frq(scals,sym2,1); %将尺度转换为频率 coefs=cwt(data6,scals,sym2);%求小波系数 figure; %画小波时频图 surf(t,f,abs(coefs),edgecolor,none); colorbar;xlabel(时间 t/s);ylabel(频率 f/Hz);title(YHNB.BHZ信号sym2小波时频图); WVD： clc; clear all; close all; fid=fopen(NACB.BHR,r,ieee-le); A=fread(fid,70,1,float32); B=f

29、read(fid,40,1,int32); C=char(fread(fid,1,192,char); HR1=fread(fid,float32); A(A=-12345.0)=NaN; B(B=-12345)=NaN; fclose(fid); WVD,t,f=WignerVille(HR1); set(gcf,Position,20 100 400 300); mesh(t,f,abs(WVD); colorbar; axis(min(t) max(t) min(f) max(f); ylabel(f/Hz); xlabel(t/s); title(NACB.BHR信号WVD); fid

30、=fopen(NACB.BHZ,r,ieee-le); A=fread(fid,70,1,float32); B=fread(fid,40,1,int32); C=char(fread(fid,1,192,char); HR2=fread(fid,float32); A(A=-12345.0)=NaN; B(B=-12345)=NaN; fclose(fid); WVD,t,f=WignerVille(HR2); mesh(t,f,abs(WVD); colorbar; axis(min(t) max(t) min(f) max(f); ylabel(f/Hz); xlabel(t/s); t

31、itle(NACB.BHZ信号WVD); fid=fopen(NACB.BHT,r,ieee-le); A=fread(fid,70,1,float32); B=fread(fid,40,1,int32); C=char(fread(fid,1,192,char); HR3=fread(fid,float32); A(A=-12345.0)=NaN; B(B=-12345)=NaN; fclose(fid); WVD,t,f=WignerVille(HR3); mesh(t,f,abs(WVD); colorbar; axis(min(t) max(t) min(f) max(f); ylab

32、el(f/Hz); xlabel(t/s); title(NACB.BHT信号WVD); fid=fopen(YHNB.BHZ,r,ieee-le); A=fread(fid,70,1,float32); B=fread(fid,40,1,int32); C=char(fread(fid,1,192,char); HR4=fread(fid,float32); A(A=-12345.0)=NaN; B(B=-12345)=NaN; fclose(fid); WVD,t,f=WignerVille(HR4); mesh(t,f,abs(WVD); colorbar; axis(min(t) ma

33、x(t) min(f) max(f); ylabel(f/Hz); xlabel(t/s); title(YHNB.BHZ信号WVD); fid=fopen(YHNB.BHT,r,ieee-le); A=fread(fid,70,1,float32); B=fread(fid,40,1,int32); C=char(fread(fid,1,192,char); HR5=fread(fid,float32); A(A=-12345.0)=NaN; B(B=-12345)=NaN; fclose(fid); WVD,t,f=WignerVille(HR5); mesh(t,f,abs(WVD);

34、colorbar; axis(min(t) max(t) min(f) max(f); ylabel(f/Hz); xlabel(t/s); title(YHNB.BHT信号WVD); fid=fopen(YHNB.BHR,r,ieee-le); A=fread(fid,70,1,float32); B=fread(fid,40,1,int32); C=char(fread(fid,1,192,char); HR6=fread(fid,float32); A(A=-12345.0)=NaN; B(B=-12345)=NaN; fclose(fid); WVD,t,f=WignerVille(H

35、R6); mesh(t,f,abs(WVD); colorbar; axis(min(t) max(t) min(f) max(f); ylabel(f/Hz); xlabel(t/s); title(YHNB.BHR信号WVD); 函数WignerVille： function WVD,t,f = WignerVille( Sig, SampFreq,FreqBins ) % Sig :输入信号 % FreqBins：频率轴划分区间数（默认为512） % SampFreq：信号的采样频率 % WVD ：计算结果 if (nargin<1), error(At least one input required !); end; Sig=hilbert(real(Sig); %计算信号的解析信号 SigLen=length(Sig); %获取信号的长度 if (nargin<3), FreqBin