訊號處理 (scipy.signal)

訊號處理工具箱目前包含一些濾波函數、一套有限的濾波器設計工具，以及一些用於 1D 和 2D 資料的 B-spline 內插演算法。雖然 B-spline 演算法在技術上可以歸類在內插類別下，但由於它們僅適用於等間隔資料，並且大量使用濾波器理論和轉移函數形式來提供快速的 B-spline 轉換，因此在此處包含它們。為了理解本節，您需要理解 SciPy 中的訊號是實數或複數的陣列。

B-splines

B-spline 是有限域上連續函數的近似，以 B-spline 係數和節點表示。如果節點間隔相等，間距為 \(\Delta x\)，則 1D 函數的 B-spline 近似是有限基底展開。

\[y\left(x\right)\approx\sum_{j}c_{j}\beta^{o}\left(\frac{x}{\Delta x}-j\right).\]

在二維中，節點間距為 \(\Delta x\) 和 \(\Delta y\)，函數表示為

\[z\left(x,y\right)\approx\sum_{j}\sum_{k}c_{jk}\beta^{o}\left(\frac{x}{\Delta x}-j\right)\beta^{o}\left(\frac{y}{\Delta y}-k\right).\]

在這些表達式中，\(\beta^{o}\left(\cdot\right)\) 是階數為 \(o\) 的空間有限 B-spline 基底函數。等間隔節點和等間隔資料點的要求，允許開發快速（逆濾波）演算法，以從樣本值 \(y_{n}\) 確定係數 \(c_{j}\)。與一般的 spline 內插演算法不同，這些演算法可以快速找到大型影像的 spline 係數。

透過 B-spline 基底函數表示一組樣本的優點是，可以從 spline 係數相對容易地計算連續域運算子（導數、重新取樣、積分等），這些運算子假設資料樣本是從底層連續函數中提取的。例如，spline 的二階導數是

\[y{}^{\prime\prime}\left(x\right)=\frac{1}{\Delta x^{2}}\sum_{j}c_{j}\beta^{o\prime\prime}\left(\frac{x}{\Delta x}-j\right).\]

使用 B-spline 的性質，即

\[\frac{d^{2}\beta^{o}\left(w\right)}{dw^{2}}=\beta^{o-2}\left(w+1\right)-2\beta^{o-2}\left(w\right)+\beta^{o-2}\left(w-1\right),\]

可以看出

\[y^{\prime\prime}\left(x\right)=\frac{1}{\Delta x^{2}}\sum_{j}c_{j}\left[\beta^{o-2}\left(\frac{x}{\Delta x}-j+1\right)-2\beta^{o-2}\left(\frac{x}{\Delta x}-j\right)+\beta^{o-2}\left(\frac{x}{\Delta x}-j-1\right)\right].\]

如果 \(o=3\)，則在樣本點

\begin{eqnarray*} \Delta x^{2}\left.y^{\prime}\left(x\right)\right|_{x=n\Delta x} & = & \sum_{j}c_{j}\delta_{n-j+1}-2c_{j}\delta_{n-j}+c_{j}\delta_{n-j-1},\\ & = & c_{n+1}-2c_{n}+c_{n-1}.\end{eqnarray*}

因此，可以從 spline 擬合中輕鬆計算二階導數訊號。如果需要，可以找到平滑 spline，以降低二階導數對隨機誤差的敏感度。

精明的讀者可能已經注意到，資料樣本透過卷積運算子與節點係數相關，因此簡單地與取樣的 B-spline 函數進行卷積，即可從 spline 係數中恢復原始資料。卷積的輸出可能會根據邊界處理方式而改變（隨著資料集中維度的數量增加，這變得越來越重要）。訊號處理子套件中與 B-spline 相關的演算法假設鏡像對稱邊界條件。因此，spline 係數是基於該假設計算的，並且可以透過假設它們也是鏡像對稱的，從 spline 係數中準確地恢復資料樣本。

目前，該套件提供函數，用於從一維和二維的等間隔樣本中確定二階和三階立方 spline 係數 (qspline1d, qspline2d, cspline1d, cspline2d)。對於大的 \(o\)，B-spline 基底函數可以很好地近似為零均值高斯函數，其標準差等於 \(\sigma_{o}=\left(o+1\right)/12\)

\[\beta^{o}\left(x\right)\approx\frac{1}{\sqrt{2\pi\sigma_{o}^{2}}}\exp\left(-\frac{x^{2}}{2\sigma_{o}}\right).\]

一個針對任意 \(x\) 和 \(o\) 計算此高斯函數的函數也可用 (gauss_spline )。以下程式碼和圖示使用 spline 濾波來計算浣熊臉部的邊緣影像（平滑 spline 的二階導數），這是一個由命令 scipy.datasets.face 返回的陣列。命令 sepfir2d 用於將可分離的 2D FIR 濾波器與鏡像對稱邊界條件應用於 spline 係數。此函數非常適合從 spline 係數重建樣本，並且比 convolve2d 更快，後者卷積任意 2D 濾波器並允許選擇鏡像對稱邊界條件。

>>> import numpy as np
>>> from scipy import signal, datasets
>>> import matplotlib.pyplot as plt

>>> image = datasets.face(gray=True).astype(np.float32)
>>> derfilt = np.array([1.0, -2, 1.0], dtype=np.float32)
>>> ck = signal.cspline2d(image, 8.0)
>>> deriv = (signal.sepfir2d(ck, derfilt, [1]) +
...          signal.sepfir2d(ck, [1], derfilt))

或者，我們可以執行

laplacian = np.array([[0,1,0], [1,-4,1], [0,1,0]], dtype=np.float32)
deriv2 = signal.convolve2d(ck,laplacian,mode='same',boundary='symm')

>>> plt.figure()
>>> plt.imshow(image)
>>> plt.gray()
>>> plt.title('Original image')
>>> plt.show()

>>> plt.figure()
>>> plt.imshow(deriv)
>>> plt.gray()
>>> plt.title('Output of spline edge filter')
>>> plt.show()

濾波

濾波是對輸入訊號進行某種修改的任何系統的通用名稱。在 SciPy 中，訊號可以被認為是 NumPy 陣列。有不同種類的濾波器用於不同種類的運算。有兩大類濾波運算：線性與非線性。線性濾波器始終可以簡化為將展平的 NumPy 陣列乘以適當的矩陣，從而產生另一個展平的 NumPy 陣列。當然，這通常不是計算濾波器的最佳方法，因為所涉及的矩陣和向量可能非常龐大。例如，使用此方法濾波 \(512 \times 512\) 影像將需要將 \(512^2 \times 512^2\) 矩陣與 \(512^2\) 向量相乘。僅嘗試使用標準 NumPy 陣列儲存 \(512^2 \times 512^2\) 矩陣就需要 \(68,719,476,736\) 個元素。每個元素 4 個位元組，這將需要 \(256\textrm{GB}\) 的記憶體。在大多數應用中，此矩陣的大多數元素為零，並且採用不同的方法來計算濾波器的輸出。

卷積/相關性

許多線性濾波器也具有位移不變性。這表示濾波操作在訊號中的不同位置是相同的，並且表示濾波矩陣可以僅從矩陣的一行（或一列）的知識建構出來。在這種情況下，可以使用傅立葉轉換完成矩陣乘法。

讓 \(x\left[n\right]\) 定義一個由整數 \(n.\) 索引的 1D 訊號。兩個 1D 訊號的完整卷積可以表示為

\[y\left[n\right]=\sum_{k=-\infty}^{\infty}x\left[k\right]h\left[n-k\right].\]

只有當我們將序列限制為可以儲存在電腦中的有限支援序列時，才能直接實作此方程式，選擇 \(n=0\) 作為兩個序列的起點，讓 \(K+1\) 為值，對於所有 \(n\geq K+1\)，\(x\left[n\right]=0\)，並且讓 \(M+1\) 為值，對於所有 \(n\geq M+1\)，\(h\left[n\right]=0\)，則離散卷積表達式為

\[y\left[n\right]=\sum_{k=\max\left(n-M,0\right)}^{\min\left(n,K\right)}x\left[k\right]h\left[n-k\right].\]

為了方便起見，假設 \(K\geq M.\) 那麼，更明確地說，此操作的輸出為

\begin{eqnarray*} y\left[0\right] & = & x\left[0\right]h\left[0\right]\\ y\left[1\right] & = & x\left[0\right]h\left[1\right]+x\left[1\right]h\left[0\right]\\ y\left[2\right] & = & x\left[0\right]h\left[2\right]+x\left[1\right]h\left[1\right]+x\left[2\right]h\left[0\right]\\ \vdots & \vdots & \vdots\\ y\left[M\right] & = & x\left[0\right]h\left[M\right]+x\left[1\right]h\left[M-1\right]+\cdots+x\left[M\right]h\left[0\right]\\ y\left[M+1\right] & = & x\left[1\right]h\left[M\right]+x\left[2\right]h\left[M-1\right]+\cdots+x\left[M+1\right]h\left[0\right]\\ \vdots & \vdots & \vdots\\ y\left[K\right] & = & x\left[K-M\right]h\left[M\right]+\cdots+x\left[K\right]h\left[0\right]\\ y\left[K+1\right] & = & x\left[K+1-M\right]h\left[M\right]+\cdots+x\left[K\right]h\left[1\right]\\ \vdots & \vdots & \vdots\\ y\left[K+M-1\right] & = & x\left[K-1\right]h\left[M\right]+x\left[K\right]h\left[M-1\right]\\ y\left[K+M\right] & = & x\left[K\right]h\left[M\right].\end{eqnarray*}

因此，長度分別為 \(K+1\) 和 \(M+1\) 的兩個有限序列的完整離散卷積，會產生長度為 \(K+M+1=\left(K+1\right)+\left(M+1\right)-1.\) 的有限序列。

1D 卷積在 SciPy 中使用函數 convolve 實作。此函數將訊號 \(x,\) \(h\) 和兩個可選標誌 ‘mode’ 和 ‘method’ 作為輸入，並傳回訊號 \(y.\)

第一個可選標誌 ‘mode’ 允許指定要傳回的輸出訊號的哪個部分。預設值 ‘full’ 傳回整個訊號。如果標誌的值為 ‘same’，則僅傳回中間 \(K\) 個值，從 \(y\left[\left\lfloor \frac{M-1}{2}\right\rfloor \right]\) 開始，以便輸出與第一個輸入具有相同的長度。如果標誌的值為 ‘valid’，則僅傳回中間 \(K-M+1=\left(K+1\right)-\left(M+1\right)+1\) 個輸出值，其中 \(z\) 取決於從 \(h\left[0\right]\) 到 \(h\left[M\right].\) 的最小輸入的所有值。換句話說，僅傳回從 \(y\left[M\right]\) 到 \(y\left[K\right]\)（包含）的值。

第二個可選標誌 ‘method’ 決定了卷積的計算方式，可以透過傅立葉轉換方法與 fftconvolve 或透過直接方法。預設情況下，它會選擇預期更快速的方法。傅立葉轉換方法的階數為 \(O(N\log N)\)，而直接方法的階數為 \(O(N^2)\)。根據大 O 常數和 \(N\) 的值，這兩種方法之一可能會更快。預設值 ‘auto’ 執行粗略計算並選擇預期更快速的方法，而值 ‘direct’ 和 ‘fft’ 強制使用其他兩種方法進行計算。

下面的程式碼顯示了 2 個序列卷積的簡單範例

>>> x = np.array([1.0, 2.0, 3.0])
>>> h = np.array([0.0, 1.0, 0.0, 0.0, 0.0])
>>> signal.convolve(x, h)
array([ 0.,  1.,  2.,  3.,  0.,  0.,  0.])
>>> signal.convolve(x, h, 'same')
array([ 2.,  3.,  0.])

這個相同的函數 convolve 實際上可以將 ND 陣列作為輸入，並將傳回兩個陣列的 ND 卷積，如下面的程式碼範例所示。相同的輸入標誌也適用於該情況。

