《數字信號處理》課程實驗1 – FFT的實現

1、按時間抽選的基-2 FFT實現原理

觀察DIT（基2）FFT的流圖（N點，N爲2的冪次），能夠總結出以下規律：
（1）共有\(L=\log_2⁡N\)級蝶形運算；
（2）輸入倒位序，輸出天然順序；
（3）第\(m\)級（\(m\)從1開始，下同）蝶形結對偶結點距離爲\(2^{m-1}\)；
（4）第\(m\)級蝶形結計算公式：
\(X_m (k)=X_{m-1} (k)+X_{m-1 } (k+2^{m-1} ) W_N^r\)
\(X_m (k+2^{m-1} )=X_{m-1} (k)-X_{m-1} (k+2^{m-1} ) W_N^r\)
\(m=1,2,…,L\)
（5）第\(m\)級不一樣旋轉因子數爲\(2^{m-1}\)；
（6）第\(m\)級每一種旋轉因子關聯\(2^{L-m}\)個蝶形結；
（7）第\(m\)級一種旋轉因子的兩次相鄰出現間隔\(2^m\)個位置；
（8）第\(m\)級第\(i\)種（\(i\)從0開始）旋轉因子\(W_N^r\)的指數\(r\)爲\(2^{L-m}i\)；
（7）每一個蝶形結運算完成後，輸出的兩節點值能夠直接放到原輸入的兩節點的存儲器中（原位運算）。

2、按時間抽選的基-2 FFT實現細節

依據如上觀察，使用Python語言編寫下列相關程序：

（1）倒位變址運算

天然序排列的二進制數，其下一個數總比前一個大1。而倒序二進制數的下一個數，是前一個數最高位加1，而後由高位向低位進位獲得的。使用Rader算法，能夠方便地計算倒位序。
Rader算法利用一個遞推關係——若是已知第\(k\)步倒位序爲\(J\)，則第\(k+1\)步倒位序計算方法爲：從\(J\)最高位看向最低位，爲1則變0；爲0則變1並馬上退出，即獲得下一步倒位序。因爲已知遞推起點第1步的序號0的倒位序也爲0，故可使用該算法求出全部倒位序。
進行倒位變址運算時，爲避免重複調換，設輸入爲\(x(n)\)，倒位序後爲\(x(m)\)，當且僅當\(m>n\)時進行對調。
下面是相關代碼：

new_index = [0]
J = 0 # J爲倒位序
for i in range(N - 1): # i爲當前數
  mask = N // 2
  while mask <= J: # J的最高位爲1
    J -= mask # J的最高位置0
    mask = mask >> 1 # 準備檢測下一位
  J += mask # 首個非1位置1
  new_index.append(int(J))
for i in range(N):
  if new_index[i] <= i:
    continue # 無需對調
  seq[i], seq[new_index[i]] = seq[new_index[i]], seq[i] # 交換

（2）旋轉因子的計算

利用歐拉公式可知：
\(W_N^k=\cos⁡(2kπ/N)-j \sin⁡(2kπ/N)\)
在\(k=0,1,…,N-1\)範圍內循環一次，便可計算全部旋轉因子。
一種優化策略是利用以下遞推關係來加速計算，遞推起點\(W_N^0=1\)：
\(W_N^{(k+1)}=W_N^k \exp⁡(-j 2π/N)\)
相關代碼以下：

