深入理解短时傅里叶变换 STFT + Python 代码详解
最近在Coursera上面听课，写博客记录一下，促进学习。Audio Signal Processing for Music ApplicationsPS: 这里再推一波郭宝龙老师的MOOC，很香。信号与线性系统分析 吴大正 郭宝龙(1) STFT 的引入1. 为什么要引入短时傅里叶变换?首先引入一个平稳与非平稳的概念：非平稳信号：对于确定信号，信号的频率成分随着时间而发生变化。比如说下面的脑电图
- 最近在Coursera上面听课，写博客记录一下，促进学习。
Audio Signal Processing for Music Applications - PS: 这里再推一波郭宝龙老师的MOOC，很香。
信号与线性系统分析 吴大正 郭宝龙
(1) STFT 的引入
1. 为什么要引入短时傅里叶变换?
-
首先引入一个平稳与非平稳的概念：
- 非平稳信号：对于确定信号，信号的频率成分随着时间而发生变化。比如说下面的脑电图信号:
- 平稳信号：对于确定信号，信号的频率成分与时间无关。比如说一个单一的正弦函数信号。
-
非平稳信号分析的挑战：
这里使用一个例子来说明，例如对以下两个函数(一个平稳、另一个非平稳) f 1 ( t ) = cos ( 2 π × 10 t ) + cos ( 2 π × 25 t ) + cos ( 2 π × 50 t ) + cos ( 2 π × 100 t ) f_1(t) = \cos(2\pi\times10t) + \cos(2\pi\times25t) + \cos(2\pi\times50t)+\cos(2\pi\times100t) f1(t)=cos(2π×10t)+cos(2π×25t)+cos(2π×50t)+cos(2π×100t)
f 2 ( t ) = { cos ( 2 π × 10 t ) 0 ≤ t ≤ 0.20 cos ( 2 π × 25 t ) 0.20 ≤ t ≤ 0.40 cos ( 2 π × 50 t ) 0.40 ≤ t ≤ 0.70 cos ( 2 π × 100 t ) 0.70 ≤ t ≤ 1.00 f_2(t) = \left\{ \begin{aligned} &\cos(2\pi\times10t)~~~0\le t\le 0.20 \\ &\cos(2\pi\times25t)~~~0.20\le t\le 0.40 \\ &\cos(2\pi\times50t)~~~0.40\le t\le 0.70 \\ &\cos(2\pi\times100t)~~~0.70\le t\le 1.00 \\ \end{aligned}\right. f2(t)=⎩⎪⎪⎪⎪⎨⎪⎪⎪⎪⎧cos(2π×10t) 0≤t≤0.20cos(2π×25t) 0.20≤t≤0.40cos(2π×50t) 0.40≤t≤0.70cos(2π×100t) 0.70≤t≤1.00
两个信号的幅频特性：四个主要的尖峰。
50Hz和100Hz分量的幅度比25Hz和10Hz分量大，这是因为高频信号比低频信号持续时间更长一些（分别为300ms和200ms）。若忽略掉因频率突变引起的毛刺（有时候他们与噪声很难区分）和两幅图中各频率分量的幅值（这些幅值可以做归一化处理），两个信号的频谱图几乎是一致的，但实际上两个时域信号的差别极大。 -
说明引入的原因：
（1）傅里叶变换的全局积分导致变换结果无法提供频率分量的时间信息。
（2）对于非平稳信号来说，傅里叶变换一般是不合适的。
（3）只有仅仅关心信号中是否包含某个频率分量而不关心它出现的时间的时候，傅里叶变换才可以用于处理非平稳信号。
概括起来就是：增加一维的自由度，使频谱可以反映其不同频率状态所处的不同时间。
2. 短时傅里叶变换 (离散的方程)
- 从时域到频域的表达式（正变换）
X l [ k ] = ∑ n = − N / 2 N / 2 − 1 w [ n ] x [ n + l H ] e − j 2 π k n / N l = 0 , 1 , . . . , L X_l[k] = \sum_{n=-N/2}^{N/2-1}w[n]x[n+lH]e^{-j2\pi kn/N}~~~l=0,1,...,L Xl[k]=n=−N/2∑N/2−1w[n]x[n+lH]e−j2πkn/N l=0,1,...,L
实质就是在时域加窗，再对加窗后的信号进行频域的分析。 w w w 表示分析窗，长度为 N， l l l 表示帧序号， H H H 表示不同帧之间的跳跃步长。 - 从频域到时域的表达式（反变换）
y [ n ] = ∑ l = 0 L − 1 S h i f t l H , n [ 1 N ∑ k = − N / 2 N / 2 − 1 X l [ k ] e j 2 π k n / N ] y[n] = \sum_{l=0}^{L-1}Shift_{lH,n}[\frac{1}{N}\sum_{k=-N/2}^{N/2-1}X_l[k] e^{j2\pi kn/N}] y[n]=l=0∑L−1ShiftlH,n[N1k=−N/2∑N/2−1Xl[k]ej2πkn/N]
实际上反变换就是首先对每一帧进行FFT的反变换得到加窗后的 x [ n ] x[n] x[n]，然后再对这些 x [ n ] x[n] x[n] 进行时移相加就可以得到还原后的原信号。
- 流程图
(2) 分析窗的说明
1. 使用中需要关注的值
- main lobe：主瓣的宽度。
- highest side lobe：最高旁瓣。
2. 列举一些常用的分析窗
-
矩形窗
- w [ n ] = 1 , n = − M / 2 , . . . , 0 , . . . , M / 2 w[n] = 1, n = -M/2,...,0,...,M/2 w[n]=1,n=−M/2,...,0,...,M/2。
-
W
[
k
]
=
sin
(
π
k
)
sin
(
π
k
/
M
)
W[k] = \displaystyle\frac{\sin(\pi k)}{\sin (\pi k /M)}
W[k]=sin(πk/M)sin(πk)。
~
取 M = 64 M=64 M=64 (即窗宽 65)，FFT长度 N = 1024 N=1024 N=1024。从频域可以看到，主瓣宽度 2 格， 旁瓣最高值 -13.3 dB。
- hanning 窗
- w [ n ] = 0.5 + 0.5 cos ( 2 π n / M ) , n = − M / 2 , . . . , 0 , . . . , M / 2 w[n] = 0.5+0.5\cos(2\pi n/M), n = -M/2,...,0,...,M/2 w[n]=0.5+0.5cos(2πn/M),n=−M/2,...,0,...,M/2。
-
W
[
k
]
=
0.5
D
[
k
]
+
0.25
(
D
[
k
−
1
]
+
D
[
k
+
1
]
)
W[k] = 0.5D[k] + 0.25(D[k-1]+D[k+1])
W[k]=0.5D[k]+0.25(D[k−1]+D[k+1])，其中
D
[
k
]
=
sin
(
π
k
)
sin
(
π
k
/
M
)
D[k] = \displaystyle\frac{\sin(\pi k)}{\sin (\pi k /M)}
D[k]=sin(πk/M)sin(πk)。
~
取 M = 64 M=64 M=64 (即窗宽 65)，FFT长度 N = 512 N=512 N=512。从频域可以看到，主瓣宽度 4 格， 旁瓣最高值 -31.5 dB。
- 汉明窗
- w [ n ] = 0.54 + 0.46 cos ( 2 π n / M ) , n = − M / 2 , . . . , 0 , . . . , M / 2 w[n] = 0.54+0.46\cos(2\pi n/M), n = -M/2,...,0,...,M/2 w[n]=0.54+0.46cos(2πn/M),n=−M/2,...,0,...,M/2。
-
W
[
k
]
=
0.54
D
[
k
]
+
0.23
(
D
[
k
−
1
]
+
D
[
k
+
1
]
)
W[k] = 0.54D[k] + 0.23(D[k-1]+D[k+1])
W[k]=0.54D[k]+0.23(D[k−1]+D[k+1])，其中
D
[
k
]
=
sin
(
π
k
)
sin
(
π
k
/
M
)
D[k] = \displaystyle\frac{\sin(\pi k)}{\sin (\pi k /M)}
D[k]=sin(πk/M)sin(πk)。
~
取 M = 64 M=64 M=64 (即窗宽 65)，FFT长度 N = 512 N=512 N=512。从频域可以看到，主瓣宽度 4 格， 旁瓣最高值 -42.7 dB。
- blackman 窗
-
w
[
n
]
=
0.42
−
0.5
cos
(
2
π
n
/
M
)
+
0.08
cos
(
4
π
n
/
M
)
,
n
=
−
M
/
2
,
.
.
.
,
0
,
.
.
.
,
M
/
2
w[n] = 0.42-0.5\cos(2\pi n/M)+0.08\cos(4\pi n/M), n = -M/2,...,0,...,M/2
w[n]=0.42−0.5cos(2πn/M)+0.08cos(4πn/M),n=−M/2,...,0,...,M/2。
-
W
[
k
]
=
0.42
D
[
k
]
−
0.25
(
D
[
k
−
1
]
+
D
[
k
+
1
]
)
+
0.04
(
D
[
k
−
2
]
+
D
[
k
+
2
]
)
其
中
D
[
k
]
=
sin
(
π
k
)
sin
(
π
k
/
M
)
W[k] = 0.42D[k] - 0.25(D[k-1]+D[k+1]) + 0.04(D[k-2]+D[k+2])其中 D[k] = \displaystyle\frac{\sin(\pi k)}{\sin (\pi k /M)}
W[k]=0.42D[k]−0.25(D[k−1]+D[k+1])+0.04(D[k−2]+D[k+2])其中D[k]=sin(πk/M)sin(πk)。
~
取 M = 64 M=64 M=64 (即窗宽 65)，FFT长度 N = 512 N=512 N=512。从频域可以看到，主瓣宽度 6 格， 旁瓣最高值 -58 dB。
-
w
[
n
]
=
0.42
−
0.5
cos
(
2
π
n
/
M
)
+
0.08
cos
(
4
π
n
/
M
)
,
n
=
−
M
/
2
,
.
.
.
,
0
,
.
.
.
,
M
/
2
w[n] = 0.42-0.5\cos(2\pi n/M)+0.08\cos(4\pi n/M), n = -M/2,...,0,...,M/2
w[n]=0.42−0.5cos(2πn/M)+0.08cos(4πn/M),n=−M/2,...,0,...,M/2。
- blackman-harris 窗
-
w
[
n
]
=
1
M
∑
l
=
0
3
α
l
cos
(
2
n
l
π
/
M
)
,
n
=
−
M
/
2
,
.
.
.
,
0
,
.
.
.
,
M
/
2
w[n] =\displaystyle\frac{1}{M}\sum_{l=0}^{3}\alpha_l\cos(2nl\pi/M), n = -M/2,...,0,...,M/2
w[n]=M1l=0∑3αlcos(2nlπ/M),n=−M/2,...,0,...,M/2，其中
α
0
=
0.35875
,
α
1
=
0.48829
,
α
2
=
0.14128
,
α
3
=
0.01168
\alpha_0=0.35875,\alpha_1=0.48829,\alpha_2=0.14128,\alpha_3=0.01168
α0=0.35875,α1=0.48829,α2=0.14128,α3=0.01168
取 M = 64 M=64 M=64 (即窗宽 65)，FFT长度 N = 512 N=512 N=512。从频域可以看到，主瓣宽度 8 格， 旁瓣最高值 -92 dB。
-
w
[
n
]
=
1
M
∑
l
=
0
3
α
l
cos
(
2
n
l
π
/
M
)
,
n
=
−
M
/
2
,
.
.
.
,
0
,
.
.
.
,
M
/
2
w[n] =\displaystyle\frac{1}{M}\sum_{l=0}^{3}\alpha_l\cos(2nl\pi/M), n = -M/2,...,0,...,M/2
w[n]=M1l=0∑3αlcos(2nlπ/M),n=−M/2,...,0,...,M/2，其中
α
0
=
0.35875
,
α
1
=
0.48829
,
α
2
=
0.14128
,
α
3
=
0.01168
\alpha_0=0.35875,\alpha_1=0.48829,\alpha_2=0.14128,\alpha_3=0.01168
α0=0.35875,α1=0.48829,α2=0.14128,α3=0.01168
(3) 不同参数对频域的影响
1. 不同类型的窗对频域的影响
使用不同的窗，频域曲线的平滑程度，不同频率分量信息的清晰程度均有区别。从下面这个例子中可以看到 Blackman Window 的效果相对较好，频域曲线相对比较平滑并且不同频率分量的信息也比较清晰。
2. 窗的长度对频域的影响
不同的窗长影响频域幅度分量信息的含量，窗越长，频域所反映出来的信息也越多。
3. 窗的长度的奇偶性对频域的影响
窗长度的奇偶性主要影响频域的相位分量的对称性，窗长为奇数时（对应M为偶数）相位是关于纵坐标轴完全对称的，因此一般我们选择窗长为奇数（对应M为偶数）。
4. FFT的长度对频域的影响
更长的 FFT 的长度（使用到Zero Padding）可以保证频域幅度曲线更加平滑，并且幅度的细节体现的更好。
5. 窗的跳跃长度对频域的影响
我们希望的情况是通过叠加，最终呈现结果是一条取值为 1 的直线，这样可以保证不会由叠加这个行为产生信号畸变与失真。
(4) 时域频域分辨率之间的关系
窗函数的宽度与信号的表示方式有关，它决定是具有良好的频率分辨率(靠的很近的频率成分可以被分离)还是良好的时间分辨率(频率变化的时间)。较宽的窗口提供更好的频率分辨率，但时间分辨率较差。较窄的窗口给出较好的时间分辨率，但频率分辨率较差。这个性质与海森堡测不准原理有关，但更具体的可以参考Gabor极限。
可以用一个例子来更好的说明这个问题，对如下这个函数采用不同的窗长进行变换：
f ( t ) = { cos ( 2 π × 10 t ) 0 ≤ t ≤ 5 cos ( 2 π × 25 t ) 5 ≤ t ≤ 10 cos ( 2 π × 50 t ) 10 ≤ t ≤ 15 cos ( 2 π × 100 t ) 15 ≤ t ≤ 20 f(t) = \left\{ \begin{aligned} &\cos(2\pi\times10t)~~~0\le t\le 5 \\ &\cos(2\pi\times25t)~~~5\le t\le 10 \\ &\cos(2\pi\times50t)~~~10\le t\le 15 \\ &\cos(2\pi\times100t)~~~15\le t\le 20 \\ \end{aligned}\right. f(t)=⎩⎪⎪⎪⎪⎨⎪⎪⎪⎪⎧cos(2π×10t) 0≤t≤5cos(2π×25t) 5≤t≤10cos(2π×50t) 10≤t≤15cos(2π×100t) 15≤t≤20
25 ms 的窗保证了比较好的时间分辨率，但是频率分辨率很差；1000 ms 的窗保证了比较好的频率分辨率，但是时间分辨率很差。
(5) 代码实现（python）
- 有几个点需要注意：
stftAnal函数输出结果的单位是 dB。stftAnal函数的处理中使用了在信号 x 的最前面与最后面补零这一技巧，目的是为了尽可能的使输出和输入的信号保持一致，减小因为进行stft而引入的噪声。DFT.dftAnal与DFT.dftSynth下面直接给出代码（博主GPS还没预习完，先偷个懒直接附上代码）。
- 从时域到频域
def stftAnal(x, w, N, H): """ x: 输入信号, w: 分析窗, N: FFT 的大小, H: hop 的大小 返回 xmX, xpX: 振幅和相位，以 dB 为单位 """ M = w.size # 分析窗的大小 hM1 = (M+1)//2 hM2 = M//2 x = np.append(np.zeros(hM2),x) # 在信号 x 的最前面与最后面补零 x = np.append(x,np.zeros(hM2)) pin = hM1 # 初始化指针，用来指示现在指示现在正在处理哪一帧 pend = x.size-hM1 # 最后一帧的位置 w = w / sum(w) # 归一化分析窗 xmX = [] xpX = [] while pin<=pend: x1 = x[pin-hM1:pin+hM2] # 选择一帧输入的信号 mX, pX = dftAnal(x1, w, N) # 计算 DFT（这个函数不是库中的） xmX.append(np.array(mX)) # 添加到 list 中 xpX.append(np.array(pX)) pin += H # 更新指针指示的位置 xmX = np.array(xmX) # 转换为 numpy 数组 xpX = np.array(xpX) return xmX, xpX - 从频域到时域
def stftSynth(mY, pY, M, H) : """ mY: 振幅谱以dB为单位, pY: 相位谱, M: 分析窗的大小, H: hop 的大小 返回 y 还原后的信号 """ hM1 = (M+1)//2 hM2 = M//2 nFrames = mY[:,0].size # 计算帧的数量 y = np.zeros(nFrames*H + hM1 + hM2) # 初始化输出向量 pin = hM1 for i in range(nFrames): # 迭代所有帧 y1 = dftSynth(mY[i,:], pY[i,:], M) # 计算IDFT（这个函数不是库中的） y[pin-hM1:pin+hM2] += H*y1 # overlap-add pin += H # pin是一个指针，用来指示现在指示现在正在处理哪一帧 y = np.delete(y, range(hM2)) # 删除头部在stftAnal中添加的部分 y = np.delete(y, range(y.size-hM1, y.size)) # 删除尾部在stftAnal中添加的部分 return y
附上dft的两段代码：
def dftAnal(x, w, N):
"""
Analysis of a signal using the discrete Fourier transform
x: input signal, w: analysis window, N: FFT size
returns mX, pX: magnitude and phase spectrum
"""
if not(UF.isPower2(N)): # raise error if N not a power of two
raise ValueError("FFT size (N) is not a power of 2")
if (w.size > N): # raise error if window size bigger than fft size
raise ValueError("Window size (M) is bigger than FFT size")
hN = (N//2)+1 # size of positive spectrum, it includes sample 0
hM1 = (w.size+1)//2 # half analysis window size by rounding
hM2 = w.size//2 # half analysis window size by floor
fftbuffer = np.zeros(N) # initialize buffer for FFT
w = w / sum(w) # normalize analysis window
xw = x*w # window the input sound
fftbuffer[:hM1] = xw[hM2:] # zero-phase window in fftbuffer
fftbuffer[-hM2:] = xw[:hM2]
X = fft(fftbuffer) # compute FFT
absX = abs(X[:hN]) # compute ansolute value of positive side
absX[absX<np.finfo(float).eps] = np.finfo(float).eps # if zeros add epsilon to handle log
mX = 20 * np.log10(absX) # magnitude spectrum of positive frequencies in dB
X[:hN].real[np.abs(X[:hN].real) < tol] = 0.0 # for phase calculation set to 0 the small values
X[:hN].imag[np.abs(X[:hN].imag) < tol] = 0.0 # for phase calculation set to 0 the small values
pX = np.unwrap(np.angle(X[:hN])) # unwrapped phase spectrum of positive frequencies
return mX, pX
def dftSynth(mX, pX, M):
"""
Synthesis of a signal using the discrete Fourier transform
mX: magnitude spectrum, pX: phase spectrum, M: window size
returns y: output signal
"""
hN = mX.size # size of positive spectrum, it includes sample 0
N = (hN-1)*2 # FFT size
if not(UF.isPower2(N)): # raise error if N not a power of two, thus mX is wrong
raise ValueError("size of mX is not (N/2)+1")
hM1 = int(math.floor((M+1)/2)) # half analysis window size by rounding
hM2 = int(math.floor(M/2)) # half analysis window size by floor
fftbuffer = np.zeros(N) # initialize buffer for FFT
y = np.zeros(M) # initialize output array
Y = np.zeros(N, dtype = complex) # clean output spectrum
Y[:hN] = 10**(mX/20) * np.exp(1j*pX) # generate positive frequencies
Y[hN:] = 10**(mX[-2:0:-1]/20) * np.exp(-1j*pX[-2:0:-1]) # generate negative frequencies
fftbuffer = np.real(ifft(Y)) # compute inverse FFT
y[:hM2] = fftbuffer[-hM2:] # undo zero-phase window
y[hM2:] = fftbuffer[:hM1]
return y
为开发者提供学习成长、分享交流、生态实践、资源工具等服务，帮助开发者快速成长。
更多推荐
39
0- 0
已为社区贡献1条内容
所有评论(0)