>>> x = np.array([[1., 1., 0., 0.], [1., 1., 0., 0.], [0., 0., 0., 0.], [0., 0., 0., 0.]])
>>> h = np.array([[1., 0., 0., 0.], [0., 0., 0., 0.], [0., 0., 1., 0.], [0., 0., 0., 0.]])
>>> signal.convolve(x, h)
array([[ 1.,  1.,  0.,  0.,  0.,  0.,  0.],
       [ 1.,  1.,  0.,  0.,  0.,  0.,  0.],
       [ 0.,  0.,  1.,  1.,  0.,  0.,  0.],
       [ 0.,  0.,  1.,  1.,  0.,  0.,  0.],
       [ 0.,  0.,  0.,  0.,  0.,  0.,  0.],
       [ 0.,  0.,  0.,  0.,  0.,  0.,  0.],
       [ 0.,  0.,  0.,  0.,  0.,  0.,  0.]])

相關性與卷積非常相似，只是減號變為加號。因此，

\[w\left[n\right]=\sum_{k=-\infty}^{\infty}y\left[k\right]x\left[n+k\right],\]

是訊號 \(y\) 和 \(x.\) 的（交叉）相關性。對於有限長度訊號，其中 \(y\left[n\right]=0\) 在範圍 \(\left[0,K\right]\) 之外，並且 \(x\left[n\right]=0\) 在範圍 \(\left[0,M\right],\) 之外，總和可以簡化為

\[w\left[n\right]=\sum_{k=\max\left(0,-n\right)}^{\min\left(K,M-n\right)}y\left[k\right]x\left[n+k\right].\]

再次假設 \(K\geq M\)，則為

\begin{eqnarray*} w\left[-K\right] & = & y\left[K\right]x\left[0\right]\\ w\left[-K+1\right] & = & y\left[K-1\right]x\left[0\right]+y\left[K\right]x\left[1\right]\\ \vdots & \vdots & \vdots\\ w\left[M-K\right] & = & y\left[K-M\right]x\left[0\right]+y\left[K-M+1\right]x\left[1\right]+\cdots+y\left[K\right]x\left[M\right]\\ w\left[M-K+1\right] & = & y\left[K-M-1\right]x\left[0\right]+\cdots+y\left[K-1\right]x\left[M\right]\\ \vdots & \vdots & \vdots\\ w\left[-1\right] & = & y\left[1\right]x\left[0\right]+y\left[2\right]x\left[1\right]+\cdots+y\left[M+1\right]x\left[M\right]\\ w\left[0\right] & = & y\left[0\right]x\left[0\right]+y\left[1\right]x\left[1\right]+\cdots+y\left[M\right]x\left[M\right]\\ w\left[1\right] & = & y\left[0\right]x\left[1\right]+y\left[1\right]x\left[2\right]+\cdots+y\left[M-1\right]x\left[M\right]\\ w\left[2\right] & = & y\left[0\right]x\left[2\right]+y\left[1\right]x\left[3\right]+\cdots+y\left[M-2\right]x\left[M\right]\\ \vdots & \vdots & \vdots\\ w\left[M-1\right] & = & y\left[0\right]x\left[M-1\right]+y\left[1\right]x\left[M\right]\\ w\left[M\right] & = & y\left[0\right]x\left[M\right].\end{eqnarray*}

SciPy 函數 correlate 實作此運算。此運算提供等效的標誌，以傳回完整 \(K+M+1\) 長度序列 (‘full’) 或與最大序列大小相同的序列，從 \(w\left[-K+\left\lfloor \frac{M-1}{2}\right\rfloor \right]\) (‘same’) 開始，或者傳回一個序列，其中值取決於最小序列 (‘valid’) 的所有值。最後一個選項傳回 \(K-M+1\) 個值 \(w\left[M-K\right]\) 到 \(w\left[0\right]\)（包含）。

函數 correlate 也可以將任意 ND 陣列作為輸入，並在輸出時傳回兩個陣列的 ND 卷積。

當 \(N=2,\) 時，correlate 和/或 convolve 可用於建構任意影像濾波器，以執行影像的模糊化、增強和邊緣偵測等動作。

>>> import numpy as np
>>> from scipy import signal, datasets
>>> import matplotlib.pyplot as plt

>>> image = datasets.face(gray=True)
>>> w = np.zeros((50, 50))
>>> w[0][0] = 1.0
>>> w[49][25] = 1.0
>>> image_new = signal.fftconvolve(image, w)

>>> plt.figure()
>>> plt.imshow(image)
>>> plt.gray()
>>> plt.title('Original image')
>>> plt.show()

>>> plt.figure()
>>> plt.imshow(image_new)
>>> plt.gray()
>>> plt.title('Filtered image')
>>> plt.show()

如上所述，在時域中計算卷積主要用於濾波，當其中一個訊號遠小於另一個訊號時 ( \(K\gg M\) )，否則線性濾波在頻域中計算效率更高，由函數 fftconvolve 提供。預設情況下，convolve 使用 choose_conv_method 估計最快的方法。

如果濾波器函數 \(w[n,m]\) 可以根據以下方式分解

\[h[n, m] = h_1[n] h_2[m],\]

則可以使用函數 sepfir2d 計算卷積。例如，我們考慮高斯濾波器 gaussian

\[h[n, m] \propto e^{-x^2-y^2} = e^{-x^2} e^{-y^2},\]

通常用於模糊化。

>>> import numpy as np
>>> from scipy import signal, datasets
>>> import matplotlib.pyplot as plt

>>> image = np.asarray(datasets.ascent(), np.float64)
>>> w = signal.windows.gaussian(51, 10.0)
>>> image_new = signal.sepfir2d(image, w, w)

>>> plt.figure()
>>> plt.imshow(image)
>>> plt.gray()
>>> plt.title('Original image')
>>> plt.show()

>>> plt.figure()
>>> plt.imshow(image_new)
>>> plt.gray()
>>> plt.title('Filtered image')
>>> plt.show()

差分方程式濾波

線性 1D 濾波器（包含卷積濾波器）的一般類別是由差分方程式描述的濾波器

\[\sum_{k=0}^{N}a_{k}y\left[n-k\right]=\sum_{k=0}^{M}b_{k}x\left[n-k\right],\]

其中 \(x\left[n\right]\) 是輸入序列，\(y\left[n\right]\) 是輸出序列。如果我們假設初始靜止狀態，使得 \(y\left[n\right]=0\) 對於 \(n<0\)，則可以使用卷積實作此類型的濾波器。但是，如果 \(a_{k}\neq0\) 對於 \(k\geq1.\)，則卷積濾波器序列 \(h\left[n\right]\) 可能是無限的。此外，此一般類別的線性濾波器允許將初始條件置於 \(y\left[n\right]\) 上，對於 \(n<0\)，導致無法使用卷積表示的濾波器。

差分方程式濾波器可以認為是根據其先前值遞迴地找到 \(y\left[n\right]\)

\[a_{0}y\left[n\right]=-a_{1}y\left[n-1\right]-\cdots-a_{N}y\left[n-N\right]+\cdots+b_{0}x\left[n\right]+\cdots+b_{M}x\left[n-M\right].\]

通常，選擇 \(a_{0}=1\) 進行正規化。SciPy 中此一般差分方程式濾波器的實作比先前的方程式所暗示的要稍微複雜一些。它的實作方式是只需要延遲一個訊號。實際的實作方程式為（假設 \(a_{0}=1\) ）

\begin{eqnarray*} y\left[n\right] & = & b_{0}x\left[n\right]+z_{0}\left[n-1\right]\\ z_{0}\left[n\right] & = & b_{1}x\left[n\right]+z_{1}\left[n-1\right]-a_{1}y\left[n\right]\\ z_{1}\left[n\right] & = & b_{2}x\left[n\right]+z_{2}\left[n-1\right]-a_{2}y\left[n\right]\\ \vdots & \vdots & \vdots\\ z_{K-2}\left[n\right] & = & b_{K-1}x\left[n\right]+z_{K-1}\left[n-1\right]-a_{K-1}y\left[n\right]\\ z_{K-1}\left[n\right] & = & b_{K}x\left[n\right]-a_{K}y\left[n\right],\end{eqnarray*}

其中 \(K=\max\left(N,M\right).\) 請注意，如果 \(K>M\) 則 \(b_{K}=0\)，如果 \(K>N.\) 則 \(a_{K}=0\)。透過這種方式，時間 \(n\) 的輸出僅取決於時間 \(n\) 的輸入和前一時間的 \(z_{0}\) 值。只要在每個時間步驟計算並儲存 \(K\) 個值 \(z_{0}\left[n-1\right]\ldots z_{K-1}\left[n-1\right]\)，就可以始終計算出此值。

差分方程式濾波器使用 SciPy 中的命令 lfilter 呼叫。此命令將向量 \(b,\) 向量 \(a,\) 訊號 \(x\) 作為輸入，並傳回使用上述方程式計算的向量 \(y\) （與 \(x\) 長度相同）。如果 \(x\) 是 ND，則濾波器沿提供的軸計算。如果需要，可以提供提供 \(z_{0}\left[-1\right]\) 到 \(z_{K-1}\left[-1\right]\) 值的初始條件，否則將假定它們都為零。如果提供了初始條件，則也會傳回中間變數的最終條件。例如，這些可用於在相同狀態下重新啟動計算。

有時，以訊號 \(x\left[n\right]\) 和 \(y\left[n\right].\) 表示初始條件更方便。換句話說，也許您有 \(x\left[-M\right]\) 到 \(x\left[-1\right]\) 的值，以及 \(y\left[-N\right]\) 到 \(y\left[-1\right]\) 的值，並且想要確定應將哪些 \(z_{m}\left[-1\right]\) 值作為初始條件傳遞到差分方程式濾波器。不難證明，對於 \(0\leq m<K,\)

\[z_{m}\left[n\right]=\sum_{p=0}^{K-m-1}\left(b_{m+p+1}x\left[n-p\right]-a_{m+p+1}y\left[n-p\right]\right).\]

使用此公式，我們可以找到初始條件向量 \(z_{0}\left[-1\right]\) 到 \(z_{K-1}\left[-1\right]\)，給定 \(y\) （和 \(x\) ）的初始條件。命令 lfiltic 執行此功能。

例如，考慮以下系統

\[y[n] = \frac{1}{2} x[n] + \frac{1}{4} x[n-1] + \frac{1}{3} y[n-1]\]

程式碼針對給定的訊號 \(x[n]\) 計算訊號 \(y[n]\)；首先針對初始條件 \(y[-1] = 0\) （預設情況），然後針對 \(y[-1] = 2\)，透過 lfiltic。

>>> import numpy as np
>>> from scipy import signal

>>> x = np.array([1., 0., 0., 0.])
>>> b = np.array([1.0/2, 1.0/4])
>>> a = np.array([1.0, -1.0/3])
>>> signal.lfilter(b, a, x)
array([0.5, 0.41666667, 0.13888889, 0.0462963])
>>> zi = signal.lfiltic(b, a, y=[2.])
>>> signal.lfilter(b, a, x, zi=zi)
(array([ 1.16666667,  0.63888889,  0.21296296,  0.07098765]), array([0.02366]))

請注意，輸出訊號 \(y[n]\) 的長度與輸入訊號 \(x[n]\) 的長度相同。

線性系統分析

線性系統由線性差分方程式描述，可以完全由係數向量 \(a\) 和 \(b\) 描述，如上所述；另一種表示形式是提供因子 \(k\)、\(N_z\) 個零點 \(z_k\) 和 \(N_p\) 個極點 \(p_k\)，分別根據其轉移函數 \(H(z)\) 描述系統，如下所示

\[H(z) = k \frac{ (z-z_1)(z-z_2)...(z-z_{N_z})}{ (z-p_1)(z-p_2)...(z-p_{N_p})}.\]

這種替代表示形式可以使用 scipy 函數 tf2zpk 獲得；反向操作由 zpk2tf 提供。

對於上面的範例，我們有

>>> b = np.array([1.0/2, 1.0/4])
>>> a = np.array([1.0, -1.0/3])
>>> signal.tf2zpk(b, a)
(array([-0.5]), array([ 0.33333333]), 0.5)

也就是說，系統在 \(z=-1/2\) 處有一個零點，在 \(z=1/3\) 處有一個極點。

scipy 函數 freqz 允許計算由係數 \(a_k\) 和 \(b_k\) 描述的系統的頻率響應。有關綜合範例，請參閱 freqz 函數的說明。

濾波器設計

時域離散濾波器可以分為有限響應 (FIR) 濾波器和無限響應 (IIR) 濾波器。FIR 濾波器可以提供線性相位響應，而 IIR 濾波器則不能。SciPy 提供了用於設計這兩種濾波器的函數。

FIR 濾波器

函數 firwin 根據視窗方法設計濾波器。根據提供的引數，該函數傳回不同的濾波器類型（例如，低通、帶通…）。