WNk = []
two_pi_div_N = 2 * PI / N # 避免屢次計算
for k in range(N // 2):
  # WN^k = cos(2kPI/N) - j sin(2kPI/N)
  WNk.append(Complex(math.cos(two_pi_div_N * k), math.sin(two_pi_div_N * -k)))

（3）蝶形結按層運算

以旋轉因子的種類爲循環標準，每一輪就算掉該種旋轉因子對應的\(2^{L-m}\)個蝶形結。結合觀察（3）~（8），相關代碼以下：

L = int(math.log2(N)) # 蝶形結層數
for m in range(1, L + 1): # m爲當前層數，從1開始
  distance = 2 ** (m - 1) # 對偶結點距離，也是該層不一樣旋轉因子的數量
  for k in range(distance): # 以結合的旋轉因子爲循環標準，每一輪就算掉該旋轉因子對應的2^(L-m)個結
    r = k * 2 ** (L - m) # 該旋轉因子對應的r
    for j in range(k, N, 2 ** m): # 2^m爲每組旋轉因子對應的分組的下標差
      right = seq[j + distance] * WNk[r]
      t1 = seq[j] + right; t2 = seq[j] - right
      seq[j] = t1; seq[j + distance] = t2

3、反變換IFFT的實現

因爲IDFT公式爲：
\(x(n)={\rm IDFT}[X(k)]=\frac {1}{N} ∑_{k=0}^{N-1} X(k) W_N^{-nk}\)
將該式取共軛：
\(x^* (n)=\frac {1}{N} [∑_{k=0}^{N-1}X(k) W_N^{-nk} ]^*=\frac {1}{N} ∑_{k=0}^{N-1}[X^* (k) W_N^{nk} ]=\frac {1}{N} {\rm DFT}[X^* (k)]\)
那麼：
\(x(n)=\frac {1}{N} {{\rm DFT}[X^* (k)]}^*\)
這意味着IFFT能夠共用FFT程序。只要將\(X(k)\)取共軛後作FFT，結果再取共軛併除以\(N\)即完成了IFFT。相關代碼以下：

def ifft(seq: list):
  # 檢查是否爲2^L點序列
  N = len(seq)
  if int(math.log2(N)) - math.log2(N) != 0:
    raise ValueError('[ifft] Not 2^L long sequence.')
  # 先對X(k)取共軛
  seq = list(map(lambda x : x.conjugate(), seq))
  # 再次利用FFT
  seq = fft_dit2(seq)
  # 再取共軛
  seq = map(lambda x : x.conjugate(), seq)
  # 最後除以N
  seq = map(lambda x : x * Complex(1 / N, 0), seq)
  return list(seq)

4、程序測試

按照要求，分別使用以下三種信號測試算法。使用matplotlib庫做圖，使用numpy庫的FFT結果與自行撰寫的FFT結果進行對照。
（1）正弦序列
產生\([-π,π]\)上的16點均勻採樣，計算相應的\(\sin\)函數值，進行FFT，結果以下，正確無誤：

繪製序列及其幅度頻譜圖，同時繪製IFFT結果以測試IFFT程序：

（2）三角波序列
從0開始向正方向，以\(π/8\)爲間隔，產生三角波的16點\(x\)值，並按\([0, 0.5, 1, 0.5, 0, 0.5, 1…]\)規律賦\(y\)值。做FFT，結果以下，正確無誤：

繪製序列及其幅度頻譜圖，同時繪製IFFT結果以測試IFFT程序：

（3）方波序列
產生\([-π,π]\)上的32點均勻採樣做爲方波的\(x\)值，並按\([0, 0, 0, 0, 1, 1, 1, 1, 0, 0, 0, 0, 1, 1, 1, 1…]\)規律賦\(y\)值（佔空比50%，週期8點）。做FFT，結果以下，正確無誤：

繪製序列及其幅度頻譜圖，同時繪製IFFT結果以測試IFFT程序：

5、補充說明

本程序的fft_dit2函數默認對輸入的\(N\)長度序列做\(N\)點FFT，即採樣點數與FFT點數相同。但有時候，須要對\(N\)長度序列做\(L(L>N)\)點FFT。爲此，程序提供了一個實用函數add_zero_to_length來把\(N\)點序列末尾補0到指定長度以支持上述FFT運算。
使用一個矩形窗序列\(R_8 (n)\)，對其做16點FFT，測試程序正確性，程序正確無誤，以下所示：

該段測試用的代碼因爲不是實驗要求，故在fft_demo.py中做註釋處理。
除了補0函數外，程序還提供了下列這些實用函數來幫助更方便地使用FFT：

附錄：所有代碼

fft.py

# coding: utf-8

# 《數字信號處理》課程實驗
# FFT與IFFT實現（一維）
# 09017227 卓旭

import math

'''
  自定義複數類
'''
class Complex:

  def __init__(self, re, im):
    self.set(re, im)

  def set(self, re, im):
    self._re = re
    self._im = im

  def get(self, what):
    return self._re if what == 're' else self._im

  def __add__(self, other):
    return Complex(
      self._re + other._re,
      self._im + other._im
    )

  def __sub__(self, other):
    return Complex(
      self._re - other._re,
      self._im - other._im
    )

  def __mul__(self, other):
    return Complex(
      (self._re * other._re) - (self._im * other._im),
      (self._re * other._im) + (self._im * other._re)
    )

  def conjugate(self):
    return Complex(
      self._re, -self._im
    )

  def __abs__(self):
    return math.sqrt(self._re ** 2 + self._im ** 2)

  def __str__(self):
    return (str(round(self._re, 3)) + ' + j' + str(round(self._im, 3))) if (self._im >= 0) \
      else (str(round(self._re, 3)) + ' - j' + str(round(-self._im, 3)))

PI=math.pi

'''
  按時間抽選的基-2快速傅里葉變換（n點）
  須要傳入list<Complex>
'''
def fft_dit2(seq: list):
  # 檢查是否爲2^L點FFT
  N = len(seq)
  if int(math.log2(N)) - math.log2(N) != 0:
    raise ValueError('[fft_dit2] Not 2^L long sequence.')

  # 輸入數據倒位序處理
  new_index = [0]
  J = 0 # J爲倒位序
  for i in range(N - 1): # i爲當前數
    mask = N // 2
    while mask <= J: # J的最高位爲1
      J -= mask # J的最高位置0
      mask = mask >> 1 # 準備檢測下一位
    J += mask # 首個非1位置1
    new_index.append(int(J))
  for i in range(N):
    if new_index[i] <= i:
      continue # 無需對調
    seq[i], seq[new_index[i]] = seq[new_index[i]], seq[i] # 交換

  # 計算全部須要的旋轉因子WN^k（k在0~N/2-1)
  # 一種優化策略是使用遞推式WN(k+1) = WN(k) * e^(-j 2PI/N)計算
  WNk = []
  two_pi_div_N = 2 * PI / N # 避免屢次計算
  for k in range(N // 2):
    # WN^k = cos(2kPI/N) - j sin(2kPI/N)
    WNk.append(Complex(math.cos(two_pi_div_N * k), math.sin(two_pi_div_N * -k)))

  # 蝶形運算
  L = int(math.log2(N)) # 蝶形結層數
  for m in range(1, L + 1): # m爲當前層數，從1開始
    # 見課本P219表4.1
    distance = 2 ** (m - 1) # 對偶結點距離，也是該層不一樣旋轉因子的數量
    for k in range(distance): # 以結合的旋轉因子爲循環標準，每一輪就算掉該旋轉因子對應的2^(L-m)個結
      r = k * 2 ** (L - m) # 該旋轉因子對應的r
      for j in range(k, N, 2 ** m): # 2^m爲每組旋轉因子對應的分組的下標差
        right = seq[j + distance] * WNk[r]
        t1 = seq[j] + right; t2 = seq[j] - right
        seq[j] = t1; seq[j + distance] = t2

  return seq

'''
  快速傅里葉變換的反變換
  須要傳入list<Complex>
'''
def ifft(seq: list):
  # 檢查是否爲2^L點序列
  N = len(seq)
  if int(math.log2(N)) - math.log2(N) != 0:
    raise ValueError('[ifft] Not 2^L long sequence.')

  # 先對X(k)取共軛
  seq = list(map(lambda x : x.conjugate(), seq))
  # 再次利用FFT
  seq = fft_dit2(seq)
  # 再取共軛
  seq = map(lambda x : x.conjugate(), seq)
  # 最後除以N
  seq = map(lambda x : x * Complex(1 / N, 0), seq)

  return list(seq)

'''
  實用函數，將實序列轉化爲list<Complex>
'''
def convert_to_complex(seq: list):
  return list(map(lambda x : Complex(x, 0), seq))

'''
  實用函數，將list<Complex>轉化爲實序列（丟棄虛部）
'''
def convert_to_real(seq: list):
  return list(map(lambda x : x.get('re'), seq))

'''
  實用函數，獲取Complex的幅度值
'''
def convert_to_amplitude(seq: list):
  return list(map(lambda x: math.sqrt(x.get('re') ** 2 + x.get('im') ** 2), seq))

'''
  實用函數，獲取Complex的相位值
'''
def convert_to_phase(seq: list):
  return list(map(lambda x: math.atan2(x.get('im'), x.get('re')), seq))

'''
  實用函數，打印FFT結果
'''
def print_fft_result(seq: list):
  toprint = '[\n'
  modder = 0
  for cplx in seq:
    toprint += str(cplx)
    toprint += '\t\t' if modder != 3 else '\n'
    modder += 1
    modder %= 4
  toprint += ']'
  return toprint

'''
  實用函數，給實序列補0到指定長度，可用於採樣點數小於FFT點數
'''
def add_zero_to_length(seq: list, n: int):
  if len(seq) == n:
    return seq
  # 若是點數不足，把seq補到n點
  if len(seq) > n:
    raise ValueError('[add_zero_to_length] n < len(seq).')
  if len(seq) < n:
    res = [*seq]
    while (len(res) < n):
      res.append(0)
  return res

---

fft_demo.py

# coding: utf-8

# 《數字信號處理》課程實驗
# FFT與IFFT實測應用（做圖）
# 09017227 卓旭

import matplotlib.pyplot as plt
import numpy as np

from fft import *

PI=math.pi

def ft_test(name, xs, ys):
  # 產生測試點
  y_points = ys
  # 繪製原圖
  plt.subplots_adjust(hspace=0.7)
  plt.subplot(311)
  plt.title('Original Singal: ' + name)
  plt.plot(xs, y_points, '.')
  # 繪製FFT結果（幅度譜）
  fft_res = fft_dit2(convert_to_complex(y_points))
  plt.subplot(312)
  plt.title('FFT Result (Amplitude Spectrum)')
  fft_res_amp = convert_to_amplitude(fft_res)
  plt.plot(xs, fft_res_amp, '.')
  max_height = np.max(fft_res_amp)
  for (idx, val) in enumerate(xs):
    plt.axvline(val, 0, fft_res_amp[idx] / max_height) # 繪製豎線
  # 繪製IFFT
  ifft_res = convert_to_real(ifft(fft_res))
  plt.subplot(313)
  plt.title('IFFT Result')
  plt.plot(xs, ifft_res, '.')
  plt.show()

if __name__ == '__main__':
  # 正弦函數測試
  xs = np.linspace(-PI, PI, 16)
  ys = [math.sin(x) for x in xs]
  print("========正弦函數========")
  print("np.fft標準結果：")
  print(np.fft.fft(ys))
  print("個人FFT結果：")
  print(print_fft_result(fft_dit2(convert_to_complex(ys))))
  ft_test('sin(x)', xs, ys)
  print("========三角波========")
  xs = [i * PI / 8 for i in range(16)]
  ys = [0, 0.5, 1, 0.5, 0, 0.5, 1, 0.5, 0, 0.5, 1, 0.5, 0, 0.5, 1, 0.5]
  print("np.fft標準結果：")
  print(np.fft.fft(ys))
  print("個人FFT結果：")
  print(print_fft_result(fft_dit2(convert_to_complex(ys))))
  ft_test('Tri(x)', xs, ys)
  print("========方波========")
  xs = np.linspace(-PI, PI, 32)
  ys = [-1,-1,-1,-1, 1, 1, 1, 1] * 4
  print("np.fft標準結果：")
  print(np.fft.fft(ys))
  print("個人FFT結果：")
  print(print_fft_result(fft_dit2(convert_to_complex(ys))))
  ft_test('Square(x)', xs, ys)
  # print("========R_8========")
  # xs = range(16)
  # ys = [1, 1] * 4
  # ys = add_zero_to_length(ys, 16)
  # print("np.fft標準結果：")
  # print(np.fft.fft(ys, 16))
  # print("個人FFT結果：")
  # print(print_fft_result(fft_dit2(convert_to_complex(ys))))
  # ft_test('R_8(x), 16 dot', xs, ys)