python 每个频段的功率谱密度 python频域分析

pyAudioKits是基于librosa和其他库的强大Python音频工作流支持。

API速查手册

通过pip安装：

pip install pyAudioKits

本项目的GitHub地址，如果这个项目帮助到了你，请为它点上一颗star，谢谢你的支持！如果你在使用过程中有任何问题，请在评论区留言或在GitHub上提issue，我将持续对该项目进行维护。

上节中对乐音和噪音的分析使用的方式均为时域分析的方法。本节我们将引入频域的概念，并介绍如何在频域分析音频信号。

傅里叶变换

连续时间信号傅里叶变换

我们说连续时间周期信号可以展开为频率为$\omega_0$的基波和频率为$(k+1)\omega_0（k=1,2,...）$的谐波。也就是说我们不仅可以使用原来的$x(t)$来表示连续时间周期信号，只要知道基波和谐波的频率以及振幅，就可以用频率和振幅来表示连续时间周期信号了。因此，我们现在假设有一个频率轴，上面有无数的频点$\omega$，就可以尝试将原本时间轴上的一维信号$x(t)$映射到频率轴上，变成$X(j\omega)$，完成时域到频域的映射。

使用最小二乘法，保证时域表示和频域表示能量相差最小。为了达成这一条件的结论是：设所有频点的值为$a_k$乘上$2\pi$。表示成频谱密度函数，有$X(j\omega)=\displaystyle\sum_{k=-∞}^∞2\pi a_k\delta(\omega-k\omega_0)$。可见周期信号的频谱密度函数是离散的冲激函数线性叠加。该函数可以通过$x(t)=\frac{1}{2\pi}\displaystyle\int_{-∞}^∞X(j\omega)e^{j\omega t}d\omega$还原为时域中的傅里叶级数表示$x(t)=\displaystyle\sum_{k=-∞}^∞a_ke^{jk\omega_0t}$

那么非周期信号能不能映射到频域呢？答案是肯定的。假如基波频率$\omega_0$无限小，即基波周期$T$无限大，则可以表示任意（无论是否周期）信号的频域特性，此时就从傅里叶级数推广到了傅里叶变换：基波频率$\omega_0$无限小，因此$k\omega_0$视为连续值$\omega$，$X(j\omega)=\displaystyle\int_{-∞}^∞x(t)e^{-j\omega t}dt$。此时得到连续时间傅里叶变换的综合方程为$x(t)=\frac{1}{2\pi}\displaystyle\int_{-∞}^∞X(j\omega)e^{j\omega t}d\omega$，分析方程为$X(j\omega)=\displaystyle\int_{-∞}^∞x(t)e^{-j\omega t}dt$。可见非周期信号的频谱密度函数是连续函数。

当连续时间周期信号的周期逐渐扩大，$\omega_0$就逐渐缩小，组成$X(j\omega)$的冲激函数就越密集。当周期无限大时，$\omega_0$无限小，$X(j\omega)$成为连续函数，就得到了连续时间非周期信号频谱密度函数。

离散时间信号傅里叶变换

可是我们在计算机中储存的信号都是数字信号，也就是离散时间信号，它们怎么映射到频域呢？

还记得模拟信号到数字信号（在我们这里也就等价于连续时间信号到离散时间信号）使用了一个采样的操作。信号在时域乘上冲激串$p(t)=\displaystyle\sum_{n=-∞}^{+∞}\delta(t-nT_s)$可以实现采样。采样后的时域表示为$x_s(t)=x(t)p(t)=\displaystyle\sum_{n=-∞}^{+∞}x(t)\delta(t-nT_s)$。采样后的频域表示为$X_s(j\omega)=\frac{1}{2\pi}X(j\omega)*P(j\omega)=\frac{1}{T_s}\displaystyle\sum_{k=-∞}^{+∞}X(j(\omega-k\omega_s))$，我们发现此时$X_s(j\omega)$是$X(j\omega)$在频谱上以采样频率为基数的周期性延拓的结果，其周期为$\omega_s=2\pi/T_s$且一个周期内就包含了所有的频谱信息。此时得到离散时间信号傅里叶变换的综合方程为$x[n]=\frac{1}{2\pi}\displaystyle\int_{\frac{2\pi}{T_s}}X(j\omega)e^{j\omega nT_s}d\omega$，分析方程为$X(j\omega)=\displaystyle \sum_{n=-∞}^{+∞}x[n]e^{-j\omega nT_s}$

