5. 頻譜分析

頻譜分析 (Spectrum Analysis) 目的是分析訊號在不同頻率下的分量，並且以圖形表示。

頻譜分析的方法透過傅立葉轉換形成傅立葉頻譜 (Fourier Spectrum) 。

本章節將透過實際例子著重說明與討論Python的FFT相關工具使用與結果解讀。

5.1 傅立葉頻譜分析弦波

若弦波定義為$x(t)=A\cos\left(\omega_0t+\phi\right)$
或
$x(t)=A\cos\left(2\pi f_0t+\phi\right)$

經由Euler’s Formula
可得到
$x(t)=\frac A2e^{j(\omega_0t+\phi)}+\frac A2e^{-j(\omega_0t+\phi)}$

經由傅立葉轉換
$
\begin{aligned}
\mathcal F\{x(t)\} & = \mathcal F{\frac A2e^{j(\omega_0t+\phi)}+\frac A2e^{-j(\omega_0t+\phi)}}\\
& = \frac A2\cdot e^{j\phi}\cdot\mathcal F{e^{j\omega_0t}}+\frac A2\cdot e^{-j\phi}\cdot\mathcal F{e^{-j\omega_0t}}\\
& = \frac A2\cdot e^{j\phi}\cdot2\mathrm{πδ}(\mathrm\omega-{\mathrm\omega}_0)+\frac A2\cdot e^{-j\phi}\cdot2\mathrm{πδ}(\mathrm\omega+{\mathrm\omega}_0)\\
& = A\pi(e^{j\phi}\delta(\mathrm\omega-{\mathrm\omega}_0)+e^{-j\phi}\delta(\mathrm\omega+{\mathrm\omega}_0))
\end{aligned}
$

振幅頻譜 (Amplitude Spectrum)

經由上述轉換公式可得到弦波包含了兩個振幅 (Amplitude)的分量，當頻率為 $\omega_0$, $-\omega_0$，振幅均為 $A\pi$，如果以圖形表示，就稱為振幅頻譜 (Amplitude Spectrum)。

相位頻譜 (Phase Spectrum)

弦波 $x(t)$ 也包含了兩個相位移 (Phase Shift) 的分相，當頻率為 $\omega_0$, $-\omega_0$，相位移分別為 $\phi$, $-\phi$，如果以圖形表示，稱為相位頻譜 (Phase Spectrum)。

5.2 傅立葉頻譜實際應用

範例 — 簡單序列分析

若數位訊號為：
$x=\{4.0,\;2.0,\;1.0,\;-3.0,\;1.5\},\;n=0,\;1,\;2,\;3,\;4$
假設取樣頻率為1000 Hz
計算其強度頻譜 (Magnitude Spectrum)，並以圖形表示

5.2.1 解法 ⇒ FFT結果直接呈現，不完整的結果

程式碼

  import numpy as np
  from numpy.fft import fft
  import matplotlib.pyplot as plt

  # 訊號點
  x = np.array([4.0, 2.0, 1.0, -3.0, 1.5])
  # 訊號點時間
  x_index = [i for i in range(0, 5, 1)]

  # 求傅立葉轉換
  y = fft(x)
  # 計算頻譜振幅
  y_m = abs(y)
  print(y)

  # 繪圖
  plt.stem(x_index, y_m, markerfmt='o')
  plt.xlabel('x')
  plt.ylabel('Magnitude')
  plt.show()

圖形結果與討論
- 結果
- 問題一：
  - x座標應該對應到頻率
- 問題二：
  - y座標振幅不合理
  - 訊號的振幅最大只有$4$，最小只有$-3$，但是頻譜分析出的振幅超過$7$，明顯不合理
- 問題三：
  - 頻譜會有正負頻率特性，此處並沒有

5.2.2 修正解法 ⇒ 修正x座標、正負頻率性質與振幅呈現

程式碼

  import numpy as np
  from numpy.fft import fft, fftshift, fftfreq
  import matplotlib.pyplot as plt

  # 訊號點
  x = np.array([4.0, 2.0, 1.0, -3.0, 1.5])
  y = fft(x)
  print(y)

  # fftfreq: 透過資料點數量與取樣頻率，產生各點頻率
  # fftshift: 平移x軸座標點，位移負頻率的振幅
  f = fftshift(fftfreq(5, 0.001))
  X = fftshift(fft(x))
  Xm = abs(X)
  # 標準化個頻率的振幅
  # 除以資料點數量
  Xm = Xm/len(x)

  plt.stem(f, Xm, markerfmt='o')
  plt.xlabel('f')
  plt.ylabel('Magnitude')

  plt.show()

