Python音频信号处理 1.短时傅里叶变换及其逆变换

短时傅里叶变换及其逆变换

本篇文章主要记录了使用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}}$
$=∑n=0N−1x[n]e−2jπknM\qquad \qquad \qquad \qquad = \sum^{N-1}_{n=0}x[n]e^{-\cfrac{2j\pi kn}{M}}$

离散傅里叶变换适用于在时域上不连续且有限的数字信号，在上述公式中，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}dt$

其中， $h∗(t−τ)\ h^*(t-\tau)$ 是一个中心为 $τ\tau$ 的窗函数

当引入了时间变量 $τ\tau$ 之后我们就可以针对不同瞬间进行频谱分析，对于每一个瞬间 $τ\tau$ 我们都可以获取信号在该时刻的频谱。

二. 使用python 进行短时傅里叶变换

在这一部分我会分享基于python的短时傅里叶变换的实现，可供参考。
首先是可能使用到的python库

import numpy as np
import matplotlib.pyplot as plt
from scipy.io import wavfile
from IPython.display import display, Audio
from 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 /=K
return 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的值与n无关。在此基础上不难证明
$K≈∑n=0Nwin−1w[n]/Nhop\ K \approx \sum^{Nwin-1}_{n=0}w[n] / Nhop$