最近复现论文需要用到短时傅里叶变换的代码，在百度谷歌都找了一番，没有找到合适的代码，于是scipy文档内的stft()函数翻译了一遍。
一、scipy.signal.stft函数介绍
scipy.signal.stft(x,fs=1.0,window='hann',nperseg=256,noverlap=None,nfft=None,
detrend=False,return_onesided=True,boundary='zeros',padded=True,axis=-1)
短时傅里叶变换可以用于量化非平稳信号频率和相位随时间的变化
参数说明：
x: 数组类型
时间序列测量值
fs: 浮点型，可选
x时间序列的采样频率，默认1.0。
window ： 字符串或元祖或数组，可选
需要使用的窗。如果window是一个字符串或元组，则传递给它window是数组类型，直接以其为窗，其长度必须是nperseg。默认为Hann窗。
nperseg ： int，可选
每个段的长度。默认为256。
noverlap ： int，可选
段之间重叠的点数。如果没有，noverlap = nperseg // 2 .默认为无。如有提前指定时，必须满足COLA约束。
nfft ： int，可选
如果需要零填充FFT，则为使用FFT的长度。如果为 None，则FFT长度为nperseg。默认为无。
detrend ： str或function或False，可选
指定如何去除每个段的趋势。如果类型参数传递给False，则不进行去除趋势。默认为False。
return_onesided ： bool，可选
如果为True，则返回实际数据的单侧频谱。如果 False返回双侧频谱。默认为 True。请注意，对于复杂数据，始终返回双侧频谱。
boundary ： str或None，可选
指定输入信号是否在两端扩展，以及如何生成新值，以使第一个窗口段在第一个输入点上居中。这具有当所采用的窗函数从零开始时能够重建第一输入点的益处。有效选项是['even', 'odd', 'constant', 'zeros', None].默认为‘zeros’,对于补零操作[1, 2, 3, 4]变成[0, 1, 2, 3, 4, 0] 当nperseg=3.
填充： bool，可选
指定输入信号在末尾是否填充零以使信号精确地拟合为整数个窗口段，以便所有信号都包含在输出中。默认为True。填充发生在边界扩展之后，如果边界不是None，则填充为True，默认情况下也是如此。
axis ： int，可选
计算STFT的轴; 默认值超过最后一个轴(即axis=-1)。
返回值
f： ndarray
采样频率。
t： ndarray
时间。
Zxx： ndarray
x的 STFT 。默认情况下，Zxx的最后一个轴对应于段时间。
from scipy import signal
import matplotlib.pyplot as plt
import numpy as np
# 产生一个测试信号，振幅为2的正弦波，其频率在3kHZ缓慢调制，振幅以指数形式下降的白噪声
fs = 10e3
N = 1e5
amp = 2 * np.sqrt(2)
noise_power = 0.01 * fs / 2
time = np.arange(N)/float(fs)
mod = 500 * np.cos(2*np.pi*0.25*time)
carrier = amp * np.sin(2*np.pi*3e3*time+mod)
noise = np.random.normal(scale=np.sqrt(noise_power), size=time.shape)
noise *= np.exp(-time/5)
x = carrier + noise
# 计算并绘制STFT的大小
f, t, Zxx = signal.stft(x, fs, nperseg=1000)
plt.pcolormesh(t, f, np.abs(Zxx), vmin = 0, vmax = amp)
plt.title('STFT Magnitude')
plt.ylabel('Frequency [Hz]')
plt.xlabel('Time [sec]')
plt.show()

