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

[tfr,t,f]=tfrstft(y,T,N,hHamming); figure; contour(t,(0:N/2-1)/N*Fs,abs(tfr(1:N/2,:))); xlabel('Time'); ylabel('Frequency(Hz)'); title(['STFT, Hamming(' num2str(L) ')']); % <三、Gabor变换;利用时频分析包进行分析 NG = 160; Q = 80; [tfr,dgr,gam] = tfrgabor(y,NG,Q); X = (1:NG)*Nfft/NG; Y = (0:N/2-1)*Fs/N; figure; contour(X,Y,abs(tfr(1:N/2,:))); xlabel('Time'); ylabel('Frequency(Hz)'); title(['Gabor Representation, N = ' num2str(NG) ' Q = ' num2str(Q)]); % <四、Wigner-Ville分布> T = 1:800; N = 8192; N = size(tfr,1); % 采用Gauss(N+1)窗
式中 gt , ( ) g ( t )e j ，且 || g ( ) || 1 ， || gt , ( ) || 1 ，窗函数 g ( ) 应取对称。 离散信号的短时傅立叶变换 当我们要在计算机上实现一个信号的短时傅立叶变换时，该信号必须是离散 的，且为有限长。设给定的信号为 x(n), n 0,1,..., L 1 ，有
% time instant(s) and number of frequency bins
[tfr,t,f] = tfrwv(y(10000:10800),T,N); figure; contour(t,(0:N/2-1)/N*Fs,abs(tfr(1:N/2,:))); xlabel('Time'); ylabel('Frequency(Hz)'); title('Wigner-Ville Distribution');
二、原理解析
2.1 短时傅里叶变换 STFT
短时傅里叶变换是基于傅里叶变换的线性时频表示。短时傅立叶变换是一种 常用的时频分析方法，它通过时间窗内的一段信号来表示某一时刻的信号特征。 对于 x(t ) L2 ( R) ，有：
STFTx (t , ) x( ) g t*, ( )d x( ) g * ( t )e j d x( ), g ( t )e j
频率是连续的。为了在计算机上实现，应将频率 离散化，令 k
STFTx (m , k ) x n( g *) n ( mN e M )
n j 2 n k
2 k ，则 M
上式将频域的一个周期 2 分成了 M 点， 显然， 上式是一个标准的 M 点 DFT， 若窗函数 g (n) 的宽度正好也是 M 点，那么上式可写成
2.2 Gabor 变换
Gabor 变换是短时 Fourier 变换中当窗函数取为高斯函数时的一种特殊情况。 Gabor 变换定义为
G f (a, b, )
f (t ) g a (t b)eit dt
其中， g a (t )
1 2 a
exp(
t2 ) ，是高斯函数，称为窗函数。Gabor 取 g (t ) 为 4a
一个高斯函数有两个原因：一是高斯函数的 Fourier 变换仍为高斯函数，这使得 Fourier 逆变换也是用窗函数局部化，同时体现了频域的局部化；二是 Gabor 变 换是最优的窗口 Fourier 变换。
2.3 Wigner-Ville 分布
对信号 s(t ) ，其 Wigner Ville 分布定义为：
采用 Windows 关机声音:“C:\WINDOWS\Media\Windows XP 关机.wav” 。这 里对该音频信号进行了重采样，且进行了时间长度的截取，最后使采样频率
Fs 8000Hz ，信号时间长度为 t 2s 。
图 1 原始音频信号时域图
2、功率谱估计
采用周期图法： P( )
a、Gabor 变换，N=80
b、Gabor 变换，N=320
图 6 分辨率理解示意图
% 通道1,取2s数据
f = Fs*(0:halfLength)/Nfft; figure; plot(f,Pyy(1:halfLength+1)); xlabel('Frequency(Hz)'); ylabel('Power Spectrum'); title('Power Spectrum Analysis'); % <二、短时傅里叶变换;利用时频分析包进行分析> L = 85; hHamming = hamming(L); T = 1:Nfft; N = 256; % time instant(s) and number of frequency bins
1 | X N (e j ) |2 N
图 2 语音信号功率谱分析
3、短时 Fourier 变换
这里加的窗为 Hamming 窗，窗宽度为 L 85 。
图 3 短时傅里叶变换
4、Gabor 变换
这里的高斯窗，宽度取为 N 160
图 4 Gabor 变换
5、Wigner-Ville 分布
这里采用整段时序信号中最前面 800 个点的信号进行分析。 从结果可以看出， Wigner-Ville 分布得到了信号分析时较高的频率分辨率。
五、分析讨论
从图 1 和图 2 的分析可知，语音信号的能量主要分布在 1.4s 以前，频率分 布主要在 1700Hz 以下。三种分解方式均是把信号的能量按时间和频率做二维分 布， 即把信号用时频二维基函数的线性组合来描述。而该基函数表示在时间和频 率区间的能量分布情况。 Wigner-Ville 分布具有对称性、时移性等许多优良的性能，对瞬时频率和群 延时有清晰的概念。其不足是由于 Wigner-Ville 分布不是信号的线性变换，因此
存在交叉项干扰。 另外，时频分析都受到 Heisenberg 测不准原理约束，一旦窗函数选定，时 频分辨率便确定下来。 这就使它对突变信号和非平稳信号的分析存在局限性，自 适应能力有所不足。 对分辨率的理解 可以用下图 6 来解释，图 a 和 b 都采用 Gabor 变换，但是窗口的宽度不一 样，图 a 采用窗口的宽度为 N 80 ，图 b 采用的窗口的宽度为 N 320 。显然， 图 a 的时间分辨率要高于图 b，在时频图中表现为沿时间轴方向（即横轴方向）, 更容易分辨出某个具体“原子”信号的时间中心。而图 b 的时间分辨率虽然低， 但是频率分辨率高于图 a，更容易分辨出频域的中心。这体现了受 Heisenberg 测不准原理约束， 我们对时间分辨率和频率分辨率只能取一个折中，不可兼得的 思想，一个的提高必然导致另一个的相对降低，反之亦然。
现代信号处理 实验七
一、作业内容
（一）信号 一段语音信号（一个词或词组，2 秒左右） ，采样频率应在 8 kHz 以上。 （二）要求 分别用短时傅立叶变换、Gabor 变换和 Wigner-Ville 分布进行分析； 根据语音信号的机理和特征，对各种时频分析的结果进行讨论，包括物理意 义（从基函数角度） 、时频分辨率以及方法的特点。
STFTx (m, k )
M 1 n 0
x(n) g
*
nk ， k 0,1,...., M 1 (n mN )WM
（1）
若 g (n) 的宽度小于 M ，那么可将其补零，使之变成 M ，若 g (n) 的宽度大于 （1）式为一标准 DFT，时域、 M ，则应增大 M 使之等于窗函数的宽度。总之， 频域的长度都是 M 点。 式中 N 的大小决定了窗函数沿时间轴移动的间距，N 越小， 上面各式中 m 的 取值越多，得到的时－频曲线越密。 MATLAB 的时频分析 Toolbox 中给出了实现（1）式的程序，即 tfrstft。
图 5 Wigner-Ville 分布
四、实验程序
下面是本次程序代码，采用了Matlab时频分析Toolbox（网上可下载） 。
clear; % <Z、音频信号读入和重采样> [y0, Fs0] = wavread('关机.wav'); Fs = 8000; Ts = 1/Fs; y = resample(y0,Fs,Fs0); % sound(y); y = y(:,1); y = y(1:2*Fs); % <一、功率谱分析> Nfft = length(y); t = (0:Nfft-1)/Fs; figure; plot(t,y); title(['Original Signal, Fs = ' num2str(Fs/1000) 'KHz']); Y = fft(y,Nfft); Pyy= Y.* conj(Y)/Nfft; halfLength = floor(Nfft/2); % fft分析 % 采用系统关机声音 % 重采样，8K Hz
W (t , ) s(t )s* (t )e j d 2 2
Wigner Ville 分布直接对信号的二次型计算傅里叶变换，满足理想时频分布

ቤተ መጻሕፍቲ ባይዱ


的大多数性质。 Wigner Ville 分布的目的是得到信号分析时较高的频率分辨率。
三、实验结果
1、原始信号
STFTx (m, e j ) x(n) g * (n mN )e jn
n
x(n), g (n mN )e jn