跳转到内容

数字信号处理/快速傅里叶变换 (FFT) 算法

来自维基教科书,开放的书籍,开放的世界

快速傅里叶变换 (FFT) 是一种用于计算离散傅里叶变换 (DFT) 的算法,使用更少的乘法次数。

库利-图基算法

[编辑 | 编辑源代码]

库利-图基算法是最常见的快速傅里叶变换 (FFT) 算法。它将较大的 DFT 分解成较小的 DFT。库利-图基算法最著名的应用是将 N 点变换分解成两个 N/2 点变换,因此它仅限于 2 的幂大小。

该算法可以通过两种不同的方式实现,使用所谓的

  • 时间抽取 (DIT)
  • 频率抽取 (DIF)

电路图

[编辑 | 编辑源代码]

下图表示一个 16 点 DIF FFT 电路图

16-point DIF FFT circuit
16 点 DIF FFT 电路

在图中,系数由

其中 N 是 FFT 的点数。它们在单位圆的下半部分均匀分布,并且 .

测试脚本

[编辑 | 编辑源代码]

以下 python 脚本允许测试上面说明的算法

#!/usr/bin/python3
import math
import cmath

# ------------------------------------------------------------------------------
# Parameters
#
coefficientBitNb = 8    # 0 for working with reals

# ------------------------------------------------------------------------------
# Functions
#
# ..............................................................................

def  bit_reversed(n, bit_nb):
    binary_number = bin(n)
    reverse_binary_number = binary_number[-1:1:-1]
    reverse_binary_number = reverse_binary_number + \
        (bit_nb - len(reverse_binary_number))*'0'
    reverse_number = int(reverse_binary_number, bit_nb-1)

    return(reverse_number)

# ..............................................................................

def  fft(x):
    j = complex(0, 1)
                                                                         # sizes
    point_nb = len(x)
    stage_nb = int(math.log2(point_nb))
                                                                       # vectors
    stored = x.copy()
    for index in range(point_nb):
        stored[index] = complex(stored[index], 0)
    calculated = [complex(0, 0)] * point_nb
                                                                  # coefficients
    coefficients = [complex(0, 0)] * (point_nb//2)
    for index in range(len(coefficients)):
        coefficients[index] = cmath.exp(-j * index/point_nb * 2*cmath.pi)
        coefficients[index] = coefficients[index] * 2**coefficientBitNb
                                                                          # loop
    for stage_index in range(stage_nb):
        # print([stored[i].real for i in range(point_nb)])
        index_offset = 2**(stage_nb-stage_index-1)
                                                           # butterfly additions
        for vector_index in range(point_nb):
            isEven = (vector_index // index_offset) % 2 == 0
            if isEven:
                operand_a = stored[vector_index]
                operand_b = stored[vector_index + index_offset]
            else:
                operand_a = - stored[vector_index]
                operand_b = stored[vector_index - index_offset]
            calculated[vector_index] = operand_a + operand_b
                                                   # coefficient multiplications
        for vector_index in range(point_nb):
            isEven = (vector_index // index_offset) % 2 == 0
            if not isEven:
                coefficient_index = (vector_index % index_offset) \
                    * (stage_index+1)
                # print(                                                  \
                #     "(%d, %d) -> %d"                                    \
                #     % (stage_index, vector_index, coefficient_index)    \
                # )
                calculated[vector_index] = \
                    coefficients[coefficient_index] * calculated[vector_index] \
                    / 2**coefficientBitNb
                if coefficientBitNb > 0:
                    calculated[vector_index] = complex(               \
                        math.floor(calculated[vector_index].real),    \
                        math.floor(calculated[vector_index].imag)     \
                    )
                                                                       # storage
        stored = calculated.copy()
                                                               # reorder results
    for index in range(point_nb):
        calculated[bit_reversed(index, stage_nb)] = stored[index]

    return calculated

# ------------------------------------------------------------------------------
# Main program
#
source = [0.0, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0]
transformed = fft(source)
amplitude = [0.0] * len(source)
amplitude_print = ''
for index in range(len(transformed)):
    amplitude[index] = abs(transformed[index])
    amplitude_print = amplitude_print + "%5.3f " % amplitude[index]
print()
print(amplitude_print)

它有一个特殊的参数,coefficientBitNb,它允许确定使用仅整数的电路的计算结果。

华夏公益教科书