《数字信号处理》课程实验1 – FFT的实现
一、按时间抽选的基-2 FFT实现原理
观察DIT（基2）FFT的流图（N点，N为2的幂次），可以总结出如下规律：
（1）共有\(L=\log_2N\)级蝶形运算；
（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)={\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)
四、程序测试
按照要求，分别使用如下三种信号测试算法。使用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 = '[\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)
