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

本篇文章主要记录了使用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 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=l=1∑Lw[n−(l−1)Nhop] K的值与n无关。在此基础上不难证明
K≈∑n=0Nwin−1w[n]/Nhop\ K \approx \sum^{Nwin-1}_{n=0}w[n] / Nhop K≈n=0∑Nwin−1w[n]/Nhop

那么通过以上的三个步骤，我们就可以从信号的短时傅里叶变换矩阵中完美重构原信号了。
下一篇文章我们将使用这些算法，使用谱减法进行声音信号的降噪处理。

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

Python音频信号处理 2.使用谱减法去除音频底噪
使用谱减法去除音频底噪上一篇文章我主要分享了短时傅立叶变换及其逆变换在python中的实现,有兴趣的可以阅读一下该篇文章,地址如下: Python音频信号处理 1.短时傅里叶变换及其逆变换那么在本 ...
Python音频信号处理库函数librosa介绍
文章目录 Python音频信号处理库函数librosa介绍(部分内容将陆续添加) 介绍安装综述(库函数结构) Core IO and DSP(核心输入输出功能和数字信号处理) Audio proc ...
【基于pyAudioKits的Python音频信号处理（一）】pyAudioKits安装与API速查手册
文章目录 pyAudioKits 基本用法创建或加载音频来自NumPy数组来自文件录音模拟 Audio对象播放绘制转为NumPy数组获取属性保存索引和切片连接合成四则运算 ...
【基于pyAudioKits的Python音频信号处理（七）】端点检测和语音识别
pyAudioKits是基于librosa和其他库的强大Python音频工作流支持. API速查手册通过pip安装: pip install pyAudioKits 本项目的GitHub地址,如果这 ...
【音频处理】短时傅里叶变换
前言上一篇博客讲了离散傅里叶变换,里面的实例是对整个信号进行计算,虽然理论上有N点傅里叶变换(本博客就不区分FFT和DFT了,因为它俩就是一个东东,只不过复杂度不同),但是我个人理解是这个N点是信号 ...
【基于pyAudioKits的Python音频信号处理项目（二）】深度学习语音识别
pyAudioKits是基于librosa和其他库的强大Python音频工作流支持. API速查手册通过pip安装: pip install pyAudioKits 本项目的GitHub地址,如果这 ...
python实时采集与处理声音信号_python之音频信号处理一
python音频信号处理,首先安装librosa模块安装好librosa模块后,进行简单的音频读取操作,包括: 1.load读取音频文件,返回音频数据与采样率 path:音频文件路径 ...
Python 音频识别以及信号处理（一）
Python 音频信号处理(一) 音频信号类型识别起始频点音频数据fft转化音频信号类型 1.输入信号为待测试模块发出的阶梯型频率的信号我们通过声卡采集,采集到的信号如下图: 我们通过exce ...
利用短时傅里叶变换（STFT）对信号进行时频谱分析和去噪声
利用短时傅里叶变换(STFT)对信号进行时频谱分析和去噪声 1.背景傅里叶变换(TF)对频谱的描绘是"全局性"的,不能反映时间维度局部区域上的特征,人们虽然从傅立叶变换能清楚地 ...

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