AcousticsML声学机器学习翻译教程一（信号处理Signal Processing）

源自：https://github.com/RAMshades/AcousticsM

信号处理

声学研究声音在环境中的传输、接收、产生及其影响。声音通常被想象为一种通过介质（如空气、液体等）振动传播的波。这种波可被可视化为与传播方向平行的压缩波，其特征由频率、波长、振幅和速度定义。声学在工程、音乐、环境科学、通信和健康等众多领域至关重要。我们可以利用信号处理来增强对声学的理解并拓展其研究。

信号处理通过多种算法实现声音的声学分析，使声学工作者能够分离、分类和操控声音。基于声学的信号处理算法与理论改变了我们在室内声学、结构声学、生物声学和水声学（仅举几例）等领域对声音的认知。尽管信号处理及其在声学中的应用理论已有悠久历史，本笔记本将从采样理论、数值表示和波传播基础入手，逐步介绍声音示例与可视化。通过这些技术，我们能更好地理解如何开发机器学习模型并将其应用于声学研究中。本笔记本旨在简要介绍波的表征形式及其与声学的关联。若您对声学有进一步探索的兴趣，我们建议您通过文末资源链接深入了解更多内容。笔记本涵盖的主题包括：

主题

1. 声音示例
2. 一维波
3. 离散与连续信号
4. 奈奎斯特采样定理
5. 快速傅里叶变换
6. 频谱图
7. 其他信号处理技术

您可能会问：信号处理如何影响声学领域的机器学习？这是一个至今仍在探索的有趣问题。机器学习的效果取决于数据点的数量或所使用的频率、时间等表征形式。多数情况下，二者相辅相成。虽然本笔记本不直接深入机器学习领域，但请通过这些内容熟悉相关基础概念。

创建者：Ryan A. McCarthy

导入包

在Notebook中，请运行下方单元格以加载所需的基础包。

# import packages
import numpy as np
import pandas as pd
import seaborn as sns
import matplotlib.pyplot as plt

import librosa 
import librosa.display
import IPython.display as ipd

from glob import glob

声音示例

让我们从加载、聆听和分析一段声学声音开始。下方是通过 Python 包 https://librosa.org/doc/latest/index.html# 提供的小号演奏者音频片段。

filename = librosa.ex('trumpet')
y, sr = librosa.load(filename)

声学声音的输出提供了音频时间序列和采样率。时间序列表示每个时间步长的声音振幅，其中时间步长（t）由以下公式决定：
$$t(i)=\frac{1}{S},$$
其中 S 是采样率（通常以每秒采样数为单位，除非另有说明）。因此，我们可以将采样率理解为每秒钟音频信号中包含的数据点数量（例如：若 S=180000，则一秒对应 180000 个采样点）。接下来让我们聆听并观察下方小号演奏者的时间序列可视化呈现。

# Play audio file
ipd.Audio(data=y,rate = sr)

plt.figure(figsize = (8, 2))
librosa.display.waveshow(y = y, sr = sr);
plt.title("Sound Waves", fontsize = 23);
plt.ylabel('Amplitude')
plt.show()

从上述声波中我们可以注意到振幅值。绝对振幅越高，声音越响亮。我们还能观察到振幅既可为正也可为负，但这是为什么呢？

这是因为声波的表征方式——我们可以将声音视为围绕基线振荡的正弦函数。基线越高，表示环境背景噪音越大。这种表征方式有助于理解声波相互作用时产生的干涉现象。接下来让我们探索正弦波的具体形态。

一维波 (1D Waves)

对波（或声音）y 最简单的表征是通过随时间 t 变化的一维振荡正弦函数来实现，其表达式为：

$$y(t)=A\sin (2\pi ft+\varphi )$$
其中 A 表示波在一个周期内最高点（波峰）与最低点（波谷）之间的最大振幅，f 是频率（单位为赫兹），$\varphi$ 是相位偏移。改变振幅会影响信号的高度，而改变频率则会调整信号随时间振荡的速率。当我们调整相位时，会沿着 x 轴（在一维情况下）将信号向右或向左平移。下文将通过示例演示：我们先改变正弦波的相位，再调整正弦波与余弦波的频率和振幅。

x = np.linspace(0, 1, 201)
steps = np.arange(0,2*np.pi,np.pi/2)
names = ['0',r'$\pi/2$',r'$3\pi/2$',r'$2\pi$']
for en,phi in enumerate(steps):
    plt.subplot(2, 2, en+1)
    plt.plot(x, np.sin(2*np.pi*x + phi))
    plt.xlabel('Angle [rad]')
    plt.ylabel('Amplitude')
    plt.axis('tight')
    plt.title(f'sin(2*pi ft + {names[en]})')
plt.tight_layout()
plt.show()

当相位为 $\pi /2$（即 $90^o$）时，我们生成的是余弦波。而当正弦波的相位达到 $3\pi /2$（即 $270^o$）时，该正弦波会呈现完全反射的特性。接下来我们将观察调整信号频率和振幅时产生的变化。

# Create x data
x = np.linspace(-np.pi, np.pi, 201)

