短时傅里叶变换

原文地址：http://blog.sina.com.cn/s/blog_82a927880102uwi1.html

小波和短时要一起学，至于原因往下看便知。

短时傅立叶变换基本思想是将信号加滑动时间窗，并对窗内信号做傅立叶变换，得到信号的时变频谱。因而它的时间分辨率和频率分辨率受Heisenberg测不准原理约束，一旦窗函数选定，时频分辨率便确定下来。这就使它对突变信号和非平稳信号的分析存在局限性，因而不是一种动态的分析方法, 不能敏感地反映信号的突变，只适用于对缓变信号的分析。
    Wigner一Ville分布定义为信号中心协方差函数的傅立叶变换，它具有许多优良的性能，如对称性、时移性、组合性、复共扼关系等，不会损失信号的幅值与相位信息，对瞬时频率和群延时有清晰的概念。
其不足是不能保证非负性，尤其是对多分量信号或具有复杂调制规律的信号会产生严重的交叉项干扰，这是二次型时频分布的固有结果，大量的交叉项会淹没或严重干扰信号的自项，模糊信号的原始特征。后续的有人对Cohen类中的核函数进行改造，提出了伪winger—ville分布、修正平滑伪Winger—Ville分布等各种各样的新型时频分布，对交叉项干扰的抑制起了较大的作用，但是不含有交叉项干扰且具有Winger—Ville分布聚集性的时频分布是不存在的。
    小波变换通过伸缩和平移运算对信号进行多尺度分解,能够有效地从信号中获取各种时频信息, 它在时域和频域同时具有良好的局部化性质，具有多分辨率分析特性。
但小波分解的结果依赖小波基函数，而各小波基函数的适用范围很不一致，这就造成了小波基选择问题，如果母小波选择不当，则应用效果会大受影响；小波分析不具有自适应性，一旦选择了小波基和分解尺度，则用它来分析具多频率成分的数据时，所得结果只能反映某一固定频带内的信号,所以要选择不同的小波基。

 

闲聊篇：

1、Fourier Transform 缺陷----FT局域化特性分析

   FT在平稳信号分析和处理中有着突出贡献的基本原因在于，人们利用它可以把复杂的时间信号和空间信号转换到频率域中，然后用频谱特性去分析和表示时域信号的特性。

    FT 正变换告诉我们：从时间（模拟）信号中提取频谱信息，就是使用（－∞，∞）的时间信息来计算单个频率的频谱（频域过程() Fω的任一频率组成部分的值），这是由在（－∞，∞）上的时域过程决定的。故可知，傅里叶变换对彼此间的刻画是“全局性”的，不能反映各自局部区域上的特征,人们虽然从傅立叶变换能清楚地看到一个信息包含的每一个频率的多少，但很难看出不同信号的发射时间和发射的延续时间，缺少时间信息使得傅立叶分析变得脆弱而容易失误。

    伊利诺依斯大学教授Y. Me y e r曾说：“若你记录1小时长的信息而在最后5分钟出错，这一错误就会毁了整个傅立叶变换。相位的错误是灾难性的，如果在相位上哪怕犯了一个错误，你最后就会发现你所干的事与最初的信号无关了。”

    实际上，对常见的不平稳信号，如语音信号、音乐信号、核探测的脉冲信号、以及核医学的图像信号等，它们的频域特性是随时间变化的，人们需要了解某些局部时段上所对应的主要频率特性是什么，也需要了解某些频率的信息出现在哪些时段上，也就是需要了解时－频局部化要求。对于这种时－频局部化要求，傅立叶变换是无能为力的。

    解决方案：短时傅里叶变换、小波变换。其各有优劣，本篇聊点短时傅里叶变换，最后会引出一点小波变换。