原理
短时傅里叶变换(Short Time Fourier Transform, STFT) 是一个用于语音信号处理的通用工具.它定义了一个非常有用的时间和频率分布类, 其指定了任意信号随时间和频率变化的复数幅度. 实际上,计算短时傅里叶变换的过程是把一个较长的时间信号分成相同长度的更短的段, 在每个更短的段上计算傅里叶变换, 即傅里叶频谱.
短时傅里叶变换通常的数学定义如下:
其中,
DTFT (Decrete Time Fourier Transform) 为离散时间傅里叶变换.  其数学公式, 如下所示:
其中,  x(n) 为在采样数 n 处的信号幅度. ω~ 的定义如下:
实现时, 短时傅里叶变换被计算为一系列加窗数据帧的快速傅里叶变换 (Fast Fourier Transform, FFT)，其中窗口随时间 “滑动” (slide) 或“跳跃” (hop) 。
Python 实现
在程序中, frame_size 为将信号分为较短的帧的大小, 在语音处理中, 通常帧大小在 20ms 到 40ms 之间. 这里设置为 25ms, 即 frame_size = 0.025;
frame_stride 为相邻帧的滑动尺寸或跳跃尺寸, 通常帧的滑动尺寸在 10ms 到 20ms 之间, 这里设置为 10ms, 即 frame_stride = 0.01. 此时, 相邻帧的交叠大小为 15ms;
窗函数采用汉明窗函数 (Hamming Function) ;
在每一帧, 进行 512 点快速傅里叶变换, 即 NFFT = 512. 具体程序如下:
#-*- coding: utf8 -*-
importnumpy as npdef calc_stft(signal, sample_rate=16000, frame_size=0.025, frame_stride=0.01, winfunc=np.hamming, NFFT=512):#Calculate the number of frames from the signal
frame_length = frame_size *sample_rate
frame_step= frame_stride *sample_rate
signal_length=len(signal)
frame_length=int(round(frame_length))
frame_step=int(round(frame_step))
num_frames= 1 + int(np.ceil(float(np.abs(signal_length - frame_length)) /frame_step))#zero padding
pad_signal_length = num_frames * frame_step +frame_length
z= np.zeros((pad_signal_length -signal_length))#Pad signal to make sure that all frames have equal number of samples
#without truncating any samples from the original signal
pad_signal =np.append(signal, z)#Slice the signal into frames from indices
indices = np.tile(np.arange(0, frame_length), (num_frames, 1)) +np.tile(np.arange(0, num_frames* frame_step, frame_step), (frame_length, 1)).T
frames= pad_signal[indices.astype(np.int32, copy=False)]#Get windowed frames
frames *=winfunc(frame_length)#Compute the one-dimensional n-point discrete Fourier Transform(DFT) of
#a real-valued array by means of an efficient algorithm called Fast Fourier Transform (FFT)
mag_frames =np.absolute(np.fft.rfft(frames, NFFT))#Compute power spectrum
pow_frames = (1.0 / NFFT) * ((mag_frames) ** 2)returnpow_framesif __name__ == ‘__main__‘:importscipy.io.wavfileimportmatplotlib.pyplot as plt#Read wav file
#"OSR_us_000_0010_8k.wav" is downloaded from http://www.voiptroubleshooter.com/open_speech/american.html
sample_rate, signal = scipy.io.wavfile.read("OSR_us_000_0010_8k.wav")#Get speech data in the first 2 seconds
signal = signal[0:int(2. *sample_rate)]#Calculate the short time fourier transform
pow_spec =calc_stft(signal, sample_rate)
plt.imshow(pow_spec)
plt.tight_layout()
plt.show()
参考资料
1. DISCRETE TIME FOURIER TRANSFORM (DTFT). https://www.dsprelated.com/freebooks/mdft/Discrete_Time_Fourier_Transform.html
2. THE SHORT-TIME FOURIER TRANSFORM. https://www.dsprelated.com/freebooks/sasp/Short_Time_Fourier_Transform.html
3. Short-time Fourier transform. https://en.wikipedia.org/wiki/Short-time_Fourier_transform
4. Speech Processing for Machine Learning: Filter banks, Mel-Frequency Cepstral Coefficients (MFCCs) and What‘s In-Between. https://haythamfayek.com/2016/04/21/speech-processing-for-machine-learning.html

一、傅里叶变换

关于傅里叶变换的基本概念在此我就不再赘述了，默认大家现在正处在理解了傅里叶但还没理解小波的道路上。（在第三节小波变换的地方我会再形象地讲一下傅里叶变换）


下面我们主要将傅里叶变换的不足。即我们知道傅里叶变化可以分析信号的频谱，那么为什么还要提出小波变换？答案就是方沁园所说的，“对非平稳过程，傅里叶变换有局限性”。看如下一个简单的信号：做完FFT（快速傅里叶变换）后，可以在频谱上看到清晰的四条线，信号包含四个频率成分。做完FFT（快速傅里叶变换）后，可以在频谱上看到清晰的四条线，信号包含四个频率成分。