此时综合方程和分析方程与采样周期$T_s$是有关的。为了使它们与采样周期$T_s$无关，我们将频域也用因子$T_s$（即采样率的倒数）归一化，记为$X(e^{j\omega})=X(j\omega T_s)$，得到离散时间信号傅里叶变换的综合方程为$x[n]=\frac{1}{2\pi}\displaystyle\int_{2\pi}X(e^{j\omega})e^{j\omega n}d\omega$，分析方程为$X(e^{j\omega})=\displaystyle \sum_{n=-∞}^{+∞}x[n]e^{-j\omega n}$。此时$X(e^{j\omega})$中，$\omega$的每$1 rad/s$实际代表频率为$\frac{1}{2\pi T_s}Hz(\frac{1}{T_s}rad/s)$。$X(e^{j\omega})$以$2\pi$为周期。

离散傅里叶变换

现在问题又来了，还记得我们为什么要对模拟信号$x(t),-∞<t<∞$采样变成$x[n],0≤n<N_{max}$吗？这是因为计算机无法储存原来的模拟信号，因为模拟信号是连续的，可能还是无限的、非因果的。那么频谱密度函数$X(e^{j\omega}),-∞<\omega<∞$实际上也面临着一样的问题。我们也需要对这个频谱密度函数进行一些“采样”操作。

令周期冲激串$\tilde p[n]=\displaystyle\sum_{r=-∞}^∞\delta[n-rN]$，则有$\tilde x[n]=x[n]*\tilde p[n]$，此时$\tilde X(e^{j\omega})=X(e^{j\omega})\tilde P(e^{j\omega})=X(e^{j\omega})\displaystyle\sum_{k=-∞}^∞\frac{2\pi}{N}\delta(\omega-\frac{2\pi k}{N})$，则有$X[k]=X(e^{j\omega})|_{\omega=\frac{2\pi k}{N}}$，从而在频域完成采样。使用$W_N$代替$e^{-j\frac{2\pi}{N}}$有：

$X[k]=\displaystyle \sum_{n=0}^{N-1} \tilde x[n]e^{-jk(\frac{2\pi}{N})n}=\sum_{n=0}^{N-1} \tilde x[n]W_N^{kn}$

$\tilde x[n]=\frac{1}{2\pi}\displaystyle\int_{2\pi}X(e^{j\omega})|_{\omega=\frac{2\pi k}{N}}e^{j\omega n}d\omega=\displaystyle\frac{1}{N}\displaystyle\sum_{k=0}^{N-1} X[k]e^{jk(\frac{2\pi}{N})n}=\frac{1}{N}\displaystyle\sum_{k=0}^{N-1} X[k]W_{N}^{-kn}$

而时域以$N$为周期进行周期延拓变为$\tilde x[n]$。原本有$x[n]=\begin{cases}x_d[n]\quad 0≤n<N_{max}\\0\quad others\end{cases}$，现在则有$\displaystyle\tilde x[n]=\sum_{r=-∞}^∞x[n-rN]$。周期延拓后的信号$\tilde x[n]$可能拥有和$x(t)$不同的频率成分。

原本$X(e^{j\omega})$中，$\omega$的每$1 rad/s$实际代表频率为$\frac{1}{2\pi T_s}Hz(\frac{1}{T_s}rad/s)$，现在$\tilde X[k]$中每个$k$实际代表$\frac{1}{T_sN}Hz(\frac{2\pi}{T_sN}rad/s)$的频率。称每个k为频谱中的一个频点，实际代表的频率为频谱分辨率。采样不能完整地表示出频谱，有些时候会遗漏频谱中的频率成分，这被称为栅栏效应。然而时频域依然满足能量相等的条件，被遗漏的频率成分的能量就会被泄露到附近的频点上，从而导致信号出现新的频率成分，这被称为频谱泄漏。

还记得$X(e^{j\omega})$的周期为$\frac{2\pi}{T_s}$吗？这意味着我们只需要取频谱中连续的$N$个频点即可得到频谱的所有信息。加上实信号的频谱关于$\omega=0$对称，我们观察频谱时只需要取频谱中连续的$N/2$个频点即可，即绘制频谱时一般只取$X[k],0≤k<\frac{N}{2}$。

由于$X[k]$为复数，频谱还分为幅度谱和相位谱，分别为$|X[k]|$和$\angle X[k] (rad/s)$。

