《数字信号处理》课程实验1 – FFT的实现

一、按时间抽选的基-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 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

三、反变换IFFT的实现

由于IDFT公式为：
(x(n)={ m 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} { m DFT}[X^* (k)])
那么：
(x(n)=frac {1}{N} {{ m 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)

四、程序测试

按照要求，分别使用如下三种信号测试算法。使用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程序：

五、补充说明

本程序的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 = '[
'
  modder = 0
  for cplx in seq:
    toprint += str(cplx)
    toprint += '		' if modder != 3 else '
'
    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)