短时傅立叶变换_Gabor变换和Wigner-Ville分布实验.pdf

《短时傅立叶变换_Gabor变换和Wigner-Ville分布实验.pdf》由会员上传分享，免费在线阅读，更多相关内容在教育资源-天天文库。

1、现代信号处理实验七一、作业内容（一）信号一段语音信号（一个词或词组，2秒左右），采样频率应在8kHz以上。（二）要求分别用短时傅立叶变换、Gabor变换和Wigner-Ville分布进行分析；根据语音信号的机理和特征，对各种时频分析的结果进行讨论，包括物理意义（从基函数角度）、时频分辨率以及方法的特点。二、原理解析2.1短时傅里叶变换STFT短时傅里叶变换是基于傅里叶变换的线性时频表示。短时傅立叶变换是一种常用的时频分析方法，它通过时间窗内的一段信号来表示某一时刻的信号特征。2对于x(t)L(R)，有：*STFT(t,)x()g()dxt,*jjx()

2、g(t)edx(),g(t)ej式中g()g(t)e，且

3、

4、g()

5、

6、1，

7、

8、g()

9、

10、1，窗函数g()应取对称。t,t,离散信号的短时傅立叶变换当我们要在计算机上实现一个信号的短时傅立叶变换时，该信号必须是离散的，且为有限长。设给定的信号为x(n),n0,1,...,L1，有j*jnSTFTmex(,)xngnmNe()()njnxngnmNe(),()式中N是在时间轴上窗函数移动的步长，是圆周频率，T，T为由ssx(t)得到x(n)的抽样间隔。该式对应傅立叶变换中的DTFT，即时间是离散的，2频率是连

11、续的。为了在计算机上实现，应将频率离散化，令k，则kM2jnkSTFTm(,)xngnmNe()*(M)xkn上式将频域的一个周期2分成了M点，显然，上式是一个标准的M点DFT，若窗函数g(n)的宽度正好也是M点，那么上式可写成M1*nkSTFTx(m,k)x(n)g(nmN)WM，k0,1,....,M1（1）n0若g(n)的宽度小于M，那么可将其补零，使之变成M，若g(n)的宽度大于M，则应增大M使之等于窗函数的宽度。总之，（1）式为一标准DFT，时域、频域的长度都是M点。式中N的大小决定了窗函数沿时间轴移动的间距，N越小，上面各式中m的取值越多

12、，得到的时－频曲线越密。MATLAB的时频分析Toolbox中给出了实现（1）式的程序，即tfrstft。2.2Gabor变换Gabor变换是短时Fourier变换中当窗函数取为高斯函数时的一种特殊情况。Gabor变换定义为itGab(,,)ftgtbe()()dtfa21t其中，g(t)exp()，是高斯函数，称为窗函数。Gabor取gt()为a2a4a一个高斯函数有两个原因：一是高斯函数的Fourier变换仍为高斯函数，这使得Fourier逆变换也是用窗函数局部化，同时体现了频域的局部化；二是Gabor变换是最优的窗口Fourier变换。2.3Wigne

13、r-Ville分布对信号st()，其WignerVille分布定义为：*jWt(,)st()(st)ed22WignerVille分布直接对信号的二次型计算傅里叶变换，满足理想时频分布的大多数性质。WignerVille分布的目的是得到信号分析时较高的频率分辨率。三、实验结果1、原始信号采用Windows关机声音:“C:WINDOWSMediaWindowsXP关机.wav”。这里对该音频信号进行了重采样，且进行了时间长度的截取，最后使采样频率Fs8000Hz，信号时间长度为ts2。图1原始音频信号时域图2、功率谱估计1j2采用周期图法：

14、P()

15、X(e)

16、NN图2语音信号功率谱分析3、短时Fourier变换这里加的窗为Hamming窗，窗宽度为L85。图3短时傅里叶变换4、Gabor变换这里的高斯窗，宽度取为N160图4Gabor变换5、Wigner-Ville分布这里采用整段时序信号中最前面800个点的信号进行分析。从结果可以看出，Wigner-Ville分布得到了信号分析时较高的频率分辨率。图5Wigner-Ville分布四、实验程序下面是本次程序代码，采用了Matlab时频分析Toolbox（网上可下载）。clear;%[y0,Fs0]=wavread('关机.wav');%

17、采用系统关机声音Fs=8000;Ts=1/Fs;y=resample(y0,Fs,Fs0);%重采样，8KHz%sound(y);y=y(:,1);y=y(1:2*Fs);%通道1,取2s数据%<一、功率谱分析>Nfft=length(y);t=(0:Nfft-1)/Fs;figure;plot(t,y);title(['OriginalSignal,Fs='num2str(Fs/1000)'KHz']);Y=fft(y,Nfft);%fft分析Pyy=Y