下面的範例分別設計了低通和帶止濾波器。

>>> import numpy as np
>>> import scipy.signal as signal
>>> import matplotlib.pyplot as plt

>>> b1 = signal.firwin(40, 0.5)
>>> b2 = signal.firwin(41, [0.3, 0.8])
>>> w1, h1 = signal.freqz(b1)
>>> w2, h2 = signal.freqz(b2)

>>> plt.title('Digital filter frequency response')
>>> plt.plot(w1, 20*np.log10(np.abs(h1)), 'b')
>>> plt.plot(w2, 20*np.log10(np.abs(h2)), 'r')
>>> plt.ylabel('Amplitude Response (dB)')
>>> plt.xlabel('Frequency (rad/sample)')
>>> plt.grid()
>>> plt.show()

請注意，firwin 預設使用正規化頻率，定義為值 \(1\) 對應於奈奎斯特頻率，而函數 freqz 的定義是值 \(\pi\) 對應於奈奎斯特頻率。

函數 firwin2 允許透過指定角頻率陣列和對應的增益，設計幾乎任意的頻率響應。

下面的範例設計了一個具有此類任意振幅響應的濾波器。

>>> import numpy as np
>>> import scipy.signal as signal
>>> import matplotlib.pyplot as plt

>>> b = signal.firwin2(150, [0.0, 0.3, 0.6, 1.0], [1.0, 2.0, 0.5, 0.0])
>>> w, h = signal.freqz(b)

>>> plt.title('Digital filter frequency response')
>>> plt.plot(w, np.abs(h))
>>> plt.title('Digital filter frequency response')
>>> plt.ylabel('Amplitude Response')
>>> plt.xlabel('Frequency (rad/sample)')
>>> plt.grid()
>>> plt.show()

請注意 y 軸的線性縮放以及 firwin2 和 freqz 中奈奎斯特頻率的不同定義（如上所述）。

IIR 濾波器

SciPy 提供了兩個函數來直接設計 IIR，iirdesign 和 iirfilter，其中濾波器類型（例如，橢圓）作為引數傳遞，以及幾個用於特定濾波器類型的更多濾波器設計函數，例如，ellip。

下面的範例設計了一個橢圓低通濾波器，分別具有定義的通帶和阻帶漣波。請注意，為了達到與上述範例中 FIR 濾波器相同的 \(\approx 60\) dB 阻帶衰減，濾波器階數要低得多（4 階）。

>>> import numpy as np
>>> import scipy.signal as signal
>>> import matplotlib.pyplot as plt

>>> b, a = signal.iirfilter(4, Wn=0.2, rp=5, rs=60, btype='lowpass', ftype='ellip')
>>> w, h = signal.freqz(b, a)

>>> plt.title('Digital filter frequency response')
>>> plt.plot(w, 20*np.log10(np.abs(h)))
>>> plt.title('Digital filter frequency response')
>>> plt.ylabel('Amplitude Response [dB]')
>>> plt.xlabel('Frequency (rad/sample)')
>>> plt.grid()
>>> plt.show()

注意

請注意，firwin 和 iirfilter 的截止頻率定義方式不同。對於 firwin，截止頻率位於半振幅處 (即 -6dB)。對於 iirfilter，截止頻率位於半功率處 (即 -3dB)。

>>> import numpy as np
>>> from matplotlib import pyplot as plt
>>> from scipy import signal as sig

>>> fs = 16000
>>> b = sig.firwin(101, 2500, fs=fs)
>>> f, h_fft = sig.freqz(b, fs=fs)
>>> h_amp = 20 * np.log10(np.abs(h_fft))
>>> _, ax = plt.subplots(layout="constrained")
>>> ax.plot(f, h_amp, label="FIR")
>>> ax.grid(True)

>>> b, a = sig.iirfilter(15, 2500, btype="low", fs=fs)
>>> f, h_fft = sig.freqz(b, a, fs=fs)
>>> h_amp = 20 * np.log10(np.abs(h_fft))
>>> ax.plot(f, h_amp, label="IIR")
>>> ax.set(xlim=[2100, 2900], ylim=[-10, 2])
>>> ax.set(xlabel="Frequency (Hz)", ylabel="Amplitude Response [dB]")
>>> ax.legend()

濾波器係數

濾波器係數可以儲存成幾種不同的格式

• ‘ba’ 或 ‘tf’ = 轉移函數係數

• ‘zpk’ = 零點、極點和整體增益

• ‘ss’ = 狀態空間系統表示法

• ‘sos’ = 二階sections的轉移函數係數

諸如 tf2zpk 和 zpk2ss 等函數可以在這些格式之間轉換。

轉移函數表示法

ba 或 tf 格式是一個 2-tuple (b, a)，代表一個轉移函數，其中 b 是一個長度為 M+1 的陣列，包含 M 階分子多項式的係數，而 a 是一個長度為 N+1 的陣列，包含 N 階分母的係數，以轉移函數變數的正向遞減冪次排列。因此，\(b = [b_0, b_1, ..., b_M]\) 和 \(a =[a_0, a_1, ..., a_N]\) 的 tuple 可以表示以下形式的類比濾波器

\[H(s) = \frac {b_0 s^M + b_1 s^{(M-1)} + \cdots + b_M} {a_0 s^N + a_1 s^{(N-1)} + \cdots + a_N} = \frac {\sum_{i=0}^M b_i s^{(M-i)}} {\sum_{i=0}^N a_i s^{(N-i)}}\]

或以下形式的離散時間濾波器

\[H(z) = \frac {b_0 z^M + b_1 z^{(M-1)} + \cdots + b_M} {a_0 z^N + a_1 z^{(N-1)} + \cdots + a_N} = \frac {\sum_{i=0}^M b_i z^{(M-i)}} {\sum_{i=0}^N a_i z^{(N-i)}}.\]

這種「正冪次」形式在控制工程中更常見。如果 M 和 N 相等 (對於雙線性轉換產生之所有濾波器皆為真)，則這恰好等同於 DSP 中偏好的「負冪次」離散時間形式

\[H(z) = \frac {b_0 + b_1 z^{-1} + \cdots + b_M z^{-M}} {a_0 + a_1 z^{-1} + \cdots + a_N z^{-N}} = \frac {\sum_{i=0}^M b_i z^{-i}} {\sum_{i=0}^N a_i z^{-i}}.\]

雖然這對於常見的濾波器來說是正確的，但請記住，這在一般情況下並不成立。如果 M 和 N 不相等，則在找到極點和零點之前，必須先將離散時間轉移函數係數轉換為「正冪次」形式。

這種表示法在高階時會產生數值誤差，因此在可能的情況下，最好使用其他格式。

零點和極點表示法

zpk 格式是一個 3-tuple (z, p, k)，其中 z 是一個長度為 M 的陣列，包含轉移函數的複數零點 \(z = [z_0, z_1, ..., z_{M-1}]\)，p 是一個長度為 N 的陣列，包含轉移函數的複數極點 \(p = [p_0, p_1, ..., p_{N-1}]\)，而 k 是一個純量增益。這些代表數位轉移函數

\[H(z) = k \cdot \frac {(z - z_0) (z - z_1) \cdots (z - z_{(M-1)})} {(z - p_0) (z - p_1) \cdots (z - p_{(N-1)})} = k \frac {\prod_{i=0}^{M-1} (z - z_i)} {\prod_{i=0}^{N-1} (z - p_i)}\]

或類比轉移函數

\[H(s) = k \cdot \frac {(s - z_0) (s - z_1) \cdots (s - z_{(M-1)})} {(s - p_0) (s - p_1) \cdots (s - p_{(N-1)})} = k \frac {\prod_{i=0}^{M-1} (s - z_i)} {\prod_{i=0}^{N-1} (s - p_i)}.\]

雖然根集合以排序的 NumPy 陣列儲存，但它們的順序並不重要：([-1, -2], [-3, -4], 1) 與 ([-2, -1], [-4, -3], 1) 是相同的濾波器。

狀態空間系統表示法

ss 格式是一個包含 4 個陣列的 tuple (A, B, C, D)，代表一個 N 階數位/離散時間系統的狀態空間，其形式為

\[\begin{split}\mathbf{x}[k+1] = A \mathbf{x}[k] + B \mathbf{u}[k]\\ \mathbf{y}[k] = C \mathbf{x}[k] + D \mathbf{u}[k]\end{split}\]

或連續/類比系統的形式

\[\begin{split}\dot{\mathbf{x}}(t) = A \mathbf{x}(t) + B \mathbf{u}(t)\\ \mathbf{y}(t) = C \mathbf{x}(t) + D \mathbf{u}(t),\end{split}\]

其中有 P 個輸入、 Q 個輸出和 N 個狀態變數，其中

• x 是狀態向量

• y 是長度為 Q 的輸出向量

• u 是長度為 P 的輸入向量

• A 是狀態矩陣，形狀為 (N, N)

• B 是輸入矩陣，形狀為 (N, P)

• C 是輸出矩陣，形狀為 (Q, N)

• D 是直通或前饋矩陣，形狀為 (Q, P)。（在系統沒有直接貫穿的情況下，D 中的所有值均為零。）

狀態空間是最通用的表示法，也是唯一允許用於多輸入、多輸出 (MIMO) 系統的表示法。對於給定的轉移函數，有多種狀態空間表示法。具體來說，「可控制標準型」和「可觀測標準型」與 tf 表示法具有相同的係數，因此也遭受相同的數值誤差。

二階sections表示法

sos 格式是一個單一的 2-D 陣列，形狀為 (n_sections, 6)，代表一系列二階轉移函數，當它們串聯級聯時，可以實現具有最小數值誤差的高階濾波器。每行對應於一個二階 tf 表示法，前三列提供分子係數，後三列提供分母係數

\[[b_0, b_1, b_2, a_0, a_1, a_2]\]

這些係數通常會被標準化，使得 \(a_0\) 始終為 1。對於浮點運算，section的順序通常並不重要；無論順序如何，濾波器輸出都將相同。

濾波器轉換

IIR 濾波器設計函數首先生成一個原型類比低通濾波器，其標準化截止頻率為 1 rad/sec。然後使用以下替換將其轉換為其他頻率和頻帶類型

類型	轉換
lp2lp	\(s \rightarrow \frac{s}{\omega_0}\)
lp2hp	\(s \rightarrow \frac{\omega_0}{s}\)
lp2bp	\(s \rightarrow \frac{s^2 + {\omega_0}^2}{s \cdot \mathrm{BW}}\)
lp2bs	\(s \rightarrow \frac{s \cdot \mathrm{BW}}{s^2 + {\omega_0}^2}\)

在此，\(\omega_0\) 是新的截止或中心頻率，而 \(\mathrm{BW}\) 是頻寬。這些在對數頻率軸上保持對稱性。

為了將轉換後的類比濾波器轉換為數位濾波器，使用了 bilinear 轉換，它進行以下替換

\[s \rightarrow \frac{2}{T} \frac{z - 1}{z + 1},\]

其中 T 是取樣時間（取樣頻率的倒數）。

其他濾波器

訊號處理套件也提供了更多濾波器。

中值濾波器

當中值濾波器雜訊明顯地是非高斯雜訊，或需要保留邊緣時，通常會應用中值濾波器。中值濾波器的工作原理是對感興趣點周圍矩形區域中的所有陣列像素值進行排序。此鄰域像素值列表的樣本中值用作輸出陣列的值。樣本中值是排序後的鄰域值列表中的中間陣列值。如果鄰域中存在偶數個元素，則將中間兩個值的平均值用作中值。適用於 N-D 陣列的通用中值濾波器是 medfilt。僅適用於 2-D 陣列的特殊版本可作為 medfilt2d 使用。

順序濾波器

中值濾波器是更通用的濾波器類別（稱為順序濾波器）的特定範例。為了計算特定像素的輸出，所有順序濾波器都使用該像素周圍區域中的陣列值。這些陣列值經過排序，然後選擇其中一個值作為輸出值。對於中值濾波器，陣列值列表的樣本中值用作輸出。通用順序濾波器允許使用者選擇排序後的值中的哪一個將用作輸出。因此，例如，可以選擇選取列表中的最大值或最小值。順序濾波器除了輸入陣列和區域遮罩之外，還需要一個額外引數，用於指定應將排序後的相鄰陣列值列表中的哪個元素用作輸出。執行順序濾波器的命令是 order_filter。

