小波前奏--短时傅里叶变换

 闲聊篇：

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开机时候的那个声音，很熟悉吧？