Python线性代数数字图像和小波分析之二

要点

1. 数学方程：数字信号和傅里叶分析，离散时间滤波器，小波分析
2. Python代码实现及应用变换过程：
   1. 读取音频和处理音频波，使用Karplus-强算法制作吉他音频
   2. 离散傅里叶计算功能和绘制图示结果
   3. 计算波形傅里叶系数
   4. 正向和反向（逆）快速傅里叶FFT 实现和使用
   5. 位反转函数
   6. 正向和反向（逆）离散余弦变换
   7. 离散时间滤波器结果绘图
   8. 一维、二维和三维正向和反向（逆）离散小波变换
   9. 张量积将小波应用于图像，二维和三维张量积实现

离散傅里叶分析数字声音

数字声音

数字声音是一个序列 x = { x i } i = 0 N − 1 \boldsymbol{x}=\left\{x_i\right\}_{i=0}^{N-1} x={xi​}i=0N−1​，对应于记录的$a$连续声音$f$的测量值以每秒 $f_s$ 测量的固定速率，即
$$x_k=f\left(k/f_s\right),\text{ 对于 }k=0,1,\dots ,N-1\text{. }$$
$f_s$ 称为采样频率或采样率。数字声音中的组成部分称为样本，连续样本之间的时间称为采样周期，通常表示为$T_s$。测量声音也称为对声音进行采样。

数字声音的质量通常通过比特率（每秒的位数）来衡量，即采样率与用于存储每个样本的位数（二进制数字）的乘积。 采样率和数字格式都会影响最终声音的质量。 这些被封装在数字声音格式中。

观察Python合成声音

首先，合成是使用各种方法创建声音的行为。 我将重点关注加法合成，即通过将不同信号相加的方式进行合成。

我将使用的库是 Numpy 和 Scipy。 Numpy 是一个用于数值计算的库，而 Scipy 是一个我将用于信号处理并将数据转换为声音文件的库。 我还将安装 matplotlib 来绘制声音信号。 Pydub 是一个用于播放音频的库。 为了安装这些包，我将使用 pip 并创建一个虚拟环境。我使用的是 Python 3.10，因此请调整您的 Python 安装说明。 在终端中，我们首先创建一个目录，然后创建一个 virtualenv 并安装必需的软件包。

请创建一个名为 sound.py 的新文件并在文本编辑器中打开它。这将是我用来合成声音的主文件。现在让我首先添加需要的导入。

import numpy as np
from scipy.io.wavfile import write
from scipy import signal
import matplotlib.pyplot as plt

让我们来回顾一下什么是声音。 声音是空气穿过空间的压力波。 它是一种称为纵波的波，意味着它垂直于传播方向移动。

声音可以表示为信号。 信号是一组随时间变化的数字，通常用方括号表示，例如 $x[n]$。 我们可以将信号设置为始终等于 1，例如 $x[n]$=1，或者创建一个像 $x[n]$ = $n$ 这样变化的线性信号。 出于我们的目的，我们将处理像正弦这样的周期性信号。 它们的形式为 $x[n]$ = $Asin(2\pi fn)$，其中 $A$ 是振幅，$f$ 是频率。 信号的 $y$ 轴是幅度或音量，$x$ 轴是当前时间。

让我们尝试使用以下方程生成正弦波。
$$y[n]=A\sin (2\pi ft[n])$$
A 是幅度，f 是信号的频率。 频率是每秒的样本数。 它决定了声音的音调，较低的频率有更多的低音，较高的频率有更多的高音。 为该声音生成波形文件的代码如下所示。

AUDIO_RATE = 44100
freq = 440
length = 1# create time values
t = np.linspace(0, length, length * AUDIO_RATE, dtype=np.float32)
# generate y values for signal
y = np.sin(2 * np.pi * freq * t)
# save to wave file
write("sine.wav", AUDIO_RATE, y)