# Plot sine and cosine and variations
plt.subplot(2, 3, 1)
plt.plot(x, np.sin(x))
plt.xlabel('Angle [rad]')
plt.ylabel('Amplitude')
plt.axis('tight')
plt.ylim(-2,2)
plt.title('sin(x)')

plt.subplot(2,3, 4)
plt.plot(x, np.cos(x))
plt.xlabel('Angle [rad]')
plt.ylabel('Amplitude')
plt.axis('tight')
plt.ylim(-2,2)
plt.title('cos(x)')

plt.subplot(2,3, 2)
plt.plot(x, 2*np.sin(x))
plt.xlabel('Angle [rad]')
plt.ylabel('Amplitude')
plt.axis('tight')
plt.ylim(-2,2)
plt.title('2*sin(x)')

plt.subplot(2,3, 5)
plt.plot(x, 2*np.cos(x))
plt.xlabel('Angle [rad]')
plt.ylabel('Amplitude')
plt.axis('tight')
plt.ylim(-2,2)
plt.title('2*cos(x)')

# change the frequency
plt.subplot(2, 3, 3)
plt.plot(x, np.sin(x*2))
plt.xlabel('Angle [rad]')
plt.ylabel('Amplitude')
plt.axis('tight')
plt.ylim(-2,2)
plt.title('sin(2*x)')

plt.subplot(2,3, 6)
plt.plot(x, np.cos(x*2))
plt.xlabel('Angle [rad]')
plt.ylabel('Amplitude')
plt.axis('tight')
plt.ylim(-2,2)
plt.title('cos(2*x)')

plt.tight_layout()
plt.show()

当我们将振幅调整为原来的2倍时，整体振幅便会翻倍。同样地，若将频率加倍，振荡次数也会翻倍。这个规律在我们持续提高频率或增大振幅时依然成立。接下来，我们将探讨这些波动现象的数学推导原理。

欧拉公式

正弦波最实用且便捷的表征形式之一是通过欧拉公式实现，其定义为：
$$e^{i\theta }=\cos (\theta )+i\sin (\theta )$$
其中 $e$ 是自然对数的底数，$i=\sqrt{(-1)}$，$\theta$ 为任意实数。该公式的优势在于允许将信号表示为复数形式，从而能够分析信号的各个分量（如相位角、频率和振幅）。若熟悉三角函数知识（抱歉这部分内容未在本笔记本中涵盖），我们可以推导出如下三角恒等式：
$$\cos (\theta )=\frac{e^{i\theta }+e^{-i\theta }}{2}$$
$$\sin (\theta )=\frac{e^{i\theta }-e^{-i\theta }}{2i}$$
我们将证明过程留给读者自行探索，但这些方程的重要价值在于：它们揭示了如何将三角函数改写为自然指数函数的求和形式（如下所示）。这一概念将为我们理解信号的频率表征（即快速傅里叶变换）提供重要基础。

x = np.linspace(0, 2*np.pi, 201)
plt.subplot(2, 2, 1)
plt.plot(x, np.sin(x))
plt.xlabel('Angle [rad]')
plt.ylabel('Amplitude')
plt.axis('tight')
plt.title('sin(x)')

plt.subplot(2,2, 2)
plt.plot(x, np.cos(x))
plt.xlabel('Angle [rad]')
plt.ylabel('Amplitude')
plt.axis('tight')
plt.title('cos(x)')

plt.subplot(2, 2, 3)
plt.plot(x,np.imag( (np.exp(1j*x) - np.exp(-1j*x))/2))
plt.xlabel('Angle [rad]')
plt.ylabel('Amplitude')
plt.axis('tight')
plt.title('Imag(e$^{ix}$)')

plt.subplot(2,2, 4)
plt.plot(x, np.real( (np.exp(1j*x) + np.exp(-1j*x))/2))
plt.xlabel('Angle [rad]')
plt.ylabel('Amplitude')
plt.axis('tight')
plt.title('Real(e$^{ix}$)')

plt.tight_layout()
plt.show()

角度关系

基于上述内容，我们可以利用欧拉公式在复平面中以单位圆形式展示正弦与余弦的关系。复平面的角度范围从0到$2\pi$（即0到360度），沿逆时针方向变化。若调整波的振幅，单位圆的半径会按比例相应变化。在单位圆的每个点上，我们可以观察到正弦和余弦如何随时间协同变化。如需进一步探讨或阅读，建议参考https://mathworld.wolfram.com/UnitCircle.html对单位圆的阐述。

x = np.linspace(0,2*np.pi,300)
re = np.cos(x)
im = np.sin(x)
plt.plot(re,im)
plt.plot(np.cos(0),np.sin(0),'rx')
plt.plot(np.cos(np.pi/2),np.sin(np.pi/2),'kx')
plt.plot(np.cos(np.pi),np.sin(np.pi),'bx')
plt.axhline(y=0,color='black',linewidth=0.5)
plt.axvline(x=0,color='black',linewidth=0.5)
plt.title('Complex Plane')
plt.ylabel('Imaginary')
plt.xlabel('Real')
plt.legend(['_','x=0','x=pi/2','x=pi'])
ax = plt.gca()
ax.set_aspect('equal', adjustable='box')
plt.show()

