最近在Coursera上面听课，写博客记录一下，促进学习。

Audio Signal Processing for Music ApplicationsPS: 这里再推一波郭宝龙老师的MOOC，很香。

信号与线性系统分析 吴大正 郭宝龙

(1) STFT 的引入

1. 为什么要引入短时傅里叶变换?

首先引入一个平稳与非平稳的概念：

非平稳信号：对于确定信号，信号的频率成分随着时间而发生变化。比如说下面的脑电图信号: 平稳信号：对于确定信号，信号的频率成分与时间无关。比如说一个单一的正弦函数信号。

非平稳信号分析的挑战：

这里使用一个例子来说明，例如对以下两个函数(一个平稳、另一个非平稳) f1(t)=cos⁡(2π×10t)+cos⁡(2π×25t)+cos⁡(2π×50t)+cos⁡(2π×100t)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)

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.00f_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. 短时傅里叶变换 (离散的方程)

从时域到频域的表达式（正变换）

Xl[k]=∑n=−N/2N/2−1w[n]x[n+lH]e−j2πkn/Nl=0,1,...,LX_l[k] = \sum_{n=-N/2}^{N/2-1}w[n]x[n+lH]e^{-j2\pi kn/N}~~~l=0,1,...,LXl​[k]=n=−N/2∑N/2−1​w[n]x[n+lH]e−j2πkn/Nl=0,1,...,L

实质就是在时域加窗，再对加窗后的信号进行频域的分析。www 表示分析窗，长度为 N，lll 表示帧序号，HHH 表示不同帧之间的跳跃步长。从频域到时域的表达式（反变换）

