一、按時間抽選的基-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)={\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)
轉載請註明出處，本文鏈接：https://www.uj5u.com/qita/77811.html
標籤：其他
上一篇：CF750D New Year and Fireworks
下一篇：單鏈表操作（帶頭節點）