def periodicExtension(p,N): #模拟进行3次周期延拓的函数。p为需要周期延拓的信号，N则是周期延拓的周期
    L = 3 * N
    p = ak.concatenate([ak.create_Single_Freq_Audio(0,0,p.sr,L/p.sr),p,ak.create_Single_Freq_Audio(0,0,p.sr,L/p.sr)])
    p1 = p[0:L]
    for i in range(N,len(p.samples)-L,N):
        p1 = p1 + p[i:i+L]
    return p1

import pyAudioKits.analyse as aly

p1 = ak.create_Single_Freq_Audio(0.1,8,256,1)   #创建振幅为0.1，频率为8Hz，采样率为256Hz，持续时间1s的正弦信号
p1.plot()

aly.FFT(p1,16).plot()   #进行16点傅里叶变换，并默认绘制幅度谱

x[n]采样自8Hz的正弦波x(t)，在幅度谱上却体现不出8Hz的峰值。因为此时频谱分辨率为$256/16=16Hz$，大于信号频率8Hz，因此出现栅栏效应和频谱泄漏。

periodicExtension(p1,16).plot()

同时在时域我们也可以看到，$x[n]$在以16点为周期进行周期延拓得到的$\tilde x[n]$不再以0.125s为周期。

aly.FFT(p1,32).plot()

使用32点傅里叶变换时，此时频谱分辨率为$256/32=8Hz$，正好等于信号频率8Hz。因此我们可以看到8Hz的峰值被体现了出来，不会观察到栅栏效应和频谱泄漏。

periodicExtension(p1,32).plot()

从时域上我们可以看到$\tilde x[n]$依然以0.125s为周期。

使用不同的频率轴刻度来观察真实频率$f$、以角频率表示的真实频率$\omega_s$、归一化频率$\omega$、频点的对应关系。

aly.FFT(p1,32).plot(xlabel="frequency/(rad/s)") #使用角频率表示的频率

aly.FFT(p1,32).plot(xlabel="normalized frequency/(rad/s)")  #归一化角频率，0~采样率/2 就对应 0~π

aly.FFT(p1,32).plot(xlabel="freq point")	#频点。32点傅里叶变换在[0,采样率/2)之间有16个频点

再来看信号中有两个频率成分（8Hz和32Hz的情况）。

p2 = p1 + ak.create_Single_Freq_Audio(0.1,32,256,1)
p2.plot()

aly.FFT(p2,16).plot()

使用16点FFT时，频谱分辨率为$256/16=16Hz$，小于信号频率32Hz且可以被32Hz整除，但大于信号频率8Hz。因此在32Hz处不会出现频谱泄漏和栅栏效应，但在8Hz处会。

aly.FFT(p2,32).plot()

使用32点FFT时，频谱分辨率为$256/32=8Hz$，小于信号频率且信号频率是频谱分辨率的整数倍，此时可以分辨出8Hz和32Hz的峰值，不会出现栅栏效应和频谱泄漏。

aly.FFT(p2,64).plot()

使用64点FFT时，频谱分辨率为$256/64=4Hz$，频谱进一步细化。

aly.FFT(p2).plot()

默认情况下，使用全256点进行FFT，频谱分辨率为$256/256=1Hz$，不会产生频谱泄漏和栅栏效应。

aly.FFT(p2,512).plot()

periodicExtension(p2,512).plot()

使用512点FFT相当于在256个样本后补256个0再进行周期延拓。此时信号中出现新的频率成分，因此产生频谱泄漏和栅栏效应。

p3 = p2[0:48]
aly.FFT(p3).plot()

取前48个样本点对信号进行截断，此时使用全48点进行FFT，频谱分辨率为$256/48=5.3Hz$。$5.3Hz$可以被32Hz整除，但不可以被8Hz整除。因此在32Hz处不会出现频谱泄漏和栅栏效应，但在8Hz处会。

p2.plot()

periodicExtension(p3,48).plot()

周期延拓的结果也显示了新的频率成分。

aly.FFT(p3,32).plot()

此时使用32点FFT依然不会观察到频谱泄漏和栅栏效应。

periodicExtension(p3,32).plot()

周期延拓时不会出现新的频率成分。

aly.FFT(p3,256).plot()

由于本身信号被截断为48点，再使用256点FFT也只是对补0后的信号进行FFT。因此必然出现新的频率成分。