維納濾波器

維納濾波器是一種用於影像去噪的簡單去模糊濾波器。這不是影像重建問題中常用的維納濾波器，而是一種簡單的局部均值濾波器。令 \(x\) 為輸入訊號，則輸出為

\[\begin{split}y=\left\{ \begin{array}{cc} \frac{\sigma^{2}}{\sigma_{x}^{2}}m_{x}+\left(1-\frac{\sigma^{2}}{\sigma_{x}^{2}}\right)x & \sigma_{x}^{2}\geq\sigma^{2},\\ m_{x} & \sigma_{x}^{2}<\sigma^{2},\end{array}\right.\end{split}\]

其中 \(m_{x}\) 是均值的局部估計值，而 \(\sigma_{x}^{2}\) 是變異數的局部估計值。這些估計值的視窗是一個可選的輸入參數（預設值為 \(3\times3\) ）。參數 \(\sigma^{2}\) 是一個閾值雜訊參數。如果未給出 \(\sigma\)，則將其估計為局部變異數的平均值。

希爾伯特濾波器

希爾伯特轉換從實數訊號建構複數值的解析訊號。例如，如果 \(x=\cos\omega n\)，則 \(y=\textrm{hilbert}\left(x\right)\) 將傳回 (邊緣附近除外) \(y=\exp\left(j\omega n\right).\) 在頻域中，希爾伯特轉換執行

\[Y=X\cdot H,\]

其中 \(H\) 對於正頻率為 \(2\)，對於負頻率為 \(0\)，對於零頻率為 \(1\)。

類比濾波器設計

函數 iirdesign、 iirfilter 和特定濾波器類型的濾波器設計函數 (例如， ellip) 都具有 analog 旗標，該旗標也允許類比濾波器的設計。

下面的範例設計了一個類比 (IIR) 濾波器，通過 tf2zpk 獲得極點和零點，並將它們繪製在複數 s 平面上。在振幅響應中可以清楚地看到 \(\omega \approx 150\) 和 \(\omega \approx 300\) 的零點。

>>> import numpy as np
>>> import scipy.signal as signal
>>> import matplotlib.pyplot as plt

>>> b, a = signal.iirdesign(wp=100, ws=200, gpass=2.0, gstop=40., analog=True)
>>> w, h = signal.freqs(b, a)

>>> plt.title('Analog filter frequency response')
>>> plt.plot(w, 20*np.log10(np.abs(h)))
>>> plt.ylabel('Amplitude Response [dB]')
>>> plt.xlabel('Frequency')
>>> plt.grid()
>>> plt.show()

>>> z, p, k = signal.tf2zpk(b, a)

>>> plt.plot(np.real(z), np.imag(z), 'ob', markerfacecolor='none')
>>> plt.plot(np.real(p), np.imag(p), 'xr')
>>> plt.legend(['Zeros', 'Poles'], loc=2)

>>> plt.title('Pole / Zero Plot')
>>> plt.xlabel('Real')
>>> plt.ylabel('Imaginary')
>>> plt.grid()
>>> plt.show()

\[% LaTeX Macros to make the LaTeX formulas more readable: \newcommand{\IC}{{\mathbb{C}}} % set of complex numbers \newcommand{\IN}{{\mathbb{N}}} % set of natural numbers \newcommand{\IR}{{\mathbb{R}}} % set of real numbers \newcommand{\IZ}{{\mathbb{Z}}} % set of integers \newcommand{\jj}{{\mathbb{j}}} % imaginary unit \newcommand{\e}{\operatorname{e}} % Euler's number \newcommand{\dd}{\operatorname{d}} % infinitesimal operator \newcommand{\abs}[1]{\left|#1\right|} % absolute value \newcommand{\conj}[1]{\overline{#1}} % complex conjugate \newcommand{\conjT}[1]{\overline{#1^T}} % transposed complex conjugate \newcommand{\inv}[1]{\left(#1\right)^{\!-1}} % inverse % Since the physics package is not loaded, we define the macros ourselves: \newcommand{\vb}[1]{\mathbf{#1}} % vectors and matrices are bold % new macros: \newcommand{\rect}{\operatorname{rect}} % rect or boxcar function \newcommand{\sinc}{\operatorname{sinc}} % sinc(t) := sin(pi*t) / (pi*t)\]

頻譜分析

頻譜分析是指研究訊號的傅立葉轉換 [1]。根據上下文，對於傅立葉轉換的各種頻譜表示，存在各種名稱，例如頻譜、頻譜密度或週期圖。[2] 本節通過固定持續時間的連續時間正弦波訊號範例來說明最常見的表示形式。然後討論在該正弦波的取樣版本上使用離散傅立葉轉換 [3]。

單獨的小節專門介紹頻譜的相位，估計功率頻譜密度，不使用 (periodogram) 和使用平均 (welch) 以及非均勻間隔訊號 (lombscargle)。

請注意，傅立葉級數的概念密切相關，但在一個關鍵點上有所不同：傅立葉級數具有由離散頻率諧波組成的頻譜，而在本節中，頻譜在頻率上是連續的。

連續時間正弦訊號

考慮一個振幅為 \(a\)、頻率為 \(f_x\) 和持續時間為 \(\tau\) 的正弦訊號，即，

(1)\[x(t) = a \sin(2 \pi f_x t)\, \rect(\frac{t}{\tau}-\frac{1}{2}) = \left(\frac{a}{2\jj} \e^{\jj 2 \pi f_x t} - \frac{a}{2\jj} \e^{-\jj 2 \pi f_x t} \right) \rect(\frac{t}{\tau}-\frac{1}{2})\ .\]

由於 \(\rect(t)\) 函數在 \(|t|<1/2\) 時為 1，在 \(|t|>1/2\) 時為零，因此它將 \(x(t)\) 限制在間隔 \([0, \tau]\) 內。用複數指數表示正弦波顯示了它的兩個週期性成分，頻率為 \(\pm f_x\)。我們假設 \(x(t)\) 是電壓訊號，因此它的單位為 \(\text{V}\)。

在訊號處理中，絕對平方的積分 \(|x(t)|^2\) 用於定義訊號的能量和功率，即，

(2)\[E_x := \int_0^\tau |x(t)|^2 \dd t\ = \frac{1}{2}|a|^2\tau\ , \qquad P_x := \frac{1}{\tau}E_x = \frac{1}{2}|a|^2\ .\]

功率 \(P_x\) 可以解釋為單位時間間隔的能量 \(E_x\)。單位方面，對 \(t\) 積分會導致乘以秒數。因此，\(E_x\) 的單位為 \(\text{V}^2\text{s}\)，而 \(P_x\) 的單位為 \(\text{V}^2\)。

將傅立葉轉換應用於 \(x(t)\)，即，

(3)\[X(f) = \int_{\IR} x(t)\, \e^{-2\jj\pi f t}\, \dd t = \frac{a \tau}{2\jj} \Big( \sinc\!\big(\tau (f-f_x)\big) - \sinc\!\big(\tau (f+f_x)\big) \Big)\e^{-\jj\pi\tau f}\ ,\]

會產生兩個以 \(\pm f_x\) 為中心的 \(\sinc(f) := \sin(\pi f) /(\pi f)\) 函數。幅度（絕對值）\(|X(f)|\) 在 \(\pm f_x\) 處有兩個最大值，值為 \(|a|\tau/2\)。在下面的圖中可以看到，\(X(f)\) 並不集中在 \(\pm f_x\) 的主瓣周圍，而是包含高度與 \(1/(\tau f)\) 成比例遞減的旁瓣。這種所謂的「頻譜洩漏」[4] 是由將正弦波限制在有限間隔內引起的。請注意，訊號持續時間 \(\tau\) 越短，洩漏越高。為了獨立於訊號持續時間，可以使用所謂的「振幅頻譜」\(X(f)/\tau\) 而不是頻譜 \(X(f)\)。它在 \(f\) 處的值對應於複數指數 \(\exp(\jj2\pi f t)\) 的振幅。

由於帕塞瓦爾定理，能量也可以通過其傅立葉轉換 \(X(f)\) 計算得出

(4)\[ E_X := \int_\IR \abs{X(f)}^2 \dd f = E_x\]

也是如此。例如，可以通過直接計算證明，方程式 (4) 中 \(X(f)\) 的能量為 \(|a|^2\tau/2\)。因此，頻帶 \([f_a, f_b]\) 中的訊號功率可以使用以下公式確定

(5)\[P_X^{a,b} = \frac{1}{\tau} \int_{f_a}^{f_b} \abs{X(f)}^2 \dd f\ .\]

因此，函數 \(|X(f)|^2\) 可以定義為所謂的「能量頻譜密度」，而 \(S_{xx}(f) := |X(f)|^2 / \tau\) 可以定義為 \(x(t)\) 的「功率頻譜密度」(PSD)。除了 PSD 之外，還使用了所謂的「振幅頻譜密度」\(X(f) / \sqrt{\tau}\)，它仍然包含相位資訊。它的絕對平方是 PSD，因此它與訊號的均方根 (RMS) 值 \(\sqrt{P_x}\) 的概念密切相關。

總之，本小節介紹了五種表示頻譜的方法

正弦訊號 \(x(t)\) (方程式 (1)) 的頻譜表示形式比較，單位為 \(\text{V}\)：

	頻譜	振幅頻譜	能量頻譜密度	功率頻譜密度 (PSD)	振幅頻譜密度
定義	\(X(f)\)	\(X(f) / \tau\)	\(|X(f)|^2\)	\(|X(f)|^2 / \tau\)	\(X(f) / \sqrt{\tau}\)
在 \(\pm f_x\) 處的幅度	\(\frac{1}{2}|a|\tau\)	\(\frac{1}{2}|a|\)	\(\frac{1}{4}|a|^2\tau^2\)	\(\frac{1}{4}|a|^2\tau\)	\(\frac{1}{2}|a|\sqrt{\tau}\)
單位	\(\text{V} / \text{Hz}\)	\(\text{V}\)	\(\text{V}^2\text{s} / \text{Hz}\)	\(\text{V}^2 / \text{Hz}\)	\(\text{V} / \sqrt{\text{Hz}}\)

請注意，上表所示的單位並非明確，例如，\(\text{V}^2\text{s} / \text{Hz} = \text{V}^2\text{s}^2 = \text{V}^2/ \text{Hz}^2\)。當使用振幅頻譜的 \(|X(f) / \tau|\) 絕對值時，它被稱為幅度頻譜。此外，請注意，表示形式的命名方案並不一致，並且在文獻中有所不同。

對於實數值訊號，通常使用所謂的「單邊」頻譜表示形式。它僅使用非負頻率（由於 \(X(-f)= \conj{X}(f)\) 如果 \(x(t)\in\IR\)）。有時，負頻率的值會添加到它們的正頻率對應項中。然後，振幅頻譜允許讀取 \(x(t)\) 在 \(f_x\) 處的完整（非一半）振幅正弦波，並且 PSD 中間隔的面積代表其完整（非一半）功率。請注意，對於振幅頻譜密度，正值不會加倍，而是乘以 \(\sqrt{2}\)，因為它是 PSD 的平方根。此外，沒有用於命名加倍頻譜的規範方法。

下圖顯示了四個正弦訊號 \(x(t)\) (方程式 (1)) 的三種不同頻譜表示形式，它們具有不同的振幅 \(a\) 和持續時間 \(\tau\)。為了減少混亂，頻譜以 \(f_x\) 為中心，並彼此相鄰繪製

import matplotlib.pyplot as plt
import numpy as np

aa = [1, 1, 2, 2]  # amplitudes
taus = [1, 2, 1, 2]  # durations

fg0, axx = plt.subplots(3, 1, sharex='all', tight_layout=True, figsize=(6., 4.))
axx[0].set(title=r"Spectrum $|X(f)|$", ylabel="V/Hz")
axx[1].set(title=r"Magnitude Spectrum $|X(f)/\tau|$ ", ylabel=r"V")
axx[2].set(title=r"Amplitude Spectral Density $|X(f)/\sqrt{\tau}|$",
           ylabel=r"$\operatorname{V} / \sqrt{\operatorname{Hz}}$",
           xlabel="Frequency $f$ in Hertz",)

x_labels, x_ticks = [], []
f = np.linspace(-2.5, 2.5, 400)
for c_, (a_, tau_) in enumerate(zip(aa, taus), start=1):
    aZ_, f_ = abs(a_ * tau_ * np.sinc(tau_ * f) / 2), f + c_ * 5
    axx[0].plot(f_, aZ_)
    axx[1].plot(f_, aZ_ / tau_)
    axx[2].plot(f_, aZ_ / np.sqrt(tau_))
    x_labels.append(rf"$a={a_:g}$, $\tau={tau_:g}$")
    x_ticks.append(c_ * 5)

axx[2].set_xticks(x_ticks)
axx[2].set_xticklabels(x_labels)
plt.show()

請注意，根據表示形式的不同，峰值的高度會有所不同。只有幅度頻譜的解釋很直接：第二個圖中 \(f_x\) 處的峰值代表正弦訊號幅度 \(|a|\) 的一半。對於所有其他表示形式，都需要考慮持續時間 \(\tau\)，才能提取有關訊號振幅的資訊。

取樣正弦訊號

實際上，取樣訊號被廣泛使用。也就是說，訊號由 \(n\) 個樣本 \(x_k := x(kT)\)、 \(k=0, \ldots, n-1\) 表示，其中 \(T\) 是取樣間隔，\(\tau:=nT\) 是訊號的持續時間，而 \(f_S := 1/T\) 是取樣頻率。請注意，連續訊號需要頻寬限制為 \([-f_S/2, f_S/2]\) 以避免混疊，其中 \(f_S/2\) 稱為奈奎斯特頻率。[5] 將積分替換為總和以計算訊號的能量和功率，即，

\[E_x = T\sum_{k=0}^{n-1} \abs{x_k}^2 = \frac{1}{2}|a|^2\tau\ , \qquad P_x = \frac{1}{\tau}E_x = \frac{1}{2}|a|^2\ ,\]

產生與方程式 (2) 的連續時間情況相同的結果。離散傅立葉轉換 (DFT) 及其反轉換（使用 scipy.fft 模組中有效率的 FFT 計算來實現）由下式給出

\[X_l := \sum_{k=0}^{n-1} x_k \e^{-2\jj\pi k l / n}\ ,\qquad x_k = \frac{1}{n} \sum_{l=0}^{n-1} X_l \e^{2\jj\pi k l / n}\ .\]

DFT 可以解釋為方程式 (3) 的連續傅立葉轉換的未縮放取樣版本，即，

\[X(l\Delta f) = T X_l\ , \quad \Delta f := 1/\tau = 1/(nT)\ .\]

下圖顯示了兩個正弦訊號的幅度頻譜，它們具有單位振幅和頻率 20 Hz 和 20.5 Hz。訊號由 \(n=100\) 個樣本組成，取樣間隔為 \(T=10\) ms，從而產生持續時間 \(\tau=1\) s 和取樣頻率 \(f_S=100\) Hz。

import matplotlib.pyplot as plt
import numpy as np

from scipy.fft import rfft, rfftfreq

n, T = 100, 0.01  # number of samples and sampling interval
fcc = (20, 20.5)  # frequencies of sines
t = np.arange(n) * T
xx = (np.sin(2 * np.pi * fx_ * t) for fx_ in fcc)  # sine signals

f = rfftfreq(n, T)  # frequency bins range from 0 Hz to Nyquist freq.
XX = (rfft(x_) / n for x_ in xx)  # one-sided magnitude spectrum

fg1, ax1 = plt.subplots(1, 1, tight_layout=True, figsize=(6., 3.))
ax1.set(title=r"Magnitude Spectrum (no window) of $x(t) = \sin(2\pi f_x t)$ ",
        xlabel=rf"Frequency $f$ in Hertz (bin width $\Delta f = {f[1]}\,$Hz)",
        ylabel=r"Magnitude $|X(f)|/\tau$", xlim=(f[0], f[-1]))
for X_, fc_, m_ in zip(XX, fcc, ('x-', '.-')):
    ax1.plot(f, abs(X_), m_, label=rf"$f_x={fc_}\,$Hz")

ax1.grid(True)
ax1.legend()
plt.show()

20 Hz 訊號的解釋似乎很直接：除了 20 Hz 之外，所有值均為零。在 20 Hz 時，它是 0.5，這對應於輸入訊號振幅的一半，符合方程式 (1)。另一方面，20.5 Hz 訊號的峰值沿著頻率軸分散。方程式 (3) 顯示，這種差異是由於 20 Hz 是 1 Hz 頻格寬度的倍數，而 20.5 Hz 不是。下圖通過將連續頻譜覆蓋在取樣頻譜之上來說明了這一點

import matplotlib.pyplot as plt
import numpy as np

from scipy.fft import rfft, rfftfreq

n, T = 100, 0.01  # number of samples and sampling interval
tau = n*T
t = np.arange(n) * T
fcc = (20, 20.5)  # frequencies of sines
xx = (np.sin(2 * np.pi * fc_ * t) for fc_ in fcc)  # sine signals

f = rfftfreq(n, T)  # frequency bins range from 0 Hz to Nyquist freq.
XX = (rfft(x_) / n for x_ in xx)  # one-sided FFT normalized to magnitude

i0, i1 = 15, 25
f_cont = np.linspace(f[i0], f[i1], 501)

fg1, axx = plt.subplots(1, 2, sharey='all', tight_layout=True,
                        figsize=(6., 3.))
for c_, (ax_, X_, fx_) in enumerate(zip(axx, XX, fcc)):
    Xc_ = (np.sinc(tau * (f_cont - fx_)) +
           np.sinc(tau * (f_cont + fx_))) / 2
    ax_.plot(f_cont, abs(Xc_), f'-C{c_}', alpha=.5, label=rf"$f_x={fx_}\,$Hz")
    m_line, _, _, = ax_.stem(f[i0:i1+1], abs(X_[i0:i1+1]), markerfmt=f'dC{c_}',
                             linefmt=f'-C{c_}', basefmt=' ')
    plt.setp(m_line, markersize=5)

    ax_.legend(loc='upper left', frameon=False)
    ax_.set(xlabel="Frequency $f$ in Hertz", xlim=(f[i0], f[i1]),
            ylim=(0, 0.59))

axx[0].set(ylabel=r'Magnitude $|X(f)/\tau|$')
fg1.suptitle("Continuous and Sampled Magnitude Spectrum ", x=0.55, y=0.93)
fg1.tight_layout()
plt.show()

頻率的輕微變化會產生外觀顯著不同的幅度頻譜，這對於許多實際應用來說顯然是不希望有的行為。可以使用以下兩種常見技術來改善頻譜

所謂的「零填充」通過在訊號末尾附加零來減小 \(\Delta f\)。為了將頻率過取樣 q 倍，請將參數 n=q*n_x 傳遞給 fft / rfft 函數，其中 n_x 是輸入訊號的長度。

第二種技術稱為視窗化，即將輸入訊號與合適的函數相乘，這樣通常可以抑制次要瓣，但會以加寬主瓣為代價。視窗化的 DFT 可以表示為

(6)\[X^w_l := \sum_{k=0}^{n-1} x_k w_k\e^{-2\jj\pi k l / n}\ ,\]

其中 \(w_k\)、 \(k=0,\ldots,n-1\) 是取樣視窗函數。為了計算先前小節中給出的頻譜表示形式的取樣版本，需要使用以下標準化常數

\[c^\text{amp}:= \abs{\sum_{k=0}^{n-1} w_k}\ ,\qquad c^\text{den} := \sqrt{\sum_{k=0}^{n-1} \abs{w_k}^2}\]

需要使用。第一個常數確保頻譜中的峰值與該頻率下訊號的振幅一致。例如，幅度頻譜可以用 \(|X^w_l / c^\text{amp}|\) 表示。第二個常數保證方程式 (5) 中定義的頻率間隔的功率是一致的。由於不禁止複數值視窗，因此需要絕對值。

下圖顯示了將 hann 視窗和三倍過取樣應用於 \(x(t)\) 的結果

import matplotlib.pyplot as plt
import numpy as np

from scipy.fft import rfft, rfftfreq
from scipy.signal.windows import hann

n, T = 100, 0.01  # number of samples and sampling interval
tau = n*T
q = 3  # over-sampling factor
t = np.arange(n) * T
fcc = (20, 20.5)  # frequencies of sines
xx = [np.sin(2 * np.pi * fc_ * t) for fc_ in fcc]  # sine signals
w = hann(n)
c_w = abs(sum(w))  # normalize constant for window

f_X = rfftfreq(n, T)  # frequency bins range from 0 Hz to Nyquist freq.
XX = (rfft(x_ * w) / c_w for x_ in xx)  # one-sided amplitude spectrum
# Oversampled spectrum:
f_Y = rfftfreq(n*q, T)  # frequency bins range from 0 Hz to Nyquist freq.
YY = (rfft(x_ * w, n=q*n) / c_w for x_ in xx)  # one-sided magnitude spectrum

i0, i1 = 15, 25
j0, j1 = i0*q, i1*q

fg1, axx = plt.subplots(1, 2, sharey='all', tight_layout=True,
                        figsize=(6., 3.))
for c_, (ax_, X_, Y_, fx_) in enumerate(zip(axx, XX, YY, fcc)):
    ax_.plot(f_Y[j0:j1 + 1], abs(Y_[j0:j1 + 1]), f'.-C{c_}',
             label=rf"$f_x={fx_}\,$Hz")
    m_ln, s_ln, _, = ax_.stem(f_X[i0:i1 + 1], abs(X_[i0:i1 + 1]), basefmt=' ',
                              markerfmt=f'dC{c_}', linefmt=f'-C{c_}')
    plt.setp(m_ln, markersize=5)
    plt.setp(s_ln, alpha=0.5)

    ax_.legend(loc='upper left', frameon=False)
    ax_.set(xlabel="Frequency $f$ in Hertz", xlim=(f_X[15], f_X[25]),
            ylim=(0, 0.59))

axx[0].set(ylabel=r'Magnitude $|X(f)/\tau|$')
fg1.suptitle(r"Magnitude Spectrum (Hann window, $%d\times$oversampled)" % q,
             x=0.55, y=0.93)
plt.show()

現在，兩個瓣看起來幾乎相同，並且旁瓣得到了很好的抑制。20.5 Hz 頻譜的最大值也非常接近預期高度的一半。

頻譜能量和頻譜功率可以類似於方程式 (4) 進行計算，從而產生相同的結果，即，

\[E^w_X = T\sum_{k=0}^{n-1} \abs{\frac{X^w_k}{c^\text{den}}}^2 = E_x\ ,\qquad P^w_X = \frac{1}{\tau} E^w_X = \frac{1}{n} \sum_{k=0}^{n-1} \abs{\frac{X^w_k}{c^\text{den}}}^2 = P_x\ .\]

此公式不應與矩形視窗（或無視窗）的特殊情況混淆，即，\(w_k = 1\)、 \(X^w_l=X_l\)、 \(c^\text{den}=\sqrt{n}\)，這會導致

\[E_X = \frac{T}{n}\sum_{k=0}^{n-1} \abs{X_k}^2\ ,\qquad P_X = \frac{1}{n^2} \sum_{k=0}^{n-1} \abs{X_k}^2\ .\]

視窗化頻率離散功率頻譜密度

\[S^w_{xx} := \frac{1}{f_S}\abs{\frac{X^w_l}{c^\text{den}}}^2 = T \abs{\frac{X^w_l}{c^\text{den}}}^2\]

定義於頻率範圍 \([0, f_S)\) 內，並且可以解釋為每頻率間隔 \(\Delta f\) 的功率。在頻帶 \([l_a\Delta f, l_b\Delta f)\) 上積分，如同方程式 (5) 中所示，變成了總和

\[P_X^{a,b} = \Delta f\sum_{k=l_a}^{l_b-1} S^w_{xx} = \frac{1}{nT}\sum_{k=l_a}^{l_b-1} S^w_{xx}\ .\]

視窗化頻率離散能量頻譜密度 \(\tau S^w_{xx}\) 也可以類似地定義。

以上討論顯示，可以定義連續時間情況下頻譜表示的取樣版本。下表總結了這些

視窗化 DFT \(X^w_l\) 的頻譜表示比較，方程式參見 (6)，適用於單位為 \(\text{V}\) 的取樣訊號：

	頻譜	振幅頻譜	能量頻譜密度	功率頻譜密度 (PSD)	振幅頻譜密度
定義	\(\tau X^w_l / c^\text{amp}\)	\(X^w_l / c^\text{amp}\)	\(\tau T |X^w_l / c^\text{den}|^2\)	\(T |X^w_l / c^\text{den}|^2\)	\(\sqrt{T} X^w_l / c^\text{den}\)
單位	\(\text{V} / \text{Hz}\)	\(\text{V}\)	\(\text{V}^2\text{s} / \text{Hz}\)	\(\text{V}^2 / \text{Hz}\)	\(\text{V} / \sqrt{\text{Hz}}\)

請注意，對於密度，由於從積分變為求和以確定頻譜能量/功率，因此在 \(\pm f_x\) 的量值與連續時間情況不同。

雖然漢寧窗是在頻譜分析中最常用的窗函數，但還存在其他窗函數。以下繪圖顯示了 windows 子模組的各種窗函數的量值頻譜。它可以解釋為單一頻率輸入訊號的瓣形狀。請注意，僅顯示了右半部分，並且 \(y\) 軸以分貝為單位，即它是對數縮放的。

import matplotlib.pyplot as plt
import numpy as np

from scipy.fft import rfft, rfftfreq
from scipy.signal import get_window

n, n_zp = 128, 16384  # number of samples without and with zero-padding
t = np.arange(n)
f = rfftfreq(n_zp, 1 / n)

ww = ['boxcar', 'hann', 'hamming', 'tukey', 'blackman', 'flattop']
fg0, axx = plt.subplots(len(ww), 1, sharex='all', sharey='all', figsize=(6., 4.))
for c_, (w_name_, ax_) in enumerate(zip(ww, axx)):
    w_ = get_window(w_name_, n, fftbins=False)
    W_ = rfft(w_ / abs(sum(w_)), n=n_zp)
    W_dB = 20*np.log10(np.maximum(abs(W_), 1e-250))
    ax_.plot(f, W_dB, f'C{c_}-', label=w_name_)
    ax_.text(0.1, -50, w_name_, color=f'C{c_}', verticalalignment='bottom',
             horizontalalignment='left', bbox={'color': 'white', 'pad': 0})
    ax_.set_yticks([-20, -60])
    ax_.grid(axis='x')

axx[0].set_title("Spectral Leakage of various Windows")
fg0.supylabel(r"Normalized Magnitude $20\,\log_{10}|W(f)/c^\operatorname{amp}|$ in dB",
              x=0.04, y=0.5, fontsize='medium')
axx[-1].set(xlabel=r"Normalized frequency $f/\Delta f$ in bins",
            xlim=(0, 9), ylim=(-75, 3))

fg0.tight_layout(h_pad=0.4)
plt.show()

此圖顯示，窗函數的選擇通常是在主瓣的寬度和旁瓣的高度之間進行權衡。請注意，boxcar 窗對應於 \(\rect\) 函數，即不進行視窗化。此外，許多描繪的窗函數更常在濾波器設計中使用，而不是在頻譜分析中使用。

頻譜的相位

傅立葉轉換的相位（即 angle）通常用於研究訊號通過系統（如濾波器）的頻譜成分的時間延遲。在以下範例中，標準測試訊號（單位功率的脈衝）通過一個簡單的濾波器，該濾波器將輸入延遲三個樣本。輸入由 \(n=50\) 個樣本組成，取樣間隔為 \(T = 1\) 秒。該圖顯示了輸入和輸出訊號在頻率上的量值和相位

import matplotlib.gridspec as gridspec
import matplotlib.pyplot as plt
import numpy as np

from scipy import signal
from scipy.fft import rfft, rfftfreq

# Create input signal:
n = 50
x = np.zeros(n)
x[0] = n

# Apply FIR filter which delays signal by 3 samples:
y = signal.lfilter([0, 0, 0, 1], 1, x)

X, Y = (rfft(z_) / n for z_ in (x, y))
f = rfftfreq(n, 1)  # sampling interval T = 1 s

fig = plt.figure(tight_layout=True, figsize=(6., 4.))
gs = gridspec.GridSpec(3, 1)
ax0 = fig.add_subplot(gs[0, :])
ax1 = fig.add_subplot(gs[1:, :], sharex=ax0)

for Z_, n_, m_ in zip((X, Y), ("Input $X(f)$", "Output $Y(f)$"), ('+-', 'x-')):
    ax0.plot(f, abs(Z_), m_, alpha=.5, label=n_)
    ax1.plot(f, np.unwrap(np.angle(Z_)), m_, alpha=.5, label=n_)

ax0.set(title="Frequency Response of 3 Sample Delay Filter (no window)",
        ylabel="Magnitude", xlim=(0, f[-1]), ylim=(0, 1.1))
ax1.set(xlabel=rf"Frequency $f$ in Hertz ($\Delta f = {1/n}\,$Hz)",
        ylabel=r"Phase in rad")
ax1.set_yticks(-np.arange(0, 7)*np.pi/2,
               ['0', '-π/2', '-π', '-3/2π', '-2π', '-4/3π', '-3π'])

ax2 = ax1.twinx()
ax2.set(ylabel=r"Phase in Degrees", ylim=np.rad2deg(ax1.get_ylim()),
        yticks=np.arange(-540, 90, 90))
for ax_ in (ax0, ax1):
    ax_.legend()
    ax_.grid()
plt.show()

輸入具有單位量值和零相位傅立葉轉換，這是用作測試訊號的原因。輸出也具有單位量值，但相位線性下降，斜率為 \(-6\pi\)。這是預期的，因為將訊號 \(x(t)\) 延遲 \(\Delta t\) 會在其傅立葉轉換中產生額外的線性相位項，即，

\[\int_\IR x(t-\Delta t) \e^{-\jj2\pi f t}\dd t = X(f)\, e^{-\jj2\pi \Delta t f}\ .\]

請注意，在圖中，相位不限於區間 \((+\pi, \pi]\)（angle 的輸出），因此沒有任何不連續性。這是通過使用 numpy.unwrap 函數實現的。如果濾波器的傳遞函數已知，則可以使用 freqz 直接確定濾波器的頻譜響應。

具有平均的頻譜

periodogram 函數計算功率頻譜密度（scaling='density'）或平方量值頻譜（scaling='spectrum'）。為了獲得平滑的週期圖，可以使用 welch 函數。它通過將輸入訊號劃分為重疊的段來進行平滑處理，然後計算每個段的視窗化 DFT。結果是這些 DFT 的平均值。

下面的範例顯示了由 \(1.27\,\text{kHz}\) 正弦訊號（振幅 \(\sqrt{2}\,\text{V}\)）和加性高斯雜訊（頻譜功率密度平均值為 \(10^{-3}\,\text{V}^2/\text{Hz}\)）組成的訊號的平方量值頻譜和功率頻譜密度。

import matplotlib.pyplot as plt
import numpy as np
import scipy.signal as signal

rng = np.random.default_rng(73625)  # seeding for reproducibility

fs, n = 10e3, 10_000
f_x, noise_power = 1270, 1e-3 * fs / 2
t = np.arange(n) / fs
x = (np.sqrt(2) * np.sin(2 * np.pi * f_x * t) +
     rng.normal(scale=np.sqrt(noise_power), size=t.shape))

fg, axx = plt.subplots(1, 2, sharex='all', tight_layout=True, figsize=(7, 3.5))
axx[0].set(title="Squared Magnitude Spectrum", ylabel="Square of Magnitude in V²")
axx[1].set(title="Power Spectral Density", ylabel="Power Spectral Density in V²/Hz")
for ax_, s_ in zip(axx, ('spectrum', 'density')):
    f_p, P_p = signal.periodogram(x, fs, 'hann', scaling=s_)
    f_w, P_w = signal.welch(x, fs, scaling=s_)
    ax_.semilogy(f_p/1e3, P_p, label=f"Periodogram ({len(f_p)} bins)")
    ax_.semilogy(f_w/1e3, P_w, label=f"Welch's Method ({len(f_w)} bins)")
    ax_.set(xlabel="Frequency in kHz", xlim=(0, 2), ylim=(1e-7, 1.3))
    ax_.grid(True)
    ax_.legend(loc='lower center')
plt.show()

這些圖顯示，welch 函數產生更平滑的雜訊基底，但犧牲了頻率解析度。由於平滑，正弦波瓣的高度更寬，不如週期圖中的高。左圖可用於讀取波瓣的高度，即正弦波平方量值的一半 \(1\,\text{V}^2\)。右圖可用於確定雜訊基底 \(10^{-3}\,\text{V}^2/\text{Hz}\)。請注意，由於頻率解析度有限，平均平方量值頻譜的波瓣高度不完全為 1。可以利用零填充（例如，將 nfft=4*len(x) 傳遞給 welch）或通過增加段長（設定參數 nperseg）來減少段數，以增加頻率箱的數量。

隆伯-史卡格週期圖 (lombscargle)

最小平方頻譜分析 (LSSA) [6] [7] 是一種估計頻譜的方法，基於正弦波對數據樣本的最小平方擬合，類似於傅立葉分析。傅立葉分析是科學中最常用的頻譜方法，通常會增強長間隙記錄中的長週期雜訊；LSSA 減輕了這些問題。

隆伯-史卡格方法對不均勻取樣的數據執行頻譜分析，並且已知是尋找和測試微弱週期訊號顯著性的有效方法。

對於包含 \(N_{t}\) 個測量值 \(X_{j}\equiv X(t_{j})\) 的時間序列，在時間 \(t_{j}\) 取樣，其中 \((j = 1, \ldots, N_{t})\)，假設已縮放和平移，使其平均值為零，方差為單位，則頻率 \(f\) 處的歸一化隆伯-史卡格週期圖為

\[P_{n}(f) = \frac{1}{2}\left\{\frac{\left[\sum_{j}^{N_{t}}X_{j}\cos\omega(t_{j}-\tau)\right]^{2}}{\sum_{j}^{N_{t}}\cos^{2}\omega(t_{j}-\tau)}+\frac{\left[\sum_{j}^{N_{t}}X_{j}\sin\omega(t_{j}-\tau)\right]^{2}}{\sum_{j}^{N_{t}}\sin^{2}\omega(t_{j}-\tau)}\right\}.\]

這裡，\(\omega \equiv 2\pi f\) 是角頻率。頻率相關的時間偏移 \(\tau\) 由下式給出

\[\tan 2\omega\tau = \frac{\sum_{j}^{N_{t}}\sin 2\omega t_{j}}{\sum_{j}^{N_{t}}\cos 2\omega t_{j}}.\]

lombscargle 函數使用 Zechmeister 和 Kürster [8] 創建的稍微修改過的演算法來計算週期圖，該演算法允許對個別樣本進行加權，並獨立地為每個頻率計算未知的偏移量（也稱為“浮動平均值”）。

短時傅立葉轉換

本節提供有關使用 ShortTimeFFT 類的一些背景資訊：短時傅立葉轉換 (STFT) 可用於分析訊號隨時間變化的頻譜特性。它通過使用滑動窗口將訊號劃分為重疊的塊，並計算每個塊的傅立葉轉換。對於連續時間複數值訊號 \(x(t)\)，STFT 定義為 [9]

\[S(f, t) := \int_\IR x(\xi)\, \conj{w(\xi-t)}\,\e^{-\jj2\pi f \xi}\dd\xi\ ,\]

其中 \(w(t)\) 是複數值窗函數，其共軛複數為 \(\conj{w(t)}\)。它可以解釋為確定 \(x\) 與窗口 \(w\) 的純量積，該窗口通過時間 \(t\) 平移，然後通過頻率 \(f\) 調變（即頻率偏移）。為了處理取樣訊號 \(x[k] := x(kT)\)，\(k\in\IZ\)，取樣間隔為 \(T\)（為取樣頻率 fs 的倒數），需要使用離散版本，即僅在離散網格點 \(S[q, p] := S(q \Delta f, p\Delta t)\)，\(q,p\in\IZ\) 評估 STFT。它可以表示為

(7)\[S[q,p] = \sum_{k=0}^{N-1} x[k]\,\conj{w[k-p h]}\, \e^{-\jj2\pi q k / N}\ , \quad q,p\in\IZ\ ,\]

其中 p 表示 \(S\) 的時間索引，時間間隔為 \(\Delta t := h T\)，\(h\in\IN\)（參見 delta_t），它可以表示為 \(h\) 個樣本的 hop 大小。\(q\) 表示 \(S\) 的頻率索引，步長為 \(\Delta f := 1 / (N T)\)（參見 delta_f），這使其與 FFT 相容。\(w[m] := w(mT)\)，\(m\in\IZ\) 是取樣窗函數。

為了更符合 ShortTimeFFT 的實作，將方程式 (7) 重新表述為兩個步驟的過程是有意義的

1. 通過使用由 \(M\) 個樣本組成的窗口 \(w[m]\)（參見 m_num）在 \(t[p] := p \Delta t = h T\)（參見 delta_t）為中心的視窗化來提取第 \(p\) 個切片，即，

   (8)\[x_p[m] = x\!\big[m - \lfloor M/2\rfloor + h p\big]\, \conj{w[m]}\ , \quad m = 0, \ldots M-1\ ,\]

   其中整數 \(\lfloor M/2\rfloor\) 表示 M//2，即它是窗口的中點（m_num_mid）。為了符號方便，假設對於 \(k\not\in\{0, 1, \ldots, N-1\}\)，\(x[k]:=0\)。在小節 滑動窗口 中，更詳細地討論了切片的索引。

2. 然後執行 \(x_p[m]\) 的離散傅立葉轉換（即 FFT）。

   (9)\[S[q, p] = \sum_{m=0}^{M-1} x_p[m] \exp\!\big\{% -2\jj\pi (q + \phi_m)\, m / M\big\}\ .\]

   請注意，可以指定線性相位 \(\phi_m\)（參見 phase_shift），這對應於將輸入移動 \(\phi_m\) 個樣本。預設值為 \(\phi_m = \lfloor M/2\rfloor\)（根據定義，對應於 phase_shift = 0），這抑制了未移動訊號的線性相位分量。此外，FFT 可以通過用零填充 \(w[m]\) 來過度取樣。這可以通過指定 mfft 大於窗口長度 m_num 來實現——這將 \(M\) 設置為 mfft（意味著 \(m\not\in\{0, 1, \ldots, M-1\}\) 時 \(w[m]:=0\) 也成立）。

反短時傅立葉轉換（istft）是通過反轉這兩個步驟來實現的

1. 執行反離散傅立葉轉換，即，

   \[x_p[m] = \frac{1}{M}\sum_{q=0}^M S[q, p]\, \exp\!\big\{ 2\jj\pi (q + \phi_m)\, m / M\big\}\ .\]

2. 將由 \(w_d[m]\) 加權的移動切片相加，以重建原始訊號，即，

   \[x[k] = \sum_p x_p\!\big[\mu_p(k)\big]\, w_d\!\big[\mu_p(k)\big]\ ,\quad \mu_p(k) = k + \lfloor M/2\rfloor - h p\]

   對於 \(k \in [0, \ldots, n-1]\)。\(w_d[m]\) 是 \(w[m]\) 的所謂規範對偶窗口，並且也由 \(M\) 個樣本組成。

請注意，並非所有窗口和跳躍大小都必然存在反 STFT。對於給定的窗口 \(w[m]\)，跳躍大小 \(h\) 必須足夠小，以確保 \(x[k]\) 的每個樣本都至少被一個窗口切片的非零值觸及。這有時被稱為“非零重疊條件”（參見 check_NOLA）。在小節 反 STFT 和對偶窗口 中給出了一些更詳細的資訊。

滑動窗口

本小節通過範例討論了 ShortTimeFFT 中滑動窗口的索引方式：考慮一個長度為 6 的窗口，hop 間隔為 2，取樣間隔 T 為 1，例如，ShortTimeFFT (np.ones(6), 2, fs=1)。下圖示意性地描述了前四個窗口位置，也稱為時間切片

x 軸表示時間 \(t\)，它對應於樣本索引 k，由藍色方框的底行指示。y 軸表示時間切片索引 \(p\)。訊號 \(x[k]\) 從索引 \(k=0\) 開始，並以淺藍色背景標記。根據定義，第零個切片（\(p=0\)）以 \(t=0\) 為中心。每個切片的中心（m_num_mid），此處為樣本 6//2=3，以文字“mid”標記。預設情況下，stft 計算所有與訊號有部分重疊的切片。因此，第一個切片在 p_min = -1，最低樣本索引為 k_min = -5。不受切片向訊號左側伸出影響的第一個樣本索引是 \(p_{lb} = 2\)，不受邊界效應影響的第一個樣本索引是 \(k_{lb} = 5\)。lower_border_end 屬性返回元組 \((k_{lb}, p_{lb})\)。

訊號末端的行為在下方針對具有 \(n=50\) 個樣本的訊號進行了描述，如藍色背景所示

這裡最後一個切片的索引為 \(p=26\)。因此，按照 Python 關於範圍外的結束索引的約定，p_max = 27 表示第一個不接觸訊號的切片。對應的樣本索引為 k_max = 55。第一個向右伸出的切片是 \(p_{ub} = 24\)，其第一個樣本位於 \(k_{ub}=45\)。upper_border_begin 函數返回元組 \((k_{ub}, p_{ub})\)。

反 STFT 和對偶窗口

術語對偶窗口源於框架理論 [10]，其中框架是一種級數展開，可以表示給定希爾伯特空間中的任何函數。在那裡，如果對於給定希爾伯特空間 \(\mathcal{H}\) 中的所有函數 \(f\)，展開 \(\{g_k\}\) 和 \(\{h_k\}\) 是對偶框架

\[f = \sum_{k\in\IN} \langle f, g_k\rangle h_k = \sum_{k\in\IN} \langle f, h_k\rangle g_k\ , \quad f \in \mathcal{H}\ ,\]

成立，其中 \(\langle ., .\rangle\) 表示 \(\mathcal{H}\) 的純量積。所有框架都有對偶框架 [10]。

僅在離散網格點 \(S(q \Delta f, p\Delta t)\) 評估的 STFT 在文獻中稱為“Gabor 框架”[9] [10]。由於窗口 \(w[m]\) 的支持限制在有限區間，ShortTimeFFT 屬於所謂的“無痛非正交展開”類別 [9]。在這種情況下，對偶窗口始終具有相同的支持，並且可以通過反轉對角矩陣來計算。以下將概述一個粗略的推導，僅需要一些對矩陣操作的理解

由於方程式 (7) 中給出的 STFT 是 \(x[k]\) 中的線性映射，因此可以用向量-矩陣符號表示。這使我們能夠通過線性最小平方方法的正式解（如 lstsq 中）來表示逆，這將導致一個優美而簡單的結果。

我們首先重新表述方程式 (8) 的視窗化

(10)\[\begin{split} \vb{x}_p = \vb{W}_{\!p}\,\vb{x} = \begin{bmatrix} \cdots & 0 & w[0] & 0 & \cdots&&&\\ & \cdots & 0 & w[1] & 0 & \cdots&&\\ & & & & \ddots&&&\\ &&\cdots & 0 & 0 & w[M-1] & 0 & \cdots \end{bmatrix}\begin{bmatrix} x[0]\\ x[1]\\ \vdots\\ x[N-1] \end{bmatrix}\ ,\end{split}\]

其中 \(M\times N\) 矩陣 \(\vb{W}_{\!p}\) 僅在第 \((ph)\) 個次對角線上具有非零條目，即，

(11)\[\begin{split} W_p[m,k] = w[m]\, \delta_{m+ph,k}\ ,\quad \delta_{k,l} &= \begin{cases} 1 & \text{ for } k=l\ ,\\ 0 & \text{ elsewhere ,} \end{cases}\end{split}\]

其中 \(\delta_{k,l}\) 是克羅內克 Delta。方程式 (9) 可以表示為

\[\vb{s}_p = \vb{F}\,\vb{x}_p \quad\text{with}\quad F[q,m] =\exp\!\big\{-2\jj\pi (q + \phi_m)\, m / M\big\}\ ,\]

這使得第 \(p\) 個切片的 STFT 可以寫為

(12)\[\vb{s}_p = \vb{F}\vb{W}_{\!p}\,\vb{x} =: \vb{G}_p\,\vb{x} \quad\text{with}\quad s_p[q] = S[p,q]\ .\]

請注意，\(\vb{F}\) 是么正的，即，逆等於其共軛轉置，意味著 \(\conjT{\vb{F}}\vb{F} = \vb{I}\)。

為了獲得 STFT 的單個向量-矩陣方程式，將切片堆疊成一個向量，即，

\[\begin{split}\vb{s} := \begin{bmatrix} \vb{s}_0\\ \vb{s}_1\\ \vdots\\ \vb{s}_{P-1} \end{bmatrix} = \begin{bmatrix} \vb{G}_0\\ \vb{G}_1\\ \vdots\\ \vb{G}_{P-1} \end{bmatrix}\, \vb{x} =: \vb{G}\, \vb{x}\ ,\end{split}\]

其中 \(P\) 是結果 STFT 的列數。為了反轉此方程式，可以使用 Moore-Penrose 逆 \(\vb{G}^\dagger\)

(13)\[\vb{x} = \inv{\conjT{\vb{G}}\vb{G}}\, \conjT{\vb{G}} \vb{s} =: \vb{G}^\dagger \vb{s}\ ,\]

如果

\[\begin{split}\vb{D} := \conjT{\vb{G}}\vb{G} = \begin{bmatrix} \conjT{\vb{G}_0}& \conjT{\vb{G}_1}& \cdots & \conjT{\vb{G}_{P-1}} \end{bmatrix}^T \begin{bmatrix} \vb{G}_0\\ \vb{G}_1\\ \vdots\\ \vb{G}_{P-1} \end{bmatrix} = \sum_{p=0}^{P-1} \conjT{\vb{G}_p}\vb{G}_p \ .\end{split}\]

是可逆的，則存在。然後 \(\vb{x} = \vb{G}^\dagger\vb{G}\,\vb{x} = \inv{\conjT{\vb{G}}\vb{G}}\,\conjT{\vb{G}}\vb{G}\,\vb{x}\) 顯然成立。\(\vb{D}\) 始終是對角矩陣，具有非負對角線條目。當將 \(\vb{D}\) 進一步簡化為

(14)\[\vb{D} = \sum_{p=0}^{P-1} \conjT{\vb{G}_p}\vb{G}_p = \sum_{p=0}^{P-1} \conjT{(\vb{F}\,\vb{W}_{\!p})}\, \vb{F}\,\vb{W}_{\!p} = \sum_{p=0}^{P-1} \conjT{\vb{W}_{\!p}}\vb{W}_{\!p} =: \sum_{p=0}^{P-1} \vb{D}_p\]

由於 \(\vb{F}\) 是么正的。此外

(15)\[\begin{split}D_p[r,s] &= \sum_{m=0}^{M-1} \conj{W_p^T[r,m]}\,W_p[m,s] = \sum_{m=0}^{M-1} \left(\conj{w[m]}\, \delta_{m+ph,r}\right) \left(w[m]\, \delta_{m+ph,s}\right)\\ &= \sum_{m=0}^{M-1} \big|w[m]\big|^2\, \delta_{r,s}\, \delta_{r,m+ph}\ .\end{split}\]

表明 \(\vb{D}_p\) 是具有非負條目的對角矩陣。因此，對 \(\vb{D}_p\) 求和保留了該屬性。這允許進一步簡化方程式 (13)，即，

(16)\[\begin{split}\vb{x} &= \vb{D}^{-1} \conjT{\vb{G}}\vb{s} = \sum_{p=0}^{P-1} \vb{D}^{-1}\conjT{\vb{W}_{\!p}}\, \conjT{\vb{F}}\vb{s}_p = \sum_{p=0}^{P-1} (\conj{\vb{W}_{\!p}\vb{D}^{-1}})^T\, \conjT{\vb{F}}\vb{s}_p\\ &=: \sum_{p=0}^{P-1}\conjT{\vb{U}_p}\,\conjT{\vb{F}}\vb{s}_p\ .\end{split}\]

利用方程式 (11)、(14)、(15)，\(\vb{U}_p=\vb{W}_{\!p}\vb{D}^{-1}\) 可以表示為

(17)\[\begin{split}U_p[m, k] &= W[m,k]\, D^{-1}[k,k] = \left(w[m] \delta_{m+ph,k}\right) \inv{\sum_{\eta=0}^{P-1} \vb{D}_\eta[k,k]} \delta_{m+ph,k}\\ &= w[m] \inv{\sum_{\eta=0}^{P-1}\sum_{\mu=0}^{M-1} \big|w[\mu]\big|^2\,\delta_{m+ph, \mu+\eta h}} \delta_{m+ph,k}\\ &= w[m] \inv{\sum_{\eta=0}^{P-1} \big|w[m+(p-\eta)h]\big|^2} \delta_{m+ph,k} \ .\end{split}\]

這顯示 \(\vb{U}_p\) 具有與 Eq. (11) 中的 \(\vb{W}_p\) 相同的結構，亦即，僅在 \((ph)\)-th 次對角線上具有非零項。反矩陣中的總和項可以解釋為將 \(|w[\mu]|^2\) 滑動於 \(w[m]\) 之上（並包含反轉），因此只有與 \(w[m]\) 重疊的成分才會產生影響。因此，所有遠離邊界的 \(U_p[m, k]\) 都是相同的窗函數。為了規避邊界效應，\(x[k]\) 會以零填充，擴大 \(\vb{U}\)，使得所有接觸 \(x[k]\) 的切片都包含相同的對偶窗函數。

\[w_d[m] = w[m] \inv{\sum_{\eta\in\IZ} \big|w[m + \eta\, h]\big|^2}\ .\]

由於 \(w[m] = 0\) 對於 \(m \not\in\{0, \ldots, M-1\}\) 成立，因此僅需要對滿足 \(|\eta| < M/h\) 的索引 \(\eta\) 進行求和。對偶窗函數的名稱可以透過將 Eq. (12) 插入 Eq. (16) 來證明，即：

\[\vb{x} = \sum_{p=0}^{P-1} \conjT{\vb{U}_p}\,\conjT{\vb{F}}\, \vb{F}\,\vb{W}_{\!p}\,\vb{x} = \left(\sum_{p=0}^{P-1} \conjT{\vb{U}_p}\,\vb{W}_{\!p}\right)\vb{x}\ ,\]

顯示 \(\vb{U}_p\) 和 \(\vb{W}_{\!p}\) 是可互換的。因此，\(w_d[m]\) 也是一個有效的窗函數，其對偶窗函數為 \(w[m]\)。請注意，由於 \(\vb{s}\) 通常比 \(\vb{x}\) 具有更多條目，因此 \(w_d[m]\) 不是唯一的對偶窗函數。可以證明，\(w_d[m]\) 具有最小能量（或 \(L_2\) 範數）[9]，這也是它被命名為「正規對偶窗函數」的原因。

與舊版實作的比較

函數 stft、istft 和 spectrogram 早於 ShortTimeFFT 實作。本節討論舊版「傳統」實作與較新的 ShortTimeFFT 實作之間的主要差異。重新編寫的主要動機是意識到，若不破壞相容性，就無法以合理的方式整合對偶窗函數。這為重新思考程式碼結構和參數化提供了機會，從而使一些隱含的行為更加明確。

以下範例比較了具有負斜率的複數值啁啾訊號的兩種 STFT

>>> import matplotlib.pyplot as plt
>>> import numpy as np
>>> from scipy.fft import fftshift
>>> from scipy.signal import stft, istft, spectrogram, ShortTimeFFT
...
>>> fs, N = 200, 1001  # 200 Hz sampling rate for 5 s signal
>>> t_z = np.arange(N) / fs  # time indexes for signal
>>> z = np.exp(2j*np.pi*70 * (t_z - 0.2*t_z**2))  # complex-valued chirp
...
>>> nperseg, noverlap = 50, 40
>>> win = ('gaussian', 1e-2 * fs)  # Gaussian with 0.01 s standard dev.
...
>>> # Legacy STFT:
>>> f0_u, t0, Sz0_u = stft(z, fs, win, nperseg, noverlap,
...                        return_onesided=False, scaling='spectrum')
>>> f0, Sz0 = fftshift(f0_u), fftshift(Sz0_u, axes=0)
...
>>> # New STFT:
>>> SFT = ShortTimeFFT.from_window(win, fs, nperseg, noverlap, fft_mode='centered',
...                                scale_to='magnitude', phase_shift=None)
>>> Sz1 = SFT.stft(z)
...
>>> # Plot results:
>>> fig1, axx = plt.subplots(2, 1, sharex='all', sharey='all',
...                          figsize=(6., 5.))  # enlarge figure a bit
>>> t_lo, t_hi, f_lo, f_hi = SFT.extent(N, center_bins=True)
>>> axx[0].set_title(r"Legacy stft() produces $%d\times%d$ points" % Sz0.T.shape)
>>> axx[0].set_xlim(t_lo, t_hi)
>>> axx[0].set_ylim(f_lo, f_hi)
>>> axx[1].set_title(r"ShortTimeFFT produces $%d\times%d$ points" % Sz1.T.shape)
>>> axx[1].set_xlabel(rf"Time $t$ in seconds ($\Delta t= {SFT.delta_t:g}\,$s)")
...
>>> # Calculate extent of plot with centered bins since
>>> # imshow does not interpolate by default:
>>> dt2 = (nperseg-noverlap) / fs / 2  # equals SFT.delta_t / 2
>>> df2 = fs / nperseg / 2  # equals SFT.delta_f / 2
>>> extent0 = (-dt2, t0[-1] + dt2, f0[0] - df2, f0[-1] - df2)
>>> extent1 = SFT.extent(N, center_bins=True)
...
>>> kw = dict(origin='lower', aspect='auto', cmap='viridis')
>>> im1a = axx[0].imshow(abs(Sz0), extent=extent0, **kw)
>>> im1b = axx[1].imshow(abs(Sz1), extent=extent1, **kw)
>>> fig1.colorbar(im1b, ax=axx, label="Magnitude $|S_z(t, f)|$")
>>> _ = fig1.supylabel(r"Frequency $f$ in Hertz ($\Delta f = %g\,$Hz)" %
...                    SFT.delta_f, x=0.08, y=0.5, fontsize='medium')
>>> plt.show()

與舊版相比，ShortTimeFFT 產生多 3 個時間切片是主要差異。如滑動窗函數章節中所述，所有接觸訊號的切片都包含在新版本中。這具有優勢，因為 STFT 可以像 ShortTimeFFT 程式碼範例中所示進行切片和重新組裝。此外，使用所有接觸切片使 ISTFT 在窗函數在某處為零的情況下更加穩健。

請注意，具有相同時間戳記的切片會產生相同的結果（在數值精度範圍內），即

>>> np.allclose(Sz0, Sz1[:, 2:-1])
True

一般而言，這些額外的切片包含非零值。由於我們範例中存在大量重疊，因此它們非常小。例如：

>>> abs(Sz1[:, 1]).min(), abs(Sz1[:, 1]).max()
(6.925060911593139e-07, 8.00271269218721e-07)

ISTFT 可用於重建原始訊號

>>> t0_r, z0_r = istft(Sz0_u, fs, win, nperseg, noverlap,
...                    input_onesided=False, scaling='spectrum')
>>> z1_r = SFT.istft(Sz1, k1=N)
...
>>> len(z0_r), len(z)
(1010, 1001)
>>> np.allclose(z0_r[:N], z)
True
>>> np.allclose(z1_r, z)
True

請注意，舊版實作返回的訊號比原始訊號長。另一方面，新的 istft 允許明確指定重建訊號的起始索引 k0 和結束索引 k1。舊版實作中的長度差異是由於訊號長度不是切片的倍數所造成的。

此範例中新舊版本之間的進一步差異包括

• 參數 fft_mode='centered' 確保對於圖中的雙邊 FFT，零頻率在垂直方向上居中。使用舊版實作，需要使用 fftshift。fft_mode='twosided' 產生與舊版本相同的行為。

• 參數 phase_shift=None 確保兩個版本的相位相同。ShortTimeFFT 的預設值 0 產生具有額外線性相位項的 STFT 切片。

頻譜圖定義為 STFT 的絕對平方[9]。spectrogram 由 ShortTimeFFT 提供，堅持該定義，即：

>>> np.allclose(SFT.spectrogram(z), abs(Sz1)**2)
True

另一方面，舊版 spectrogram 提供了另一個 STFT 實作，其主要差異在於訊號邊界的處理方式不同。以下範例展示如何使用 ShortTimeFFT 來獲得與舊版 spectrogram 產生的 SFT 相同的結果

>>> # Legacy spectrogram (detrending for complex signals not useful):
>>> f2_u, t2, Sz2_u = spectrogram(z, fs, win, nperseg, noverlap,
...                               detrend=None, return_onesided=False,
...                               scaling='spectrum', mode='complex')
>>> f2, Sz2 = fftshift(f2_u), fftshift(Sz2_u, axes=0)
...
>>> # New STFT:
... SFT = ShortTimeFFT.from_window(win, fs, nperseg, noverlap,
...                                fft_mode='centered',
...                                scale_to='magnitude', phase_shift=None)
>>> Sz3 = SFT.stft(z, p0=0, p1=(N-noverlap)//SFT.hop, k_offset=nperseg//2)
>>> t3 = SFT.t(N, p0=0, p1=(N-noverlap)//SFT.hop, k_offset=nperseg//2)
...
>>> np.allclose(t2, t3)
True
>>> np.allclose(f2, SFT.f)
True
>>> np.allclose(Sz2, Sz3)
True

與其他 STFT 的不同之處在於，時間切片不是從 0 開始，而是從 nperseg//2 開始，即：

>>> t2
array([0.125, 0.175, 0.225, 0.275, 0.325, 0.375, 0.425, 0.475, 0.525,
       ...
       4.625, 4.675, 4.725, 4.775, 4.825, 4.875])

此外，僅返回不向右突出的切片，將最後一個切片置中於 4.875 秒，這使其比預設的 stft 參數化更短。

使用 mode 參數，舊版 spectrogram 也可以返回 ‘angle’、‘phase’、‘psd’ 或 ‘magnitude’。scaling 舊版 spectrogram 的行為並不直接，因為它取決於參數 mode、scaling 和 return_onesided。ShortTimeFFT 中並非所有組合都有直接的對應關係，因為它僅提供 ‘magnitude’、‘psd’ 或根本沒有窗函數的 scaling。下表顯示了這些對應關係

舊版 spectrogram				ShortTimeFFT
mode	scaling	return_onesided		fft_mode	scaling
psd	density	True		onesided2X	psd
psd	density	False		twosided	psd
magnitude	spectrum	True		onesided	magnitude
magnitude	spectrum	False		twosided	magnitude
complex	spectrum	True		onesided	magnitude
complex	spectrum	False		twosided	magnitude
psd	spectrum	True		—	—
psd	spectrum	False		—	—
complex	density	True		—	—
complex	density	False		—	—
magnitude	density	True		—	—
magnitude	density	False		—	—
—	—	—		*	None

當在複數值輸入訊號上使用 onesided 輸出時，舊版 spectrogram 會切換到 two-sided 模式。ShortTimeFFT 引發 TypeError，因為使用的 rfft 函數僅接受實數值輸入。有關各種參數化引起的各種頻譜表示的討論，請參閱上面的頻譜分析章節。

參考文獻

一些延伸閱讀和相關軟體

[1]

“傅立葉變換”，維基百科，https://en.wikipedia.org/wiki/Fourier_transform

[2]

“頻譜密度”，維基百科，https://en.wikipedia.org/wiki/Spectral_density

[3]

“離散傅立葉變換”，維基百科，https://en.wikipedia.org/wiki/Discrete_Fourier_transform

[4]

“頻譜洩漏”，維基百科，https://en.wikipedia.org/wiki/Spectral_leakage

[5]

“奈奎斯特-香農採樣定理”，維基百科，https://en.wikipedia.org/wiki/Nyquist-Shannon_sampling_theorem

[6]

N.R. Lomb “不均勻間隔資料的最小平方頻率分析”，《天文學和太空科學》，第 39 卷，第 447-462 頁，1976 年

[7]

J.D. Scargle “天文時間序列分析研究。II - 不均勻間隔資料的頻譜分析的統計方面”，《天文物理期刊》，第 263 卷，第 835-853 頁，1982 年

[8]

M. Zechmeister 和 M. Kürster，“廣義 Lomb-Scargle 週期圖。浮動平均值和克卜勒週期圖的新形式主義”，《天文學和天文物理學》，第 496 卷，第 577-584 頁，2009 年

[9] (1,2,3,4,5)

Karlheinz Gröchenig：“時頻分析基礎”，Birkhäuser Boston 2001 年，DOI:10.1007/978-1-4612-0003-1

[10] (1,2,3)

Ole Christensen：“幀和 Riesz 基底導論”，Birkhäuser Boston 2016 年，DOI:10.1007/978-3-319-25613-9