我们首先有音频速率常数，也称为采样率，它是声卡从麦克风或其他音频输入源采样音频的速率。 在本例中，频率将为 44,100Hz。 赫兹 (Hz) 是频率的标准度量，是大多数声卡上每秒的样本数。 其原因是人类听觉的限制和奈奎斯特采样定理。 其中指出，要播放频率最多为 f 的声音，我们至少需要 2f 个样本。 人的听觉通常敏感度为20,000Hz，这意味着采样率需要在40,000Hz以上。

这意味着我们需要在几秒内生成相当于音频剪辑长度 44,100 倍的样本数量。 在本例中，我们需要 1 秒的音频剪辑，因此需要将音频速率乘以 1。这将使我们能够生成跨越人类听觉范围的声音。

我们使用 linspace 生成时间值，范围从 0 到 1，样本数为 1 * AUDIO_RATE，即 44100 个样本。 然后使用每个时间值的 sin 计算信号。 频率为440Hz，即中A。然后我们可以使用Scipy中的write函数和相应的AUDIO_RATE将声音写入波形文件。

如果您听 sine.wav，您会听到短促的一秒正弦波声音。 这听起来像是扬声器发出蜂鸣声，如果您听不到它，请检查以确保代码正确并且您的音频设置也正确。 您的音频格式可能需要在设置中设置为 44,100Hz。

让我们使用 matplotlib 中的绘图函数绘制信号。

def plot(ts, ys, title, num_samples):plt.xlabel("t")plt.ylabel("y")plt.title(title)plt.plot(ys[:num_samples])plt.show()plot(t, y, "Sine Signal", 512)

如果我们绘制前 512 个样本，这就是我们得到的图。

方波类似，但我们使用方形而不是正弦曲线。这使用阶跃函数代替正弦函数来计算。我们可以使用正弦函数和 Heaviside 阶跃函数来实现这一点。

y = np.heaviside(np.sin(2 * np.pi * freq * t), 1.0)

我们可以使用 Scipy 的 signal.square 函数实现同样的效果。

t = np.linspace(0, length, length * AUDIO_RATE, dtype=np.float32)# generate y values for signal
y = signal.square(2 * np.pi * freq * t)
# save to wave file
write("square.wav", AUDIO_RATE, y)plot(t, y, "Square Signal", 512)

方波会产生比正弦波更刺耳的嗡嗡声。这是因为方波比只有一种音调的正弦波有更多的泛音。

离散傅里叶分析

在离散傅里叶分析中，向量 $\mathbf{x}=\left(x_0,\dots ,x_{N-1}\right)$ 表示为 $N$ 个向量的线性组合
$${\varphi }_n=\frac{1}{\sqrt{N}}\left(1,e^{2\pi in/N},e^{2\pi i2n/N},\dots ,e^{2\pi ikn/N},\dots ,e^{2\pi in(N-1)/N}\right)$$
这些向量称为归一化复指数，或 $N$ 阶纯数字音调。 $n$也称为频率指数。整个集合 F N = { ϕ n } n = 0 N − 1 \mathcal{F}_N=\left\{\phi_n\right\}_{n=0}^{N-1} FN​={ϕn​}n=0N−1​称为$N$点傅里叶基。

请注意，纯数字音可以被视为纯音的样本，在一个周期内均匀采集： $f(t)=e^{2\pi int/T}/\sqrt{N}$ 是具有频率的纯音$n/T$，其样本为
$$f(kT/N)=\frac{e^{2\pi in(kT/N)/T}}{\sqrt{N}}=\frac{e^{2\pi ink/N}}{\sqrt{N}},$$
将纯音映射到数字纯音时，索引 $n$ 对应于频率 $\nu =n/T$，$N$ 对应于一个周期内采集的样本数。由于$T{f}_s=N$，其中$f_s$是采样频率，因此我们在频率和频率指数之间有以下联系：
$$\nu =\frac{n{f}_s}{N}\text{ and }n=\frac{\nu N}{f_s}$$
以下引理表明傅立叶基是正交的，因此它确实是一个基。

离散傅里叶变换