说了这么多大家应该也已经发现了，从时域上看周期延拓带来的新的频率成分，和从频域上看由栅栏效应带来的频谱泄漏，其实是同一现象的两种解释罢了。那么有没有一个统一的方法量化这一影响呢？就是说所谓新的频率成分到底是多少？

还记得我们一开始对音频信号进行过截断吗？截断的时候用到了窗函数$w[n]$。$x[n]=x_d[n]w[n]$，所以$X(e^{j\omega})=\frac{1}{2\pi}\displaystyle\int_{2\pi}X_d(e^{j\theta})W(e^{j(\omega-\theta)})d\theta$，也就是说$x[n]$的频谱是$x_d[n]$的频谱与窗函数频谱的卷积。$w[n]=\begin{cases}1\quad 0≤n<N_{max}\\0\quad others\end{cases}$的频谱为$W(e^{j\omega})=\frac{1-e^{-j\omega N_{max}}}{1-e^{-j\omega}}$。假设$x_d[n]=sin\omega_0n$采样自正弦波$x(t)=sin\omega_0t$，其频谱为$X_d(e^{j\omega})=\frac{\pi}{j}\displaystyle\sum_{l=-∞}^{+∞}[\delta(\omega-\omega_0-2\pi l)-\delta(\omega+\omega_0-2\pi l)]$。则$X(e^{j\omega})=\frac{\pi}{j}\displaystyle\sum_{l=-∞}^{+∞}[\frac{1-e^{-j(\omega-\omega_0-2\pi l) N_{max}}}{1-e^{-j(\omega-\omega_0-2\pi l)}}+\frac{1-e^{-j(\omega+\omega_0-2\pi l) N_{max}}}{1-e^{-j(\omega+\omega_0-2\pi l)}}]$。如果采样取$\omega=\frac{2\pi k}{N}$，则有$\begin{aligned}X[k]=X(e^{j\omega})|*_{\omega=\frac{2\pi k}{N}}&=\frac{\pi}{j}\displaystyle\sum_*{l=-∞}^{+∞}[\frac{1-e^{-j(\frac{2\pi k}{N}-\frac{2\pi k_0}{N}-2\pi l) N_{max}}}{1-e^{-j(\frac{2\pi k}{N}-\frac{2\pi k_0}{N}-2\pi l)}}+\frac{1-e^{-j(\frac{2\pi k}{N}+\frac{2\pi k_0}{N}-2\pi l) N_{max}}}{1-e^{-j(\frac{2\pi k}{N}+\frac{2\pi k_0}{N}-2\pi l)}}]\\&=\frac{\pi}{j}\displaystyle\sum_{l=-∞}^{+∞}[\frac{1-e^{-j(\frac{2\pi k}{N}-\frac{2\pi}{N_0}-2\pi l) N_{max}}}{1-e^{-j(\frac{2\pi k}{N}-\frac{2\pi}{N_0}-2\pi l)}}+\frac{1-e^{-j(\frac{2\pi k}{N}+\frac{2\pi}{N_0}-2\pi l) N_{max}}}{1-e^{-j(\frac{2\pi k}{N}+\frac{2\pi}{N_0}-2\pi l)}}]\end{aligned}$，其中$k_0=\frac{\omega_0N}{2\pi}$为$X[k]$中代表$\omega_0$的频点，$N_0=\frac{2\pi}{\omega_0}$为以$x[n]$中样本点数表示的正弦波周期。这样，当$N_{max}$是$N$的整数倍且$N_{max}$是$N_0$整数倍时，$X[k]$仅在$\pm k_0+Nl,l\in Z$处有值，其余部分都是0，此时窗函数的频谱被隐藏了，我们只会看到$x_d[n]$的频谱；其他情况下，窗函数的频谱则会被显示出来。

也就是说，想要屏蔽窗函数的频谱，就必须满足：

1. 截断点数是离散傅里叶变换点数的整数倍
2. 截断点数是信号所有分量的周期的整数倍

不满足这些条件的情况下，我们就会在频谱上观测到窗函数的频谱，这就是周期延拓带来的“新的”频率成分，以及频谱泄漏产生的“新的”频率成分。第二个条件是非常苛刻的，实际应用中往往无法满足。因此观测到窗函数的频谱是非常非常正常的现象。

input1 = ak.create_Single_Freq_Audio(0.02,8,256,1)
aly.FFT(input1,256).plot(0, 256)

