《数字信号处理》课程实验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)