我们用$F_N$来表示坐标矩阵，从标准基${\mathbb{R}}^N$到傅里叶基${\mathcal{F}}_N$的变化。我们也将其称为（$N$​点）傅立叶矩阵。

矩阵 $\sqrt{N}{F}_N$ 也称为（$N$点）离散傅里叶变换，或 DFT。如果$\mathbf{x}$是$R^N$中的向量，则$\mathbf{y}=$ DFT $\mathbf{x}$称为$\mathbf{x}$的DFT系数。 （因此，DFT 系数是 ${\mathcal{F}}_N$ 中的坐标，用 $\sqrt{N}$ 缩放）。 DFT $\mathbf{x}$ 有时写为 $\hat{\mathbf{x}}$​。

请注意，我们将傅立叶矩阵和 DFT 定义为两个不同的矩阵，一个是另一个的缩放版本。究其原因，在于不同领域有不同的传统。在纯数学中，主要使用傅立叶矩阵，因为正如我们将看到的，它是酉矩阵。在信号处理中，主要使用 DFT 提供的缩放版本。我们通常将${\mathbb{R}}^N$中的给定向量写为$\mathbf{x}$，并为其DFT写为$\mathbf{y}$。在应用领域中，傅里叶基向量也称为合成向量，因为它们可以用来“合成”向量 $\mathbf{x}$，其权重由傅里叶基中的坐标提供，即
$$\mathbf{x}=y_0{\varphi }_0+y_1{\varphi }_1+\cdots +y_{N-1}{\varphi }_{N-1}\text{. }$$
此方程也称为合成方程。

Python离散傅里叶变换计算并绘制振幅谱

离散傅里叶变换可以将均匀间隔的信号序列变换为需要求和为时域信号的所有正弦波的频率信息。它定义为：
$$X_k=\sum _{n=0}^{N-1}{x}_n\cdot e^{-i2\pi kn/N}=\sum _{n=0}^{N-1}{x}_n[\cos (2\pi kn/N)-i\cdot \sin (2\pi kn/N)]$$
让我们看看如何使用它。

import matplotlib.pyplot as plt
import numpy as npplt.style.use('seaborn-poster')
%matplotlib inline# sampling rate
sr = 100
# sampling interval
ts = 1.0/sr
t = np.arange(0,1,ts)freq = 1.
x = 3*np.sin(2*np.pi*freq*t)freq = 4
x += np.sin(2*np.pi*freq*t)freq = 7   
x += 0.5* np.sin(2*np.pi*freq*t)plt.figure(figsize = (8, 6))
plt.plot(t, x, 'r')
plt.ylabel('Amplitude')plt.show()

编写一个函数 FT(x)，它接受一个参数，x - 输入一维实值信号。该函数将计算信号的 DFT 并返回 DFT 值。将此函数应用于我们上面生成的信号并绘制结果。

def FT(x):N = len(x)n = np.arange(N)k = n.reshape((N, 1))e = np.exp(-2j * np.pi * k * n / N)X = np.dot(e, x)return XX = FT(x)N = len(X)
n = np.arange(N)
T = N/sr
freq = n/T plt.figure(figsize = (8, 6))
plt.stem(freq, abs(X), 'b', \markerfmt=" ", basefmt="-b")
plt.xlabel('Freq (Hz)')
plt.ylabel('Amplitude |X(freq)|')
plt.show()

绘制 DFT 结果的前半部分，我们可以看到频率为 1 Hz、4 Hz 和 7 Hz 的 3 个清晰峰值，幅度为 3、1、0.5，符合预期。 这就是我们如何使用离散傅里叶变换将任意信号分解为简单的正弦波来分析它。

数字图像小波分析

小波是具有与傅里叶基不同性质的函数基，因此它们可以用来解决上述的一些缺点。 与傅立叶基相反，小波基不是固定的：存在多种此类基，用于不同的应用。

Python线性代数傅里叶分析和动态系统模拟分析之一

参阅一：计算思维

参阅二：亚图跨际