对于采样率为$256Hz$的音频，取傅里叶变换的点数$N=256$以计算幅度谱，则每个频点代表$1Hz$的真实频率，幅度谱的频点$k=0$到$k=255$就涵盖了$0Hz$到$255Hz$上所有的功率信息。我们发现幅度谱上有两个峰值，其中左边的峰值出现在单音的频率$8Hz$对应的频点上。实信号的幅度谱关于$\omega=0$偶对称，因此在对应$-8Hz$的频点上也会出现一个峰值，而这个峰值经过以采样率$256Hz$为周期的周期延拓后，就出现在了$256 - 8 = 248 Hz$对应的频点上，即功率谱上右边的峰值。

import numpy as np

aly.FFT(input1,256).plot(0, 2*np.pi, xlabel="normalized frequency/(rad/s)")

若使用归一化频率来作为横坐标，则我们发现幅度谱关于$\omega=\pi$对称。事实上，如果我们绘制$\omega\in (-∞,∞)$的幅度谱，则还会发现幅度谱以$\omega=2\pi$为周期，且关于$\omega=n\pi, n\in \mathbb Z$偶对称。

import numpy as np

aly.FFT(input1,256).plot(0, 2*np.pi, plot_type = "phase", xlabel="normalized frequency/(rad/s)")

实信号的相位谱也具有对称的特性，且也以$\omega=2\pi$为周期，但它是关于$\omega=n\pi, n\in \mathbb Z$奇对称的。

能量谱和功率谱

除了频谱之外，能量谱和功率谱也是能够在频域描述信号特点的函数。能量谱和功率谱的理论基础是帕塞瓦尔定理，即时域内信号能量等于频域内信号能量。对于能量有限信号而言，根据帕塞瓦尔定理，其能量为$E=\displaystyle\int_{-∞}^∞x^2(t)dt=\int_{-∞}^∞|X(j\omega)|^2d\omega$，而能量谱密度就为$G(j\omega)=|X(j\omega)|^2$。对于功率有限信号而言，首先将$x(t)$截短为长度等于T的一个截短信号$x_T(t),-T/2<t<T/2$，使其成为能量有限信号。根据帕塞瓦尔定理，其能量为$E=\displaystyle\int_{-∞}^{∞}x_T^2(t)dt=\int_{-∞}^{∞}|X_T(j\omega)|^2d\omega$，而信号功率为$P=\displaystyle\lim_{T\rightarrow ∞}\frac{1}{T}\int_{-∞}^{∞}|X_T(j\omega)|^2d\omega$，则有功率谱密度$P(j\omega)=\displaystyle\lim_{T\rightarrow ∞}\frac{1}{T}|X_T(j\omega)|^2$。

由于我们对音频信号进行了加窗截断，因此所有信号都是能量和功率有限信号。我们仅关心其功率谱密度，计算式为$P[k]=\frac{1}{N_{max}}|X[k]|^2$。

import matplotlib.pyplot as plt

input1 = ak.create_Single_Freq_Audio(0.02,8,256,1)
aly.PSD(input1, 256).plot(0,32,xlabel="frequency/Hz")

功率谱也可以使用增益的形式来表示。

aly.PSD(input1, 256, dB=True).plot(0,32,xlabel="frequency/Hz")

当信号中同时有两个频率成分时，两个频率成分的功率信息都会被反映到功率谱上。

input2 = ak.create_Single_Freq_Audio(0.01,24,256,1)
aly.PSD(input1 + input2, 256).plot(0,32,xlabel="frequency/Hz")

其中由于input2的振幅是input1的一半，所以$16Hz$上的功率谱密度是$8Hz$上的四分之一。

噪音的频域特性

和确定性周期信号以及具有明显音高的乐音不同，噪音是没有明显音高，即没有明显周期性的。乐音可以在幅度谱和功率谱上体现出明显的峰值，而噪音的频率成分则呈现连续而随机的分布。

白噪音的频率成分均匀分布在整个频率范围内。

monotone_without_wn = ak.create_Single_Freq_Audio(0.02,440,44100,1)
monotone_with_wn = monotone_without_wn.addWgn(10)   #在单音上加上信噪比为10dB的高斯白噪音
wn = monotone_with_wn - monotone_without_wn #减去单音，只留下白噪音
aly.FFT(wn).plot(), aly.PSD(wn).plot()