python中短时傅里叶变换模块短时傅里叶变换原理

转载

goody 2023-12-08 17:05:04

文章标签 python中短时傅里叶变换模块小波分析相关傅里叶变换傅立叶变换时域 文章分类 Python 后端开发

闲聊篇：

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

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

FT 正变换告诉我们：从时间（模拟）信号中提取频谱信息

，就是使用（－∞，∞）的时间信息来计算单个频率的频谱（频域过程() Fω的任一频率组成部分的值），这是由在（－∞，∞）上的时域过程

决定的。故可知，傅里叶变换对彼此间的刻画是“全局性”的，不能反映各自局部区域上的特征,人们虽然从傅立叶变换能清楚地看到一个信息包含的每一个频率的多少，但很难看出不同信号的发射时间和发射的延续时间，缺少时间信息使得傅立叶分析变得脆弱而容易失误。

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

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);