一切没有问题。但是，如果是非平稳信号呢？


如上图，最上边的是频率始终不变的平稳信号。而下边两个则是频率随着时间改变的非平稳信号，它们同样包含和最上信号相同频率的四个成分。

做FFT后，我们发现这三个时域上有巨大差异的信号，频谱却非常一致。尤其是下边两个非平稳信号，我们从频域上无法区分它们，因为它们包含的四个频率的信号的成分确实是一样的，只是出现的先后顺序不同。


可见，傅里叶变换处理非平稳信号有天生缺陷。它只能获取一段信号总体上包含哪些频率的成分，但是对各成分出现的时刻并无所知。因此时域相差很大的两个信号，可能频谱图一样。


然而平稳信号大多是人为制造出来的，自然界的大量信号几乎都是非平稳的，所以在比如生物医学信号分析等领域的papers中，基本看不到单纯傅里叶变换这样naive的方法。
上图所示的是一个正常人的事件相关电位。对于这样的非平稳信号，只知道包含哪些频率成分是不够的，我们还想知道上图所示的是一个正常人的事件相关电位。对于这样的非平稳信号，只知道包含哪些频率成分是不够的，我们还想知道各个成分出现的时间。知道信号频率随时间变化的情况，各个时刻的瞬时频率及其幅值——这也就是时频分析。



二、短时傅里叶变换（Short-time Fourier Transform, STFT）

一个简单可行的方法就是——加窗。我又要套用方沁园同学的描述了，“把整个时域过程分解成无数个等长的小过程，每个小过程近似平稳，再傅里叶变换，就知道在哪个时间点上出现了什么频率了。”这就是短时傅里叶变换。

看图：
时域上分成一段一段做FFT，不就知道频率成分随着时间的变化情况了吗！时域上分成一段一段做FFT，不就知道频率成分随着时间的变化情况了吗！

用这样的方法，可以得到一个信号的时频图了：
——此图像来源于“THE WAVELET TUTORIAL” ——此图像来源于“THE WAVELET TUTORIAL”

图上既能看到300Hz, 200 Hz, 100 Hz, 50 Hz四个频域成分，还能看到出现的时间。两排峰是对称的，所以大家只用看一排就行了。


是不是棒棒的？时频分析结果到手。但是STFT依然有缺陷。


使用STFT存在一个问题，我们应该用多宽的窗函数？

窗太宽太窄都有问题：
窗太窄，窗内的信号太短，会导致频率分析不够精准，频率分辨率差。窗太宽，时域上又不够精细，时间分辨率低。窗太窄，窗内的信号太短，会导致频率分析不够精准，频率分辨率差。窗太宽，时域上又不够精细，时间分辨率低。

（这里插一句，这个道理可以用海森堡不确定性原理来解释。类似于我们不能同时获取一个粒子的动量和位置，我们也不能同时获取信号绝对精准的时刻和频率。这也是一对不可兼得的矛盾体。我们不知道在某个瞬间哪个频率分量存在，我们知道的只能是在一个时间段内某个频带的分量存在。 所以绝对意义的瞬时频率是不存在的。）


看看实例效果吧：
——此图像来源于“THE
 WAVELET TUTORIAL” ——此图像来源于“THE WAVELET TUTORIAL”

上图对同一个信号（4个频率成分）采用不同宽度的窗做STFT，结果如右图。用窄窗，时频图在时间轴上分辨率很高，几个峰基本成矩形，而用宽窗则变成了绵延的矮山。但是频率轴上，窄窗明显不如下边两个宽窗精确。


所以窄窗口时间分辨率高、频率分辨率低，宽窗口时间分辨率低、频率分辨率高。对于时变的非稳态信号，高频适合小窗口，低频适合大窗口。然而STFT的窗口是固定的，在一次STFT中宽度不会变化，所以STFT还是无法满足非稳态信号变化的频率的需求。


三、小波变换


