短时傅里叶变换及其逆变换
本篇文章主要记录了使用python进行短时傅里叶变换,分析频谱,以及通过频谱实现在频域内降低底噪的代码及分析,希望可以给同样在学习信号处理的大家一点帮助,也希望大家对我的文章多提意见建议。
一. 短时傅里叶变换与离散傅里叶变换
在这篇文章中我们主要运用了短时傅里叶变换,要想清楚地理解短时傅里叶变换,首先必须要了解离散傅里叶变换(Discrete Fourier Transform,DFT)。
1.离散傅里叶变换
离散傅里叶的定义:
∀k∈[0,M−1],X[k]=∑n=0N−1x[n]e−2jπnfkFe=∑n=0N−1x[n]e−2jπnkFeMFe\\{\forall} k\in [0,M-1] ,X[k]=\sum^{N-1}_{n=0}x[n] e^{-\cfrac{2j\pi nf_k} {F_e}} = \sum^{N-1}_{n=0}x[n]e^{-\cfrac{2j\pi nk\frac{F_e}{M}}{Fe}} ∀k∈[0,M−1],X[k]=n=0∑N−1x[n]e−Fe2jπnfk=n=0∑N−1x[n]e−Fe2jπnkMFe
=∑n=0N−1x[n]e−2jπknM\qquad \qquad \qquad \qquad = \sum^{N-1}_{n=0}x[n]e^{-\cfrac{2j\pi kn}{M}}=∑n=0N−1x[n]e−M2jπkn
离散傅里叶变换适用于在时域上不连续且有限的数字信号,在上述公式中,x[n] 就是我们在时域中的初始数字信号,Fe 对应这个信号的采样频率。在离散傅里叶变换中,首先初始数字信号本身是离散的,在上式中初始信号x[n]是在时域内的一段有限信号,N代表了该段数字信号一共包含N个采样点,即其在时域上的长度为1/Fe * N。同时离散傅里叶变换所得的结果X[k]在频域上也是离散的,简而言之,是将频域[0,Fe]等分成了M份 :
X[k]也可以理解为包含了M个复数值的向量。在离散傅里叶变换中,采样点的个数N(时域上的取样长度),以及频域上的采样个数M都是可调整的,两个采样个数的选择对于DFT的解析度和精确度会有影响,这里就不过多展开。
我们在实际应用离散傅里叶变换时,会发现,有的时域上完全不同的信号,他们的离散傅里叶变换频谱却是一致的,因此我们引入一个新的 短时傅里叶变换(STFT)。
2.短时傅里叶变换
短时傅里叶变换的定义 :
X(τ,f)=∫Rx(t)h∗(t−τ)e−2jπftdtX(\tau,f)=\int_Rx(t)h^*(t-\tau)e^{-2j\pi ft}dtX(τ,f)=∫Rx(t)h∗(t−τ)e−2jπftdt
其中,h∗(t−τ)\ h^*(t-\tau)h∗(t−τ) 是一个中心为τ\tauτ的窗函数
当引入了时间变量τ\tauτ之后我们就可以针对不同瞬间进行频谱分析,对于每一个瞬间τ\tauτ我们都可以获取信号在该时刻的频谱。
二. 使用python 进行短时傅里叶变换
在这一部分我会分享基于python的短时傅里叶变换的实现,可供参考。
首先是可能使用到的python库
import numpy as npimport matplotlib.pyplot as pltfrom scipy.io import wavfilefrom IPython.display import display, Audiofrom numpy import log10
接下来就是短时傅里叶变换的python实现
def TFCT(trame, Fe, Nfft,fenetre,Nwin,Nhop):L = round((len(trame) - len(fenetre))/Nhop)+1M = Nfftxmat = np.zeros((M,L))print('xmat',xmat.shape)print(Nwin+Nhop)for j in range(L):xmat[:,j] = np.fft.fft(trame[j*Nhop:Nwin+Nhop*j]*fenetre,Nfft) x_temporel = np.linspace(0,(1/Fe)*len(trame),len(trame))x_frequentiel = np.linspace(0, Fe,Nfft)return xmat,x_temporel,x_frequentiel
上述函数解释:
参数部分
trame和Fe: 初始的数字信号和它的采样频率
Nfft: 上文提到的离散傅里叶变换中,频域的采样个数M
fenetre: 短时傅里叶变换中使用的窗函数,在接下来的实现中我都使用了汉明窗np.hamming。
Nwin: 窗函数的长度(包含点的个数)
Nhop: 窗函数每次滑动的长度,一般规定为Nwin/2,窗函数长度的一半
首先创建一个M行L列的矩阵xmat,该矩阵的每一行代表一个0-Fe的频率,单位为Hz,每一列对应该段被窗函数截取的信号的FFT快速傅里叶变换。
三. 使用overlapp-add算法进行短时傅里叶变换的逆变换重构原信号
在这一部分中,我们使用了overlapp-add算法来进行短时傅里叶变换的逆变换。
下面是该部分的全部代码,之后会逐步解释算法的实现 :
def ITFD(xmat,Fe,Nfft,Nwin,Nhop):window = np.hamming(Nwin)Te = 1/Feyvect = np.zeros(Nfft + (xmat.shape[1]-1)*Nhop,dtype=complex)t_vecteur = np.arange(0,Te*len(yvect),Te)index = 0K = 0L = xmat.shape[1]yl = np.zeros(xmat.shape,dtype=complex)for j in range(L):yl[:,j] = np.fft.ifft(xmat[:,j])# 平移和求和for k in range(L):yvect[Nhop*k:Nfft+Nhop*k] += yl[:,k]# 标准化幅值for n in range(Nwin-1):K += window[n]K /= Nhopyvect /=Kreturn t_vecteur, yvect
该算法的实现需要三步。
1. 快速傅里叶逆变换
yl = np.zeros(xmat.shape,dtype=complex)for j in range(L):yl[:,j] = np.fft.ifft(xmat[:,j])
第一步,对上一部分得出的矩阵xmat进行快速傅里叶变换的逆变换,得出同样规格M行L列的矩阵yl。
2. 对各列进行平移并叠加
# 平移和求和for k in range(L):yvect[Nhop*k:Nfft+Nhop*k] += yl[:,k]
对yl矩阵的每一列平移 (l-1)Nhop,l ∈\in∈ [1,L],例如第一列不变,第二列平移Nhop,第三列平移2Nhop,以此类推。之后将所有列的转置,叠加到总长度为Nfft +(L-1)*Nhop的向量yvect中。
3. 标准化
# 标准化幅值for n in range(Nwin-1):K += window[n]K /= Nhopyvect /=Kreturn t_vecteur, yvect
window[n] (w[n]) 是长度为Nwin的窗函数,在选取窗函数的时候,我们总满足规则 K=∑l=1Lw[n−(l−1)Nhop]\ K=\sum^L_{l=1}w[n-(l-1)Nhop]K=l=1∑Lw[n−(l−1)Nhop] K的值与n无关。在此基础上不难证明
K≈∑n=0Nwin−1w[n]/Nhop\ K \approx \sum^{Nwin-1}_{n=0}w[n] / NhopK≈n=0∑Nwin−1w[n]/Nhop
那么通过以上的三个步骤,我们就可以从信号的短时傅里叶变换矩阵中完美重构原信号了。
下一篇文章我们将使用这些算法,使用谱减法进行声音信号的降噪处理。