离散信号与连续信号

在了解一维波和波动理论的基础知识后，我们需要探讨音频领域的一些关键差异。首先，什么是离散信号和连续信号？为何要讨论这个话题？

连续信号是指在给定时间区间内（可能是无限区间）所有时间点上都存在的信号值。例如在空间中传播的光波或传入人耳的声波都属于连续信号，这类信号也被称为模拟信号。我们可以将模拟信号视为能够无限采样并完全还原的信号。但在现实世界中，我们可能无法获得如此高分辨率的数据。实际上，采集器接收到的模拟信号会被转换为数字信号——即用数字或比特编码的序列。这类信号称为离散信号，它们是在给定时间区间内离散点采样的数值。WAV音频文件就是典型实例。每秒能够编码的采样数量称为采样率（前文已讨论）。经验表明：采样率越高，可还原的频率范围就越大。但为何这一点如此重要呢？

机器学习中的离散信号

虽然通过机器学习处理连续信号极具价值，但现实中我们通过仪器和采样来测量声音在不同介质（例如水或空气）中的传播。当我们对信号进行采样时，可能会丢失某些重要或非重要信息。为更深入理解这一点，下文将通过余弦波示例演示采样点数量增减对信号表征的影响。

fs = 2
x = np.linspace(0, 2*np.pi, 401)
y = np.linspace(0,2*np.pi,4)

plt.plot(x,np.cos(fs*x))
plt.plot(y,np.cos(fs*y))
plt.xlabel('Samples')
plt.ylabel('Amplitude')
plt.show()

从以上示例可以看出，改变采样数量会导致信号信息的丢失。这将造成无法准确识别余弦波并错误改变其频率的现象，即所谓的混叠效应。混叠是一种因信号欠采样而产生的测量误差。为提升离散信号质量，我们需要增加采样点数；然而更高采样率意味着需要更大的设备存储空间或配备支持高速采样的设备。

奈奎斯特定理

在给定采样率下，我们能够成功还原哪些频率？恢复信号所需的最小采样数是多少？此时我们可以开始引入相关理论，特别是奈奎斯特定理。该定理表述如下：
$$S>2B$$
其中 B 表示信号带宽，$S$ 为采样率。通俗地说，在给定采样率 S（通常受仪器限制）下可还原的最高频率为 $\frac{S}{2}$。反之，若需要还原特定频率，仪器的采样率必须达到目标频率的两倍。让我们通过实际示例来观察这一原理。

fs = 2 # sampling rate
L = 4 # length of signal
x = np.linspace(0, L*np.pi, 401)
y = np.linspace(0,L*np.pi,fs*L+1)

plt.plot(x,np.cos(fs*x))
plt.plot(y,np.cos(fs*y))
plt.xlabel('Samples')
plt.ylabel('Amplitude')
plt.show()

以信号频率两倍的速率进行采样，我们可以获得相对粗略但准确的信号轮廓。

波叠加原理

在了解一维波分量、三角关系及采样定理后，让我们开始探究这些波的相互作用机制。假设有两个独立波同时传播，波的组合可视为这些波的线性求和。这种叠加会导致波产生干涉现象，形成建设性干涉和/或破坏性干涉。当两个波相位一致时会发生建设性干涉，导致振幅增大；反之，当两个波相位相反时会发生破坏性干涉，导致振幅减小。下面我们通过具体示例进行观察。

# Constructive interference example
N = 400
# sample spacing
T = 1.0 / 1000.0
x = np.linspace(0.0, N*T, N, endpoint=False)
y = np.real(np.exp(50.0 * 1.j * 2.0*np.pi*x)) #+ 0.5*np.exp(-80.0 * 1.j * 2.0*np.pi*x)
plt.subplot(3,1,1)
plt.plot(x,y)
plt.title('Frequency 1')

y1 = np.real(2*np.exp(50.0 * 1.j * 2.0*np.pi*x))
plt.subplot(3,1,2)
plt.plot(x,y1)
plt.title('Frequency 2')

y2 = y+y1
plt.subplot(3,1,3)
plt.plot(x,y2)
plt.title('Frequency 1 + Frequency 2')

plt.suptitle('Constructive Interference')
plt.tight_layout()

plt.show()

# Destructive interference example
x = np.linspace(0, 2*np.pi, 201)
y = np.real(np.cos(6*x)) 
plt.subplot(3,1,1)
plt.plot(x,y)
plt.title('Frequency 1')

y1 = np.real(np.cos(6*x-np.pi)) # pi is 180 degrees out of phase
plt.subplot(3,1,2)
plt.plot(x,y1)
plt.title('Frequency 2')

y2 = y+y1
plt.subplot(3,1,3)
plt.plot(x,y2)
plt.title('Frequency 1 + Frequency 2')
plt.ylim(-1,1)

plt.suptitle('Destructive Interference')
plt.tight_layout()