那么你可能会想到，让窗口大小变起来，多做几次STFT不就可以了吗？！没错，小波变换就有着这样的思路。

但事实上小波并不是这么做的（关于这一点，方沁园同学的表述“小波变换就是根据算法，加不等长的窗，对每一小部分进行傅里叶变换”就不准确了。小波变换并没有采用窗的思想，更没有做傅里叶变换。）

至于为什么不采用可变窗的STFT呢，我认为是因为这样做冗余会太严重，STFT做不到正交化，这也是它的一大缺陷。


于是小波变换的出发点和STFT还是不同的。STFT是给信号加窗，分段做FFT；而小波直接把傅里叶变换的基给换了——将无限长的三角函数基换成了有限长的会衰减的小波基。这样不仅能够获取频率，还可以定位到时间了~


【解释】

来我们再回顾一下傅里叶变换吧，没弄清傅里叶变换为什么能得到信号各个频率成分的同学也可以再借我的图理解一下。

傅里叶变换把无限长的三角函数作为基函数：


这个基函数会伸缩、会平移（其实是两个正交基的分解）。缩得窄，对应高频；伸得宽，对应低频。然后这个基函数不断和信号做相乘。某一个尺度（宽窄）下乘出来的结果，就可以理解成信号所包含的当前尺度对应频率成分有多少。于是，基函数会在某些尺度下，与信号相乘得到一个很大的值，因为此时二者有一种重合关系。那么我们就知道信号包含多少该频率的成分。

（看，这两种尺度能乘出一个大的值，所以信号包含较多的这两个频率成分，在频谱上这两个频率会出现两个峰）（看，这两种尺度能乘出一个大的值，所以信号包含较多的这两个频率成分，在频谱上这两个频率会出现两个峰）


以上，就是粗浅意义上傅里叶变换的原理。





如前边所说，小波做的改变就在于，将无限长的三角函数基换成了有限长的会衰减的小波基。
这就是为什么它叫“小波”，因为是很小的一个波嘛~这就是为什么它叫“小波”，因为是很小的一个波嘛~



从公式可以看出，不同于傅里叶变换，变量只有频率ω，小波变换有两个变量：尺度a（scale）和平移量 τ（translation）。尺度a控制小波函数的伸缩，平移量τ控制小波函数的平移。尺度就对应于频率（反比），平移量 τ就对应于时间。




当伸缩、平移到这么一种重合情况时，也会相乘得到一个大的值。这时候和傅里叶变换不同的是，这不仅可以知道信号有这样频率的成分，而且知道它在时域上存在的具体位置。


而当我们在每个尺度下都平移着和信号乘过一遍后，我们就知道信号在每个位置都包含哪些频率成分。


看到了吗？有了小波，我们从此再也不害怕非稳定信号啦！从此可以做时频分析啦！


做傅里叶变换只能得到一个频谱，做小波变换却可以得到一个时频谱！
↑：时域信号 ↑：时域信号
↑：傅里叶变换结果↑：傅里叶变换结果

——此图像来源于“THE WAVELET TUTORIAL” ——此图像来源于“THE WAVELET TUTORIAL”

↑：小波变换结果


小波还有一些好处：

1. 我们知道对于突变信号，傅里叶变换存在吉布斯效应，我们用无限长的三角函数怎么也拟合不好突变信号：
然而衰减的小波就不一样了：然而衰减的小波就不一样了：




2. 小波可以实现正交化，短时傅里叶变换不能。



以上，就是小波的意义。


-----------------------------------------------------------------------------------------------------------


以上只是用形象地给大家展示了一下小波的思想，希望能对大家的入门带来一些帮助。毕竟如果对小波一无所知，直接去看那些堆砌公式、照搬论文语言的教材，一定会痛苦不堪。

在这里推荐几篇入门读物，都是以感性介绍为主，易懂但并不深入，对大家初步理解小波会很有帮助。文中有的思路和图也选自于其中：

1. THE WAVELET TUTORIAL （强烈推荐，点击链接：INDEX TO SERIES OF TUTORIALS TO WAVELET TRANSFORM BY ROBI POLIKAR）

2. WAVELETS：SEEING THE FOREST AND THE TREES

