广州网站建设信科便宜,高端ppt模板,网站维护费用计入什么科目,做网站 提要求短时傅里叶变换及其逆变换 本篇文章主要记录了使用python进行短时傅里叶变换#xff0c;分析频谱#xff0c;以及通过频谱实现在频域内降低底噪的代码及分析#xff0c;希望可以给同样在学习信号处理的大家一点帮助#xff0c;也希望大家对我的文章多提意见建议。 一. 短…短时傅里叶变换及其逆变换 本篇文章主要记录了使用python进行短时傅里叶变换分析频谱以及通过频谱实现在频域内降低底噪的代码及分析希望可以给同样在学习信号处理的大家一点帮助也希望大家对我的文章多提意见建议。 一. 短时傅里叶变换与离散傅里叶变换 在这篇文章中我们主要运用了短时傅里叶变换要想清楚地理解短时傅里叶变换首先必须要了解离散傅里叶变换(Discrete Fourier TransformDFT)。 1.离散傅里叶变换 离散傅里叶的定义: ∀k∈[0,M−1],X[k]∑n0N−1x[n]e−2jπnfkFe∑n0N−1x[n]e−2jπnkFeMFe\\{\forall} k\in [0,M-1] ,X[k]\sum^{N-1}_{n0}x[n] e^{-\cfrac{2j\pi nf_k} {F_e}} \sum^{N-1}_{n0}x[n]e^{-\cfrac{2j\pi nk\frac{F_e}{M}}{Fe}} ∀k∈[0,M−1],X[k]n0∑N−1​x[n]e−Fe​2jπnfk​​n0∑N−1​x[n]e−Fe2jπnkMFe​​​ ∑n0N−1x[n]e−2jπknM\qquad \qquad \qquad \qquad \sum^{N-1}_{n0}x[n]e^{-\cfrac{2j\pi kn}{M}}∑n0N−1​x[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)∫R​x(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(NwinNhop)for j in range(L):xmat[:,j] np.fft.fft(trame[j*Nhop:NwinNhop*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,dtypecomplex)t_vecteur np.arange(0,Te*len(yvect),Te)index 0K 0L xmat.shape[1]yl np.zeros(xmat.shape,dtypecomplex)for j in range(L):yl[:,j] np.fft.ifft(xmat[:,j])# 平移和求和for k in range(L):yvect[Nhop*k:NfftNhop*k] yl[:,k]# 标准化幅值for n in range(Nwin-1):K window[n]K / Nhopyvect /Kreturn t_vecteur, yvect该算法的实现需要三步。 1. 快速傅里叶逆变换 yl np.zeros(xmat.shape,dtypecomplex) 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:NfftNhop*k] yl[:,k]对yl矩阵的每一列平移 (l-1)Nhopl ∈\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, yvectwindow[n] (w[n]) 是长度为Nwin的窗函数在选取窗函数的时候我们总满足规则 K∑l1Lw[n−(l−1)Nhop]\ K\sum^L_{l1}w[n-(l-1)Nhop] Kl1∑L​w[n−(l−1)Nhop] K的值与n无关。在此基础上不难证明 K≈∑n0Nwin−1w[n]/Nhop\ K \approx \sum^{Nwin-1}_{n0}w[n] / Nhop K≈n0∑Nwin−1​w[n]/Nhop 那么通过以上的三个步骤我们就可以从信号的短时傅里叶变换矩阵中完美重构原信号了。 下一篇文章我们将使用这些算法使用谱减法进行声音信号的降噪处理。