Python線性代數數字圖像和小波分析之二

要點

數學方程：數字信號和傅里葉分析，離散時間濾波器，小波分析
Python代碼實現及應用變換過程：
1. 讀取音頻和處理音頻波，使用Karplus-強算法制作吉他音頻
2. 離散傅里葉計算功能和繪制圖示結果
3. 計算波形傅里葉系數
4. 正向和反向（逆）快速傅里葉FFT 實現和使用
5. 位反轉函數
6. 正向和反向（逆）離散余弦變換
7. 離散時間濾波器結果繪圖
8. 一維、二維和三維正向和反向（逆）離散小波變換
9. 張量積將小波應用于圖像，二維和三維張量積實現

離散傅里葉分析數字聲音

數字聲音

數字聲音是一個序列 $\boldsymbol{x}=\left\{x_i\right\}_{i=0}^{N-1}$ ，對應于記錄的 $a$ 連續聲音 $f$ 的測量值以每秒 $f_s$ 測量的固定速率，即
$x_k=f\left(k / f_s\right), \quad \text { 對于 } k=0,1, \ldots, 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]$ = $A s in (2 π f n)$ ，其中 $A$ 是振幅， $f$ 是頻率。信號的 $y$ 軸是幅度或音量， $x$ 軸是當前時間。

讓我們嘗試使用以下方程生成正弦波。
$\sin (2 \pi f t[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)

方波會產生比正弦波更刺耳的嗡嗡聲。這是因為方波比只有一種音調的正弦波有更多的泛音。

離散傅里葉分析

在離散傅里葉分析中，向量 $\boldsymbol{x}=\left(x_0, \ldots, x_{N-1}\right)$ 表示為 $N$ 個向量的線性組合
$\phi_n=\frac{1}{\sqrt{N}}\left(1, e^{2 \pi i n / N}, e^{2 \pi i 2 n / N}, \ldots, e^{2 \pi i k n / N}, \ldots, e^{2 \pi i n(N-1) / N}\right)$
這些向量稱為歸一化復指數，或 $N$ 階純數字音調。 $n$ 也稱為頻率指數。整個集合 $\mathcal{F}_N=\left\{\phi_n\right\}_{n=0}^{N-1}$ 稱為 $N$ 點傅里葉基。

請注意，純數字音可以被視為純音的樣本，在一個周期內均勻采集： $f(t)=e^{2 \pi i n t / T} / \sqrt{N}$ 是具有頻率的純音 $n / T$ ，其樣本為
$N)=\frac{e^{2 \pi i n(k T / N) / T}}{\sqrt{N}}=\frac{e^{2 \pi i n k / 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。如果 $\boldsymbol{x}$ 是 $R^N$ 中的向量，則 $\boldsymbol{y}=$ DFT $\boldsymbol{x}$ 稱為 $\boldsymbol{x}$ 的DFT系數。（因此，DFT 系數是 $\mathcal{F}_N$ 中的坐標，用 $\sqrt{N}$ 縮放）。 DFT $\boldsymbol{x}$ 有時寫為 $\hat{\boldsymbol{x}}$ ?。

請注意，我們將傅立葉矩陣和 DFT 定義為兩個不同的矩陣，一個是另一個的縮放版本。究其原因，在于不同領域有不同的傳統。在純數學中，主要使用傅立葉矩陣，因為正如我們將看到的，它是酉矩陣。在信號處理中，主要使用 DFT 提供的縮放版本。我們通常將 $\mathbb{R}^N$ 中的給定向量寫為 $\boldsymbol{x}$ ，并為其DFT寫為 $\boldsymbol{y}$ 。在應用領域中，傅里葉基向量也稱為合成向量，因為它們可以用來“合成”向量 $\boldsymbol{x}$ ，其權重由傅里葉基中的坐標提供，即
$\boldsymbol{x}=y_0 \phi_0+y_1 \phi_1+\cdots+y_{N-1} \phi_{N-1} \text {. }$
此方程也稱為合成方程。

Python離散傅里葉變換計算并繪制振幅譜

離散傅里葉變換可以將均勻間隔的信號序列變換為需要求和為時域信號的所有正弦波的頻率信息。它定義為：
$X_k=\sum_{n=0}^{N-1} x_n \cdot e^{-i 2 \pi k n / N}=\sum_{n=0}^{N-1} x_n[\cos (2 \pi k n / N)-i \cdot \sin (2 \pi k n / 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，符合預期。這就是我們如何使用離散傅里葉變換將任意信號分解為簡單的正弦波來分析它。