3. A Really Friendly Guide to Wavelets

4. Conceptual wavelets


但是真正理解透小波变换，这些还差得很远。比如你至少还要知道有一个“尺度函数”的存在，它是构造“小波函数”的关键，并且是它和小波函数一起才构成了小波多分辨率分析，理解了它才有可能利用小波做一些数字信号处理；你还要理解离散小波变换、正交小波变换、二维小波变换、小波包……这些内容国内教材上讲得也很糟糕，大家就一点一点啃吧~有问题欢迎私信我。水平有限，但一定帮助。


第一次在知乎写这么长的回答，多数图都是用MATLAB和PPT自己画出来的，都是利用实验室搬完砖之余的时间一点点弄的，欢迎分享，如转载还请跟我说一声哈~











-------------------------------------------------------------------------------------------------------------------------

2015.3.26

评论中的一些问题的回答：

1. 关于海森堡不确定性原理

不确定性原理，或者叫测不准原理，最早出自量子力学，意为在微观世界，粒子的位置与动量不可同时被确定。但是这个原理并不局限于量子力学，有很多物理量都有这样的特征，比如能量和时间、角动量和角度。体现在信号领域就是时域和频域。不过更准确一点的表述应该是：一个信号不能在时空域和频域上同时过于集中；一个函数时域越“窄”，它经傅里叶变换的频域后就越“宽”。

如果有兴趣深入研究一下的话，这个原理其实非常耐人寻味。信号处理中的一些新理论在根本上都和它有所相连，比如压缩感知。如果你剥开它复杂的数学描述，最后会发现它在本质上能实现就源于不确定性原理。而且大家不觉得这样一些矛盾的东西在哲学意义上也很奇妙吗，世界观感觉就此被改变了。。



2. 关于正交化

什么是正交化？为什么说小波能实现正交化是优势?

简单说，如果采用正交基，变换域系数会没有冗余信息，等于是用最少的数据表达最大的信息量，利于数值压缩等领域。JPEG2000压缩就是用正交小波变换。

比如典型的正交基：二维笛卡尔坐标系的（1,0）、（0,1），用它们表达一个信号显然非常高效，计算简单。而如果用三个互成120°的向量表达，则会有信息冗余，有重复表达。

但是并不意味着正交一定优于不正交。比如如果是做图像增强，有时候反而希望能有一些冗余信息，更利于对噪声的抑制和对某些特征的增强。


3. 关于瞬时频率

原问题：图中时刻点对应一频率值，一个时刻点只有一个信号值，又怎么能得到他的频率呢？

很好的问题。如文中所说，绝对意义的瞬时频率其实是不存在的。单看一个时刻点的一个信号值，当然得不到它的频率。我们只不过是用很短的一段信号的频率作为该时刻的频率，所以我们得到的只是时间分辨率有限的近似分析结果。这一想法在STFT上体现得很明显。小波等时频分析方法，如用衰减的基函数去测定信号的瞬时频率，思想也类似。


4. 关于小波变换的缺点

这要看和谁比了。

A.作为图像处理方法，和多尺度几何分析方法（超小波）比：

对于图像这种二维信号的话，二维小波变换只能沿2个方向进行，对图像中点的信息表达还可以，但是对线就比较差，这时候ridgelet（脊波）, curvelet（曲波）等多尺度几何分析方法就更有优势了。

B. 作为时频分析方法，和HHT比：

相比于HHT等时频分析方法，小波依然没脱离海森堡测不准原理的束缚，某种尺度下，不能在时间和频率上同时具有很高的精度；以及小波是非适应性的，基函数选定了就不改了。

知识有限，暂时想到的有这些，欢迎补充。

网友方沁园对[短时傅里叶变换]能不能通俗的讲解下傅立叶分析和小波分析之间的关系？给出的答复：
两句话：
傅里叶变换：知道一段时间内，信号的各个频率分量有分别多少。
小波变换：知道一段时间内，信号的各个频率分量有分别多少，以及他们都是什么时候出现的。


解释：

傅里叶变换是把时域的函数比作无数不同频率的余弦函数的叠加，计算每个频率分量有多少，形成频谱。
出现了问题： 如果每个频率分量都在整个时间域上持续，没问题。