y[n]=∑l=0L−1ShiftlH,n[1N∑k=−N/2N/2−1Xl[k]ej2πkn/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−1​ShiftlH,n​[N1​k=−N/2∑N/2−1​Xl​[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/2w[n] = 1, n = -M/2,...,0,...,M/2w[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=64M=64M=64 (即窗宽 65)，FFT长度 N=1024N=1024N=1024。从频域可以看到，主瓣宽度 2 格， 旁瓣最高值 -13.3 dB。hanning 窗 w[n]=0.5+0.5cos⁡(2πn/M),n=−M/2,...,0,...,M/2w[n] = 0.5+0.5\cos(2\pi n/M), n = -M/2,...,0,...,M/2w[n]=0.5+0.5cos(2πn/M),n=−M/2,...,0,...,M/2。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])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=64M=64M=64 (即窗宽 65)，FFT长度 N=512N=512N=512。从频域可以看到，主瓣宽度 4 格， 旁瓣最高值 -31.5 dB。汉明窗 w[n]=0.54+0.46cos⁡(2πn/M),n=−M/2,...,0,...,M/2w[n] = 0.54+0.46\cos(2\pi n/M), n = -M/2,...,0,...,M/2w[n]=0.54+0.46cos(2πn/M),n=−M/2,...,0,...,M/2。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])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=64M=64M=64 (即窗宽 65)，FFT长度 N=512N=512N=512。从频域可以看到，主瓣宽度 4 格， 旁瓣最高值 -42.7 dB。blackman 窗 w[n]=0.42−0.5cos⁡(2πn/M)+0.08cos⁡(4πn/M),n=−M/2,...,0,...,M/2w[n] = 0.42-0.5\cos(2\pi n/M)+0.08\cos(4\pi n/M), n = -M/2,...,0,...,M/2w[n]=0.42−0.5cos(2πn/M)+0.08cos(4πn/M),n=−M/2,...,0,...,M/2。

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)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=64M=64M=64 (即窗宽 65)，FFT长度 N=512N=512N=512。从频域可以看到，主瓣宽度 6 格， 旁瓣最高值 -58 dB。blackman-harris 窗 w[n]=1M∑l=03αlcos⁡(2nlπ/M),n=−M/2,...,0,...,M/2w[n] =\displaystyle\frac{1}{M}\sum_{l=0}^{3}\alpha_l\cos(2nl\pi/M), n = -M/2,...,0,...,M/2w[n]=M1​l=0∑3​αl​cos(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=64M=64M=64 (即窗宽 65)，FFT长度 N=512N=512N=512。从频域可以看到，主瓣宽度 8 格， 旁瓣最高值 -92 dB。

(3) 不同参数对频域的影响

1. 不同类型的窗对频域的影响

使用不同的窗，频域曲线的平滑程度，不同频率分量信息的清晰程度均有区别。从下面这个例子中可以看到 Blackman Window 的效果相对较好，频域曲线相对比较平滑并且不同频率分量的信息也比较清晰。

2. 窗的长度对频域的影响

不同的窗长影响频域幅度分量信息的含量，窗越长，频域所反映出来的信息也越多。

3. 窗的长度的奇偶性对频域的影响

窗长度的奇偶性主要影响频域的相位分量的对称性，窗长为奇数时（对应M为偶数）相位是关于纵坐标轴完全对称的，因此一般我们选择窗长为奇数（对应M为偶数）。

4. FFT的长度对频域的影响

更长的 FFT 的长度（使用到Zero Padding）可以保证频域幅度曲线更加平滑，并且幅度的细节体现的更好。

5. 窗的跳跃长度对频域的影响

我们希望的情况是通过叠加，最终呈现结果是一条取值为 1 的直线，这样可以保证不会由叠加这个行为产生信号畸变与失真。

(4) 时域频域分辨率之间的关系

窗函数的宽度与信号的表示方式有关，它决定是具有良好的频率分辨率(靠的很近的频率成分可以被分离)还是良好的时间分辨率(频率变化的时间)。较宽的窗口提供更好的频率分辨率，但时间分辨率较差。较窄的窗口给出较好的时间分辨率，但频率分辨率较差。这个性质与海森堡测不准原理有关，但更具体的可以参考Gabor极限。

可以用一个例子来更好的说明这个问题，对如下这个函数采用不同的窗长进行变换：

f(t)={cos⁡(2π×10t)0≤t≤5cos⁡(2π×25t)5≤t≤10cos⁡(2π×50t)10≤t≤15cos⁡(2π×100t)15≤t≤20f(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)//2hM2 = 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-addpin += 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 transformx: 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 tworaise ValueError("FFT size (N) is not a power of 2")if (w.size > N):# raise error if window size bigger than fft sizeraise ValueError("Window size (M) is bigger than FFT size")hN = (N//2)+1 # size of positive spectrum, it includes sample 0hM1 = (w.size+1)//2 # half analysis window size by roundinghM2 = w.size//2 # half analysis window size by floorfftbuffer = np.zeros(N) # initialize buffer for FFTw = w / sum(w) # normalize analysis windowxw = x*w # window the input soundfftbuffer[:hM1] = xw[hM2:]# zero-phase window in fftbufferfftbuffer[-hM2:] = xw[:hM2] X = fft(fftbuffer) # compute FFTabsX = abs(X[:hN]) # compute ansolute value of positive sideabsX[absX<np.finfo(float).eps] = np.finfo(float).eps # if zeros add epsilon to handle logmX = 20 * np.log10(absX) # magnitude spectrum of positive frequencies in dBX[:hN].real[np.abs(X[:hN].real) < tol] = 0.0 # for phase calculation set to 0 the small valuesX[: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 frequenciesreturn mX, pXdef dftSynth(mX, pX, M):"""Synthesis of a signal using the discrete Fourier transformmX: magnitude spectrum, pX: phase spectrum, M: window sizereturns y: output signal"""hN = mX.size # size of positive spectrum, it includes sample 0N = (hN-1)*2 # FFT sizeif not(UF.isPower2(N)): # raise error if N not a power of two, thus mX is wrongraise ValueError("size of mX is not (N/2)+1")hM1 = int(math.floor((M+1)/2))# half analysis window size by roundinghM2 = int(math.floor(M/2))# half analysis window size by floorfftbuffer = np.zeros(N) # initialize buffer for FFTy = np.zeros(M) # initialize output arrayY = np.zeros(N, dtype = complex) # clean output spectrumY[:hN] = 10**(mX/20) * np.exp(1j*pX)# generate positive frequenciesY[hN:] = 10**(mX[-2:0:-1]/20) * np.exp(-1j*pX[-2:0:-1]) # generate negative frequenciesfftbuffer = np.real(ifft(Y)) # compute inverse FFTy[:hM2] = fftbuffer[-hM2:]# undo zero-phase windowy[hM2:] = fftbuffer[:hM1]return y