圖形結果與討論
- 結果
  
  修正x座標
- 因為題目的資料是1000 Hz，可以透過 numpy 的 fftfreq產生相對應的頻率座標
  共軛性質呈現
- 共軛性質需要將FFT的結果做位移，可以透過 numpy 的 fftshift 完成
  振幅呈現修改
- ~~前面講解 $cos$ 函數的頻譜分析時有看到，頻率若為 $\omega$，振幅的大小為 $A\pi$，則頻率若為 $f$，振幅的大小為 $\frac A2$。~~
- 但這是因為 $cos$ 函數此時只有一個頻率，如果有N個頻率點時，每個頻率點的振幅會變成是$\frac{A_k}2N$，因此在畫FFT的Magnitude Spectrum時，要把振幅除以$N$，才會是正確的振幅
- 在畫FFT的Magnitude Spectrum時，要把振幅除以$N$，才會是正確的振幅
- 除以 $N$ 的原因，簡單講是Fourier Series到Fourier Transform的數學式定義上的些微差異造成的，詳情參考 (8. 傅立葉番外篇)(待續)
另外一個問題：
- 明明就是有在變動的訊號，為什麼在0 Hz的地方可以看到那麼大的分量

5.2.3 修正解法 ⇒ 去除直流成分，將 0 Hz處的振幅去除

程式碼

  import numpy as np
  from numpy.fft import fft, fftshift, fftfreq
  import matplotlib.pyplot as plt

  # 訊號點
  x = np.array([4.0, 2.0, 1.0, -3.0, 1.5])
  # 去除直流成分
  x = x - np.mean(x)
  y = fft(x)
  print(y)

  # fftfreq: 透過資料點數量與取樣頻率，產生各點頻率
  # fftshift: 平移x軸座標點，位移負頻率的振幅
  f = fftshift(fftfreq(5, 0.001))
  X = fftshift(fft(x))
  Xm = abs(X)
  # 標準化個頻率的振幅
  # 除以資料點數量
  Xm = Xm/len(x)

  plt.stem(f, Xm, markerfmt='o')
  plt.xlabel('f')
  plt.ylabel('Magnitude')

  plt.show()

圖形結果與討論
- 結果
  
  在實務上如果要去除直流成分都是將全部訊號減去平均值
  這個作法的原因可以有兩種解釋
- 有查到一個說法：一段訊號的平均值就是直流成分的 "估計量” or “不偏估計量”，所以平均值可以代表直流成分
- 另外一個解釋法可以參考 (8. 傅立葉番外篇) (待續)

範例 — 諧波訊號頻譜分析

若諧波定義為：
$x(t)=10+10\cos\left(2\pi\cdot100\cdot t\right)+2\cos\left(2\mathrm\pi\cdot200\cdot\mathrm t\right)$
其中包含直流分量 $A_0=10$
交流分量由兩個不同振福、不同頻率的弦波構成
取樣頻率 $f_s=1000Hz$
求傅立葉頻譜分析

解法

程式碼

  import numpy as np
  from numpy.fft import fft, fftshift, fftfreq
  import matplotlib.pyplot as plt

  t = np.linspace( 0, 1, 1000, endpoint = False )
  x = 10 + 10 * np.cos( 2 * np.pi * 100 * t ) + 2 * np.cos( 2 * np.pi * 200 * t )

  f = fftshift( fftfreq( 1000, 0.001 ) )
  X = fftshift( fft( x ) ) 
  Xm = abs( X )/len(X)

  plt.plot( f, Xm )
  plt.xlabel( 'f (Hz)' )
  plt.ylabel( 'Magnitude' )

  plt.show( )

圖形結果

5.3 傅立葉頻譜分析圖特性總結

傅立葉頻譜分析會呈現正負頻率的結果
- 有正、負頻率，且頻譜資料點數量各一半
- 在實務上只要看正頻率即可
傅立葉頻譜分析可以看到的頻率成分只有取樣頻率的一半
- 符合Nyquist-Shannon 取樣定理
- 使用 $f_s$ 的取樣頻率，代表原始訊號的頻率最大只有 $\frac{f_s}2$，因此傅立葉頻譜分析只能看到取樣頻率的一半的頻率成分
資料點數量就是傅立葉頻譜的解析度
- 資料點太少，轉換成傅立葉頻譜時頻率的跨度就會很大
- 增加資料點的數量，或是降低取樣頻率，都可以增加傅立葉頻譜的解析度