plt.show()

可以观察到：当波频率相同且相位一致时，振幅会产生叠加效应；而在后一种情况下，两个波相位相反时会产生完全相消干涉，导致振幅归零。接下来让我们观察频率变化时的情形。

# Varying frequencies
N = 400
# sample spacing
T = 1.0 / 1000.0
x = np.linspace(0.0, N*T, N, endpoint=False)
y = np.real(np.exp(50.0 * 1.j * 2.0*np.pi*x)) #+ 0.5*np.exp(-80.0 * 1.j * 2.0*np.pi*x)
plt.subplot(3,1,1)
plt.plot(x,y)
plt.title('Frequency 1')

y1 = np.real(0.5*np.exp(-80.0 * 1.j * 2.0*np.pi*x))
plt.subplot(3,1,2)
plt.plot(x,y1)
plt.title('Frequency 2')

y2 = y+y1
plt.subplot(3,1,3)
plt.plot(x,y2)
plt.title('Frequency 1 + Frequency 2')

plt.suptitle('Adding Different Frequencies')

plt.tight_layout()
plt.show()

两个波叠加后形成了一个形态奇特的信号，但仍有可能还原出原始的线性组合信号。

卷积还是相关？

在探讨如何识别和还原构成信号的频率成分之前，我们首先需要了解两个数学运算。信号处理中两个重要操作是卷积（convolution）与相关（correlation）。卷积通过整合两个信号的重叠部分来生成第三个信号，其运算过程包含将其中一个信号沿y轴翻转后平移，再与另一个信号进行积分。数学表达式如下：

$$(f\ast g)(t)={\int }_{-\inf }^{\inf }f(t)g(x-t)dt$$
本质上，该运算考察的是一个信号与另一个信号所有平移版本的重叠量之和。卷积在统计学、概率论以及（本文讨论的声学应用领域）中都是强大工具。下文我们将展示两个信号卷积的实例。

from scipy import signal

x = np.arange(256)/256
segment1 = np.sin(2*np.pi*x)
segment2 = segment1 + np.random.normal(0,2,len(segment1))
plt.subplot(3,1,1)
plt.plot(segment1)
plt.title('Signal 1')

plt.subplot(3,1,2)
plt.plot(segment2)
plt.title('Signal 2')

plt.subplot(3,1,3)
plt.plot(np.convolve(segment1,segment2))
plt.title('Convolution')

plt.tight_layout()

除了卷积之外，另一种实用的数学操作是互相关（cross-correlation）。尽管与卷积非常相似，互相关通过位移函数来识别两个信号之间的统计关系。其运算通过以下公式实现：

$$(f\star g)(t)={\int }_{-\infty }^{\infty }f(t)g(x-t)dx$$
与卷积类似，该操作将一个信号直接平移覆盖另一个信号（无需进行翻转）。这种方法通常应用于长时间序列中以识别特定特征，具体示例如下所示。

x = np.arange(256)/256
segment1 = np.sin(2*np.pi*x)
segment2 = segment1 + np.random.normal(0,2,len(segment1))
plt.subplot(3,1,1)
plt.plot(segment1)
plt.title('Signal 1')

plt.subplot(3,1,2)
plt.plot(segment2)
plt.title('Signal 2')

plt.subplot(3,1,3)
corr = signal.correlate(segment1,segment2)
lags = signal.correlation_lags(len(segment1), len(segment2))
plt.plot(lags,corr)
plt.title('Correlation')

plt.tight_layout()

快速傅里叶变换

如何解析上述信号并观察其频率成分的变化？幸运的是，我们可以通过一种有效的数学工具——快速傅里叶变换（FFT）——将信号从时域转换到频域进行分析。对于采样信号而言，FFT可视为离散傅里叶变换（DFT）的一种实现方式；本notebook中将统一使用FFT指代该变换。虽然不会深入探讨FFT的全部理论背景，但我们将阐述其基本原理与应用。给定长度为N的信号$x[n]$，FFT的数学定义如下：
$$X[k]=\sum _{n=0}^{N-1}x[n]e^{\frac{-j2\pi kn}{N}}$$
其中k = 0,…,n-1表示FFT的变换长度。这是FFT的离散形式。可以看出该信号通过自然指数函数及其不同频率分量进行表征（如前文所示）。当我们对周期性正弦波施加该变换时，将得到如下结果。
x[n]‌：时域离散信号的第n个采样点。
k: 频率分量，比如k=10，表示为10的频率
X[k]‌：频域第k个频率分量的复数值，表示幅度和相位。
N‌：总采样点数，通常为2的幂次（如1024）。

# Parameters:
fs = 64  # Sampling rate (samples per second)
T = 1  # Duration of the signal (seconds)
t = np.arange(0, T, 1/fs)  # Time vector
f = 10  # Frequency of the sine wave (Hz)
A = 4 # Amplitude
x = A*np.sin(2 * np.pi * f * t)

# Compute the FFT
X = np.fft.fft(x)
freqs = np.fft.fftfreq(len(x), 1/fs)

