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

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

本篇文章主要记录了使用python进行短时傅里叶变换,分析频谱,以及通过频谱实现在频域内降低底噪的代码及分析,希望可以给同样在学习信号处理的大家一点帮助,也希望大家对我的文章多提意见建议。

一. 短时傅里叶变换与离散傅里叶变换

在这篇文章中我们主要运用了短时傅里叶变换,要想清楚地理解短时傅里叶变换,首先必须要了解离散傅里叶变换(Discrete Fourier Transform,DFT)。

1.离散傅里叶变换

离散傅里叶的定义:
∀ k ∈ [ 0 , M − 1 ] , X [ k ] = ∑ n = 0 N − 1 x [ n ] e − 2 j π n f k F e = ∑ n = 0 N − 1 x [ n ] e − 2 j π n k F e M F e \\{\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,M1],X[k]=n=0N1x[n]eFe2jπnfk=n=0N1x[n]eFe2jπnkMFe
= ∑ n = 0 N − 1 x [ n ] e − 2 j π k n M \qquad \qquad \qquad \qquad = \sum^{N-1}_{n=0}x[n]e^{-\cfrac{2j\pi kn}{M}} =n=0N1x[n]eM2jπ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 ) = ∫ R x ( t ) h ∗ ( t − τ ) e − 2 j π f t d t X(\tau,f)=\int_Rx(t)h^*(t-\tau)e^{-2j\pi ft}dt X(τ,f)=Rx(t)h(tτ)e2jπ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)+1
  M = Nfft
  xmat = 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

上述函数解释:
参数部分
trameFe : 初始的数字信号和它的采样频率
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/Fe
  yvect = np.zeros(Nfft + (xmat.shape[1]-1)*Nhop,dtype=complex)
  t_vecteur = np.arange(0,Te*len(yvect),Te)
  index = 0
  K = 0
  L = 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 /= Nhop
  yvect /=K
  return 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 /= Nhop
	yvect /=K
return t_vecteur, yvect

window[n] (w[n]) 是长度为Nwin的窗函数,在选取窗函数的时候,我们总满足规则   K = ∑ l = 1 L w [ n − ( l − 1 ) N h o p ] \ K=\sum^L_{l=1}w[n-(l-1)Nhop]  K=l=1Lw[n(l1)Nhop] K的值与n无关。在此基础上不难证明
  K ≈ ∑ n = 0 N w i n − 1 w [ n ] / N h o p \ K \approx \sum^{Nwin-1}_{n=0}w[n] / Nhop  Kn=0Nwin1w[n]/Nhop

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

已标记关键词 清除标记
相关推荐
©️2020 CSDN 皮肤主题: 点我我会动 设计师:白松林 返回首页