2、短时傅里叶变换（窗式傅里叶变换）

    基本思想：把非平稳过程看成是一系列短时平稳信号的叠加，短时性可通过在时间上加窗实现。。通过该方法，人们至少可以说，无论发生了什么，它一定是发生在信号的某个特定部分。

    在傅立叶积分中，使用时间窗口函数与信号相乘，实现在u附近的开窗和平移，然后进行傅立叶变换。在线性空间有一个可测的、平方可积的函数,对其进行窗式傅立叶变换：

    窗口傅立叶变换的逆变换式：

3、小波引出

    由窗口傅立叶变换对函数（信号）进行的分析，相当于用一个形状、大小和放大倍数相同的“放大镜”在时－频相平面上移动去观察某固定长度时间内的频率特性。这里的问题是：尽管窗式傅立叶变换能解决变换函数的局域化问题，但是，其窗口的大小和形状是固定的，即窗口没有自适应性。这意味着什么？

    而在实际问题中，对于高频谱信息，由于波形相对要窄，时间间隔要相对的小，以求给出比较好的精度，进而更好地确定峰值和断点，或者说需要用窄的时域窗来反映信息的高频成分；而对于低频谱信息，由于波形相对是宽的，时间段要相对的宽才能给出完整的信号信息，或者说必须用较宽的时域窗来反映信息的低频成分。而用短时傅里叶变换，如果你选择一扇宽窗子，低频成分可以看得清楚，在高频部分确定时间时就很糟糕；若你选一扇窄窗子，在高频可以很好确定时间，但在低频的频率就可能装不进去。

    这样，真正合适的做法是“放大镜”的长宽是可以变化的，正是为了实现这样的目的，人们引进了小波变换。

 实践篇：

    主要有两个程序，分别：短时傅里叶变换、逆短时傅里叶变换

Short-Time Fourier Transform

1: function [stft, f, t] = stft(x, wlen, h, nfft, fs)

2: if size(x,2) > 1

3: x = x';

4: end

5: xlen = length(x);

6: win = hamming(wlen, 'periodic');

7: rown = ceil((1+nfft)/2);

8: coln = 1+fix((xlen-wlen)/h);

9: stft = zeros(rown, coln);

10: indx = 0;

11: col = 1;

12: % perform STFT

13: while indx + wlen <= xlen

14: xw = x(indx+1:indx+wlen).*win;

15: X = fft(xw, nfft);

16: stft(:,col) = X(1:(rown));

17: indx = indx + h;

18: col = col + 1;

19: end

20: t = (wlen/2:h:xlen-wlen/2-1)/fs;

21: f = (0:rown-1)*fs/nfft;

22: end

 

Inverse Short-Time Fourier Transform  

1: function [x, t] = istft(stft, h, nfft, fs)

2: coln = size(stft, 2);

3: xlen = nfft + (coln-1)*h;

4: x = zeros(1, xlen);

5: win = hamming(nfft, 'periodic');

6: if rem(nfft, 2)

7: for b = 0:h:(h*(coln-1))

8: X = stft(:, 1 + b/h);

9: X = [X; conj(X(end:-1:2))];

10: xprim = real(ifft(X));

11: x((b+1):(b+nfft)) = x((b+1):(b+nfft)) + (xprim.*win)';

12: end

13: else

14: for b = 0:h:(h*(coln-1))

15: X = stft(:, 1+b/h);

16: X = [X; conj(X(end-1:-1:2))]; 

17: xprim = real(ifft(X));

18: x((b+1):(b+nfft)) = x((b+1):(b+nfft)) + (xprim.*win)';

19: end

20: end

21: W0 = sum(win.^2);

22: x = x.*h/W0;

23: actxlen = length(x);

24: t = (0:actxlen-1)/fs;

25: end

 

测试信号1：

    这里先是选取比较简单的正弦信号做测试

参数：

1: fs = 48000;

2: t = 0:1/fs:1-1/fs;

3: x = 10*sin(2*pi*t*10);

测试信号2：

    这里选择的是以前我们windows开机时候的那个声音，很熟悉吧？