# Plot time data:
plt.subplot(2,1,1)
plt.plot(t,x,'*k')        
plt.plot(t,x) 
plt.title('Sinusoid')
plt.xlabel('Time (samples)') 
plt.ylabel('Amplitude')

# Plot the results
plt.subplot(2,1,2)
plt.stem(freqs, np.abs(X))
plt.xlabel('Frequency (Hz)')
plt.ylabel('Magnitude')
plt.title('FFT of a Sine Wave')

plt.tight_layout()
plt.show()

在上例中，离散傅里叶变换成功提取了10Hz频率的正弦波。幅度谱中出现在归一化频率位置的两个峰值符合我们的预期，因为：

$$DFT(\sin (\omega n))=A\sum _{n=0}^{N-1}\frac{e^{j0n}}{2}=A\frac{N}{2}$$
其中N为采样点数，A表示振幅。

傅里叶逆变换

使用傅里叶变换的优势在于我们可以通过逆变换还原出时间序列信号。从频域表征转换回时域信号的逆向傅里叶变换计算公式如下：

$$x[n]=\frac{1}{N}\sum _{k=0}^{N-1}X[k]e^{\frac{2\pi jkn}{N}}$$

应用离散傅里叶逆变换可精确还原我们最初创建的时域信号。

# Parameters:
fs = 64  # Sampling rate (samples per second)
T = 1  # Duration of the signal (seconds)
t = np.arange(0, T, 1/fs)  # Time vector
f = 10  # Frequency of the sine wave (Hz)
A = 4 # Amplitude
x = A*np.sin(2 * np.pi * f * t)

# Compute the FFT
X = np.fft.fft(x)
freqs = np.fft.fftfreq(len(x), 1/fs)

# Plot time data:
plt.subplot(3,1,1)
plt.plot(t,x,'*k')        
plt.plot(t,x) 
plt.title('Sinusoid')
plt.xlabel('Time (samples)') 
plt.ylabel('Amplitude')

# Plot the results
plt.subplot(3,1,2)
plt.stem(freqs, np.abs(X))
plt.xlabel('Frequency (Hz)')
plt.ylabel('Magnitude')
plt.title('FFT of a Sine Wave')

inversesignal = np.fft.ifft(X)
# Plot the results of Ifft
plt.subplot(3,1,3)
plt.plot(t,inversesignal,'*k')        
plt.plot(t,inversesignal) 
plt.title('IFFT of FFT')
plt.xlabel('Time (samples)') 
plt.ylabel('Amplitude')

plt.tight_layout()
plt.show()

在理解快速傅里叶变换（FFT）的基础原理后，让我们进一步探讨其应用中的复杂性问题。当我们对非周期性信号实施FFT时，会观察到频谱泄漏现象。频谱泄漏会导致特定频率的能量扩散至整个频谱范围。下文将通过示例演示频谱泄漏的具体表现形态。

fs = 64  # Sampling rate (samples per second)
T = .75  # Duration of the signal (seconds)
t = np.arange(0, T, 1/fs)  # Time vector
f = 10  # Frequency of the sine wave (Hz)
A = 4 # Amplitude
x = A*np.sin(2 * np.pi * f * t)

# Compute the FFT
X = np.fft.fft(x)
freqs = np.fft.fftfreq(len(x), 1/fs)

# Plot time data:
plt.subplot(2,1,1)
plt.plot(t,x,'*k')        
plt.plot(t,x) 
plt.title('Sinusoid')
plt.xlabel('Time (samples)') 
plt.ylabel('Amplitude')

# Plot the results
plt.subplot(2,1,2)
plt.stem(freqs, np.abs(X))
plt.xlabel('Frequency (Hz)')
plt.ylabel('Magnitude')
plt.title('FFT of a Sine Wave')

plt.tight_layout()
plt.show()

频谱泄漏与加窗技术

从上例可以看出，10Hz信号的能量已泄漏到其他频率，并扩散至整个频谱范围。这种泄漏会产生旁瓣——即频域中通常低于主瓣但遍布频谱的峰值。遗憾的是，对离散信号应用FFT时无法避免这种现象，但可以通过加窗技术缓解该问题。具体而言，频谱加窗法通过在信号段两端（左右侧）施加锥化函数，确保我们仅观察窗口内的瞬态事件。实际操作中，先对信号进行加窗处理，再将数据与窗函数值相乘。这有助于减少傅里叶变换可能产生的频谱泄漏。为简洁起见，我们以布莱克曼窗（Blackman window）为例进行说明。

布莱克曼窗旨在尽可能减少泄漏并抑制旁瓣（虽不能完全消除）。其计算公式为：
$$w=a_0-a_1\cos (\frac{2\pi n}{N})+a_2\cos (\frac{4\pi n}{N})$$
其中$a_n$为用于减弱旁瓣和频谱泄漏影响的常系数。具体参数请参阅SciPy库中https://docs.scipy.org/doc/scipy/reference/generated/scipy.signal.windows.blackman.html。此时FFT公式相应调整为：
$$X[k]=\sum _{n=0}^{N-1}x[n]w[n]e^{\frac{-j2\pi kn}{N}}$$
当时域卷积（或频域乘法）施加窗函数后，我们即可获得处理后的信号。

