这两天一直在研究FFT算法,研究了两天总算整明白了快速傅里叶变换算法,今天就写一篇文章来总结一下。
FFT和IFFT的Python语言实现源代码
直接把我用了一个晚上写好的快速傅里叶变换和快速傅里叶逆变换的Python语言代码贴出,关键部分有注释,里面只用到了Python标准库cmath库,因为要计算cos、sin函数的值。直接复制到自己的Python程序中就可以直接使用了。
"""
@Author: Sam
@Function: Fast Fourier Transform
@Time: 2020.02.22 16:00
"""
from cmath import sin, cos, pi
class FFT_pack():
def __init__(self, _list=[], N=0): # _list 是传入的待计算的离散序列,N是序列采样点数,对于本方法,点数必须是2^n才可以得到正确结果
self.list = _list # 初始化数据
self.N = N
self.total_m = 0 # 序列的总层数
self._reverse_list = [] # 位倒序列表
self.output = [] # 计算结果存储列表
self._W = [] # 系数因子列表
for _ in range(len(self.list)):
self._reverse_list.append(self.list[self._reverse_pos(_)])
self.output = self._reverse_list.copy()
for _ in range(self.N):
self._W.append((cos(2 * pi / N) - sin(2 * pi / N) * 1j) ** _) # 提前计算W值,降低算法复杂度
def _reverse_pos(self, num) -> int: # 得到位倒序后的索引
out = 0
bits = 0
_i = self.N
data = num
while (_i != 0):
_i = _i // 2
bits += 1
for i in range(bits - 1):
out = out << 1
out |= (data >> i) & 1
self.total_m = bits - 1
return out
def FFT(self, _list, N, abs=True) -> list: # 计算给定序列的傅里叶变换结果,返回一个列表,结果是没有经过归一化处理的
"""参数abs=True表示输出结果是否取得绝对值"""
self.__init__(_list, N)
for m in range(self.total_m):
_split = self.N // 2 ** (m + 1)
num_each = self.N // _split
for _ in range(_split):
for __ in range(num_each // 2):
temp = self.output[_ * num_each + __]
temp2 = self.output[_ * num_each + __ + num_each // 2] * self._W[__ * 2 ** (self.total_m - m - 1)]
self.output[_ * num_each + __] = (temp + temp2)
self.output[_ * num_each + __ + num_each // 2] = (temp - temp2)
if abs == True:
for _ in range(len(self.output)):
self.output[_] = self.output[_].__abs__()
return self.output
def FFT_normalized(self, _list, N) -> list: # 计算给定序列的傅里叶变换结果,返回一个列表,结果经过归一化处理
self.FFT(_list, N)
max = 0 # 存储元素最大值
for _ in range(len(self.output)):
if max < self.output[_]:
max = self.output[_]
for _ in range(len(self.output)):
self.output[_] /= max
return self.output
def IFFT(self, _list, N) -> list: # 计算给定序列的傅里叶逆变换结果,返回一个列表
self.__init__(_list, N)
for _ in range(self.N):
self._W[_] = (cos(2 * pi / N) - sin(2 * pi / N) * 1j) ** (-_)
for m in range(self.total_m):
_split = self.N // 2 ** (m + 1)
num_each = self.N // _split
for _ in range(_split):
for __ in range(num_each // 2):
temp = self.output[_ * num_each + __]
temp2 = self.output[_ * num_each + __ + num_each // 2] * self._W[__ * 2 ** (self.total_m - m - 1)]
self.output[_ * num_each + __] = (temp + temp2)
self.output[_ * num_each + __ + num_each // 2] = (temp - temp2)
for _ in range(self.N): # 根据IFFT计算公式对所有计算列表中的元素进行*1/N的操作
self.output[_] /= self.N
self.output[_] = self.output[_].__abs__()
return self.output
def DFT(self, _list, N) -> list: # 计算给定序列的离散傅里叶变换结果,算法复杂度较大,返回一个列表,结果没有经过归一化处理
self.__init__(_list, N)
origin = self.list.copy()
for i in range(self.N):
temp = 0
for j in range(self.N):
temp += origin[j] * (((cos(2 * pi / self.N) - sin(2 * pi / self.N) * 1j)) ** (i * j))
self.output[i] = temp.__abs__()
return self.output
if __name__ == '__main__':
list = [1, 2, 3, 4, 5, 6, 7, 8, 0, 0, 0, 0, 0, 0, 0, 0]
a = FFT_pack().FFT(list, 16, False)
print(a)
下面对源代码的一些关键部分进行解释:
我直接将计算方法的调用封装成了一个类,这样可以很方便进行扩展和调用。在FFT_pack()
这个类之中,我定义了下面的几种方法
1.构造函数__init__()
,目的是初始化需要进行FFT变换的序列,对采样点数进行赋值、计算分治算法需要进行分组的层数,计算旋转因子的系数列表等。
2._reverse_pos()
方法,是为了获得位倒序后的顺序,这个方法是一个不需要外部调用的方法。
3.FFT()
方法,顾名思义,是计算给定序列的快速傅里叶变换结果。里面涉及到四个参数,在实际调用的时候需要传入三个参数-list
N
abs
。_list
参数是需要进行FFT运算的列表,N
参数是参加计算的点的个数,注意: 这里面N的值必须是2的m次方,例如N可以是2、4、8、16、32、……1024、2048、4096等这样的数,如果填入了别的数,无法得到正确的计算结果。abs
参数是缺省值为True的参数,当abs赋值为True或者使用缺省值的时候,返回的FFT运算结果是取绝对值以后的序列,当abs参数赋值为False时,返回的FFT运算结果就是一个复数,含有实部和虚部的值。
4.FFT_normalized()
方法,用法与FFT()
方法类似,只不过没有abs参数,方法的返回值是经过归一化处理的FFT变换结果。
5.IFFT()
方法,返回给定序列的快速傅里叶逆变换的序列。
6.DFT()
方法,返回给定序列的离散傅里叶变换序列,返回结果是经过取绝对值运算的。这个DFT()方法主要是用来与FFT方法的运算性能进行对比的,实际使用中还是使用FFT方法。
下面总结一下FFT算法
对于每一个鲽形运算最小单元
对于不同旋转因子的种类和选取方法
将原始序列进行位倒序
IFFT算法
根据IDFT算法的定义式,可以知道逆变换与正变换的差别非常小,旋转因子表每一项的方次取负值,运算之后的序列整个乘以1/N即可。