但如果有些频率分量可能一开始没有，只在某一段时间范围内出现（非平稳过程）。频谱图上看不出来。

也就是说，时域相差很大的两个信号，可能频谱图一样！

所以，对非平稳过程，傅里叶变换有局限性！！


于是出现了短时傅里叶变换STFT--把整个时域过程分解成无数个等长的小过程，每个小过程近似平稳，再傅里叶变换，就知道在哪个时间点上出现了什么频率了。


其实“等长小过程”还是有局限性。如果根据某些原理，分解成不等长的小过程，就被优化了。--小波变换。就是根据算法，加不等长的窗，对每一小部分进行傅里叶变换。结果就是更更准确的知道在哪个时间点出现什么频率了。

网友蒙面大侠对[短时傅里叶变换]能不能通俗的讲解下傅立叶分析和小波分析之间的关系？给出的答复：

建议参考matlab里的wave toolbox(命令行里输入‘wavemenu’),然后参照Matlab小波工具箱的使用1_了凡春秋和Matlab小波工具箱的使用2_了凡春秋,基本上大概可以弄懂小波分析大概是干什么的,对于具体的理论细节如:小波基的选取,尺度的确定,阀值的选定等, 可以慢慢通过相关的书籍了解.


PS：有傅里叶分析基础理解起来比较快.大概应用的思路就是:把原始信号变换到一个特定的域里(傅里叶变化就是频域，小波变化姑且称之为小波域)，在这个小波域中，可以反映出在原始信号在时域中并不明显的特征，加以阀值或其他相关的处理，再重构回时域.（以除白噪声为例，白噪声在时域和小波域里都是均匀分布的，而有用的纯净信号变换到了小波域不同尺度中的特性不同,利用这个特征设定相应阀值,再重构回时域.可以去除原始信号中的白噪声）

时频结合和多分辨率是小波分析的优势.

网友xinqing moyi对[短时傅里叶变换]能不能通俗的讲解下傅立叶分析和小波分析之间的关系？给出的答复：

谢谢

傅立叶变换 有局限性：只能知道频域上的事情

而小波变换 利用窗函数 既可以知道频域上的事 也可以知道时间域上的事


傅立叶变换 可以让你知道发生了一件事情 它有多厉害

而小波变换 不仅可以让你知道这件事有多厉害 还可以让你知道这件事情什么时候发生的

网友蒙面大侠对[短时傅里叶变换]能不能通俗的讲解下傅立叶分析和小波分析之间的关系？给出的答复：

说到底都是在搞基，只是基不同而已

网友杨硕对[短时傅里叶变换]能不能通俗的讲解下傅立叶分析和小波分析之间的关系？给出的答复：

借这个地方写个笑话，文科生也能看懂。

-----------------------------------------------------------------------------------------------------------------------------------------

我没有认真学过wavelet，因为当时有人给我推荐了一个很好的教材，文章名字是这样的：


看了几页就哭了：我真是连kids都不如啊……

------------------------------------------------------------------------------------------------------------------------------------------


这篇文章链接在此：http://gtwavelet.bme.gatech.edu/wp/kidsA.pdf


这篇文章是Google能找到citation最多的关于wavelet的tutorial，热心网友一致好评：


网友傻妹妹牌炸酱面对[短时傅里叶变换]能不能通俗的讲解下傅立叶分析和小波分析之间的关系？给出的答复：

两者基函数不同。傅里叶是固定基函数，小波则是尺度和平移可变的基函数

网友张生对[短时傅里叶变换]能不能通俗的讲解下傅立叶分析和小波分析之间的关系？给出的答复：

看了本文你还不懂傅里叶变换，那就来掐死我吧，相当好

网友荔枝对[短时傅里叶变换]能不能通俗的讲解下傅立叶分析和小波分析之间的关系？给出的答复：

麻省理工开放课程_线性代数

Lec31_基变换及图像压缩
教授讲的很有意思，也许对你有所帮助教授讲的很有意思，也许对你有所帮助

傅里叶变换、短时傅里叶变换、小波变换
https://www.cnblogs.com/warmbeast/p/7809286.html