window = signal.windows.blackman(len(x))

# Compute the FFT
X = np.fft.fft(x*window)
freqs = np.fft.fftfreq(len(X), 1/fs)

plt.stem(freqs, np.abs(X))
plt.xlabel('Frequency (Hz)')
plt.ylabel('Magnitude')
plt.title('FFT of a Sine Wave')

plt.tight_layout()
plt.show()

我们可以看到10Hz正弦信号的能量泄漏已经减小。虽然我们未讨论所有窗函数，但我们建议读者查看其他窗函数，例如链接中提供的矩形窗（Boxcar）、凯泽窗（Kaiser）或巴特利特窗（Bartlett）。

窗函数 (scipy.signal.windows) — SciPy v1.16.2 手册
窗函数 (scipy.signal.windows
get_window (window, Nx[, fftbins, xp, device]) 返回给定长度和类型的窗函数。
barthann (M[, sym, xp, device]) 返回修正的巴特利特-汉恩窗。
bartlett (M[, sym, xp, device]) 返回巴特利特窗。
blackman (M[, sym, xp, device]) 返回布莱克曼窗。
blackmanharris (M[, sym, xp, device]) 返回最小4项布莱克曼-哈里斯窗。
bohman (M[, sym, xp, device]) 返回博曼窗。
boxcar (M[, sym, xp, device]) 返回矩形窗或箱式窗。
chebwin (M, at[, sym, xp, device]) 返回多尔夫-切比雪夫窗。
cosine (M[, sym, xp, device]) 返回简单余弦形状的窗函数。
dpss (M, NW[, Kmax, sym, norm, …]) 计算离散长球序列 (DPSS)。
exponential (M[, center, tau, sym, xp, device]) 返回指数（或泊松）窗。
flattop (M[, sym, xp, device]) 返回平顶窗。
gaussian (M, std[, sym, xp, device]) 返回高斯窗。
general_cosine (M, a[, sym]) 通用余弦项加权和窗函数。
general_gaussian (M, p, sig[, sym, xp, device]) 返回广义高斯形状的窗函数。
general_hamming (M, alpha[, sym, xp, device]) 返回广义汉明窗。
hamming (M[, sym, xp, device]) 返回汉明窗。
hann (M[, sym, xp, device]) 返回汉恩窗。
kaiser (M, beta[, sym, xp, device]) 返回凯泽窗。
kaiser_bessel_derived (M, beta, *[, sym, xp, …]) 返回凯泽-贝塞尔派生窗。
lanczos (M, *[, sym, xp, device]) 返回兰索斯窗（亦称 sinc 窗）。
nuttall (M[, sym, xp, device]) 根据纳托尔返回最小4项布莱克曼-哈里斯窗。
parzen (M[, sym, xp, device]) 返回帕曾窗。
taylor (M[, nbar, sll, norm, sym, xp, device]) 返回泰勒窗。
triang (M[, sym, xp, device]) 返回三角窗。
tukey (M[, alpha, sym, xp, device]) 返回图基窗（亦称 tapered cosine 窗）。

频谱图

由于声音并非总是连续的，FFT可能无法始终捕捉特定频率或实时分离声音。为解决这个问题，我们可以对信号进行分段并对每个分段施加FFT。通过这种方法，我们可以创建频谱图。频谱图包含振幅、频率和时间信息，对于分析瞬态声音非常有用。下图展示了小号声音的频谱图示例。

filename = librosa.ex('trumpet')
y, sr = librosa.load(filename)

# Default FFT window size
n_fft = 2048 # FFT window size
hop_length = 512 # number audio of frames between STFT columns (looks like a good default)

# Short-time Fourier transform (STFT)
D = np.abs(librosa.stft(y, n_fft = n_fft, hop_length = hop_length))

# Convert an amplitude spectrogram to Decibels-scaled spectrogram.
DB = librosa.amplitude_to_db(D, ref = np.max)

# Creating the Spectogram
plt.figure(figsize = (16, 6))
plt.subplot(1,2,1)
librosa.display.specshow(DB, sr = sr, hop_length = hop_length, x_axis = 'time', y_axis = 'log',
                        cmap = 'hot')
plt.colorbar();
plt.title('Log Scale')

plt.subplot(1,2,2)
librosa.display.specshow(DB, sr = sr, hop_length = hop_length, x_axis = 'time', y_axis = 'linear',
                        cmap = 'hot')
plt.colorbar();
plt.title('Linear Scale')

plt.show()

如先前所述，声音信号在时间轴上被分割成若干存在重叠的区段，每个区段的FFT计算结果会被拼接组合。通过这种方式，我们可以观察频率随时间的变化趋势以及某些主导频率的分布。上例展示了两种频谱图：y轴采用对数刻度的频谱图与采用线性刻度的频谱图。当关注低频区域时，使用对数刻度更容易识别频率特征；而对于高频分析，线性刻度通常更具实用性。

其他信号处理技术

平滑处理

平滑处理涉及修改信号中的数据点以消除噪声样本，这有助于识别趋势并提升信号质量。可采用多种平滑技术（包括https://docs.scipy.org/doc/scipy/reference/signal.html#filtering中提及的方法）。为简洁起见，我们将演示其中一种平滑方法：移动平均法。

移动平均平滑法

该方法通过提取信号（x）的窗口（例如10个离散点），计算这些点的平均值，并将该平均值作为当前点的值来消除噪声。平滑信号（y）的计算公式为：

$$y[t]=\frac{x_t+x_{t+1}+...+x_{t+L-1}}{L}$$
其中 t 表示时间索引，L 为窗口长度。如下图所示，我们在正弦信号中加入噪声，并尝试用此方法滤除噪声。

x = np.linspace(0, 2*np.pi, 301)
signal = np.cos(6*x)
noisysignal = signal + np.random.normal(0,.2,301)

# moving mean of 3 points
N = 3
newsignal3point = np.convolve(noisysignal, np.ones(N)/N, mode='same')
# moving mean of 7 points
N = 7
newsignal7point = np.convolve(noisysignal, np.ones(N)/N, mode='same')
# moving mean of 31 points
N = 31
newsignal31point = np.convolve(noisysignal, np.ones(N)/N, mode='same')

plt.subplot(4,1,1)
plt.plot(x,signal)
plt.plot(x,noisysignal)
plt.title('Original Signal')
plt.ylim(-1.5,1.5)

plt.subplot(4,1,2)
plt.plot(x,newsignal3point)
plt.title('3 Point Mean')
plt.ylim(-1.5,1.5)

plt.subplot(4,1,3)
plt.plot(x,newsignal7point)
plt.title('7 Point Mean')
plt.ylim(-1.5,1.5)

plt.subplot(4,1,4)
plt.plot(x,newsignal31point)
plt.title('31 Point Mean')
plt.ylim(-1.5,1.5)

plt.tight_layout()
plt.show()

需注意：当使用较小窗口进行平滑处理时，信号中仍残留部分噪声。相反，增大窗口时信号基本呈现平滑状态。但随着窗口尺寸持续增大，信号振幅开始衰减，同时可能消除信号中的高频成分（若存在）。

高通、低通与带通滤波器

本笔记本最后将讨论的是通带滤波技术，具体包含三种不同类型的滤波器：高通、低通和带通滤波器。这些滤波器根据对信号频率成分的保留或抑制特性进行区分。低通滤波器允许低频通过而抑制高频；高通滤波器则允许高频通过而抑制低频；带通滤波器则只保留特定频段内的频率（例如10至1000Hz之间）。多年来已发展出多种滤波器设计（详见https://docs.scipy.org/doc/scipy/reference/signal.html#filter-design），本节我们将重点介绍巴特沃斯滤波器（Butterworth filter）。

巴特沃斯滤波器的设计目标是在通带内呈现平坦的响应特性。对于频域中的给定信号，其传递函数（H）定义为：
$$H(\omega )=\frac{1}{\sqrt{1+(\frac{\omega }{{\omega }_c}{)}^{2n}}}$$
其中$\omega$是以弧度/秒为单位的角频率（满足$f2\pi =\omega$关系），${\omega }_c$为截止频率，n表示滤波器阶数。

# import package
from scipy import signal

低通滤波器

for i in range(1,10,2):
    b, a = signal.butter(i, 10, 'low', analog=True)
    w, h = signal.freqs(b, a)
    plt.semilogx(w, 20 * np.log10(abs(h)))
plt.title('Butterworth filter frequency response')
plt.xlabel('Frequency [radians / second]')
plt.ylabel('Amplitude [dB]')
plt.margins(0, 0.1)
plt.grid(which='both', axis='both')
plt.axvline(10, color='green') # cutoff frequency
plt.legend(['1','3','5','7','9'],title='Order')
plt.show()

上图展示了截止频率为10Hz的巴特沃斯滤波器频率响应特性。在截止频率处，幅值开始以显著速率下降（具体速率取决于所用阶数）。这种下降会抑制高频成分，如下方示例所示。

低通滤波器示例

t = np.linspace(0, 1, 1000, False)  # 1 second
sig = np.sin(2*np.pi*10*t) + np.sin(2*np.pi*20*t)
fig, (ax1, ax2) = plt.subplots(2, 1, sharex=True)
ax1.plot(t, sig)
ax1.set_title('10 Hz and 20 Hz sinusoids')
ax1.axis([0, 1, -2, 2])

sos = signal.butter(10, 15, 'lp', fs=1000, output='sos')
filtered = signal.sosfilt(sos, sig)
ax2.plot(t, filtered)
ax2.set_title('After 15 Hz Low-pass filter Butterworth')
ax2.axis([0, 1, -2, 2])
ax2.set_xlabel('Time [seconds]')
plt.tight_layout()
plt.show()

从上述示例中可观察到，低频信号得到了保留。

高通滤波器

要将上述方程调整为高通滤波器，我们可以进行以下变换：
$$H(\omega )=\frac{1}{\sqrt{1+(\frac{\omega }{{\omega }_c}{)}^{-2n}}}$$

for i in range(1,10,2):
    b, a = signal.butter(i, 10, 'highpass', analog=True)
    w, h = signal.freqs(b, a)
    plt.semilogx(w, 20 * np.log10(abs(h)))
plt.title('Butterworth filter frequency response')
plt.xlabel('Frequency [radians / second]')
plt.ylabel('Amplitude [dB]')
plt.margins(0, 0.1)
plt.grid(which='both', axis='both')
plt.axvline(10, color='green') # cutoff frequency
plt.legend(['1','3','5','7','9'],title='Order')
plt.show()

巴特沃斯滤波器会衰减截止频率以下的频率响应，同时保留高频成分。下方高通滤波器示例展示了这一特性。

高通滤波器示例

t = np.linspace(0, 1, 1000, False)  # 1 second
sig = np.sin(2*np.pi*10*t) + np.sin(2*np.pi*20*t)
fig, (ax1, ax2) = plt.subplots(2, 1, sharex=True)
ax1.plot(t, sig)
ax1.set_title('10 Hz and 20 Hz sinusoids')
ax1.axis([0, 1, -2, 2])

sos = signal.butter(10, 15, 'hp', fs=1000, output='sos')
filtered = signal.sosfilt(sos, sig)
ax2.plot(t, filtered)
ax2.set_title('After 15 Hz high-pass filter')
ax2.axis([0, 1, -2, 2])
ax2.set_xlabel('Time [seconds]')
plt.tight_layout()
plt.show()

从该滤波器的效果可见，它保留了信号的高频成分，同时滤除了低频部分。基于相同原理，我们可以在下例中应用带通滤波器进行验证。

带通滤波器

for i in range(1,10,2):
    b, a = signal.butter(i, [1,20], 'bandpass', analog=True)
    w, h = signal.freqs(b, a)
    plt.semilogx(w, 20 * np.log10(abs(h)))
plt.title('Butterworth filter frequency response')
plt.xlabel('Frequency [radians / second]')
plt.ylabel('Amplitude [dB]')
plt.margins(0, 0.1)
plt.grid(which='both', axis='both')
plt.axvline(1, color='green') # cutoff frequency
plt.axvline(20, color='green') # cutoff frequency
plt.legend(['1','3','5','7','9'],title='Order')
plt.show()

带通滤波器示例

t = np.linspace(0, 1, 1000, False)  # 1 second
sig = np.sin(2*np.pi*10*t) + np.sin(2*np.pi*20*t) + np.sin(2*np.pi*5*t)
fig, (ax1, ax2) = plt.subplots(2, 1, sharex=True)
ax1.plot(t, sig)
ax1.set_title('10 Hz and 20 Hz sinusoids')
ax1.axis([0, 1, -2, 2])

sos = signal.butter(10, [7,15], 'bp', fs=1000, output='sos')
filtered = signal.sosfilt(sos, sig)
ax2.plot(t, filtered)
ax2.set_title('After 7-15 Hz Band-pass filter')
ax2.axis([0, 1, -2, 2])
ax2.set_xlabel('Time [seconds]')
plt.tight_layout()
plt.show()

移除高频与低频部分后，重构信号保留了截止频带内的成分，最终得到上图所示信号。

Citations

1. Alan V. Oppenheim and Ronald W. Schafer. 2009. Discrete-Time Signal Processing (3rd. ed.). Prentice Hall Press, USA.
2. Blackman, R.B. and Tukey, J.W., (1958) The measurement of power spectra, Dover Publications, New York.
3. Harris, Fredric J. (Jan 1978). “On the use of Windows for Harmonic Analysis with the Discrete Fourier Transform”. Proceedings of the IEEE 66 (1): 51-83. DOI:10.1109/PROC.1978.10837.
4. Pauli Virtanen, Ralf Gommers, Travis E. Oliphant, Matt Haberland, Tyler Reddy, David Cournapeau, Evgeni Burovski, Pearu Peterson, Warren Weckesser, Jonathan Bright, Stéfan J. van der Walt, Matthew Brett, Joshua Wilson, K. Jarrod Millman, Nikolay Mayorov, Andrew R. J. Nelson, Eric Jones, Robert Kern, Eric Larson, CJ Carey, İlhan Polat, Yu Feng, Eric W. Moore, Jake VanderPlas, Denis Laxalde, Josef Perktold, Robert Cimrman, Ian Henriksen, E.A. Quintero, Charles R Harris, Anne M. Archibald, Antônio H. Ribeiro, Fabian Pedregosa, Paul van Mulbregt, and SciPy 1.0 Contributors. (2020) SciPy 1.0: Fundamental Algorithms for Scientific Computing in Python. Nature Methods, 17(3), 261-272. DOI: 10.1038/s41592-019-0686-2.
5. Proakis, J. G., & Manolakis, D. G. (1996). Digital signal processing (3rd ed.): principles, algorithms, and applications. USA: Prentice-Hall, Inc.