[转载] 傅里叶分析之掐死教程

[信号与频谱]  2.9 傅里叶变换

[信号与频谱]  2.10 离散傅里叶变换

[信号与频谱] 快速傅里叶变换（FFT）算法及蝶形运算推导

迭代版FFT

import numpy as np
import matplotlib.pyplot as plt

# 辅助函数：将输入补零到2的幂次长度
defpad_to_power2(x):
    N = len(x)
    next_pow2 = 1
    L = 0
while next_pow2 < N:
        next_pow2 <<= 1# 等价于*2
        L += 1
return np.pad(x, (0, next_pow2 - N), mode='constant'), next_pow2, L

# 将一个数的二进制位逆序重排
defreverse_bits(n, L):
    result = 0
for i in range(L):
# 获取第i位的值
        bit = (n >> i) & 1
# 逆序放到对应的位置
        result |= (bit << (L - 1 - i))
return result

# 位反转置换方法二
defreverse_bits_method2(j, N):
    bit = N >> 1
while j & bit:
        j ^= bit
        bit >>= 1
    j ^= bit
return j

# 迭代版通过位逆序重排和蝶形运算实现原地计算，是实际应用的主流方式
deffft_iterative(x):
    x, N, L = pad_to_power2(x)  # 确保长度为2的幂次
    x = x.astype(np.complex128)

# 步骤1：位逆序重排（Bit-reversal permutation）
# 例如N=8：下标0(000)→0, 1(001)→4, 2(010)→2, 3(011)→6, 4(100)→1, 5(101)→5, 6(110)→3, 7(111)→7
# for i in range(1, N):
#     j = reverse_bits(i, L)
#     if(i < j):
#         x[i], x[j] = x[j], x[i]

# 位逆序方法二
    j = 0
for i in range(1, N):
        j = reverse_bits_method2(j, N)
if(i < j):
            x[i], x[j] = x[j], x[i]

# 步骤2：蝶形运算（按层计算）
# 层：从2点DFT开始，到N点DFT，步长为2^n
    n = 2
# 每一级的DFT大小m, 例如第一级n=2, 进行2点DFT，第二级n=4,进行4点DFT
while(n <= N):
# 循环遍历当前级中每个大小为n的DFT块，步长为n，所以每次循环处理一个DFT块
for i in range(0, N, n):
# 循环处理每个DFT块内的蝶形运算，每个大小为n的DFT块可以分解为n/2个蝶形运算
for k in range(n//2):   # k = 0, 1, ..., n/2 -1
                W_N = np.exp(-2j * np.pi * k / n)  # 旋转因子
                x0 = x[i + k]
                x1 = x[i + k + n//2] * W_N
                x[i + k] = x0 + x1
                x[i + k + n//2] = x0 - x1

        n <<= 1# 下一级DFT大小

return x

# 测试信号：正弦波 + 噪声
fs = 1000# 采样率
t = np.linspace(0, 1, fs, endpoint=False)  # 1秒，1000个采样点
x = 3 * np.sin(2 * np.pi * 50 * t) + 2 * np.sin(2 * np.pi * 120 * t) + 0.5 * np.random.randn(len(t))
# x = [8, 2, 3, -3, 5, 6, 4, 7]
#x = np.array(x, dtype=complex)

# 计算FFT
x_padded, N, L = pad_to_power2(x)
X_recursive = fft_recursive(x_padded)
X_iterative = fft_iterative(x)
X_numpy = np.fft.fft(x, N)  # NumPy官方FFT

# 验证误差（应接近0）
print("递归版 vs NumPy 最大误差：", np.max(np.abs(X_recursive - X_numpy)))
print("迭代版 vs NumPy 最大误差：", np.max(np.abs(X_iterative - X_numpy)))
# 绘制频谱
freq = np.fft.fftfreq(N, 1/fs)  # 频率轴
mask = freq >= 0# 只显示正频率
plt.figure(figsize=(12, 6))
plt.plot(freq[mask], 2/N * np.abs(X_iterative[mask]), label='FFT iterative')
plt.plot(freq[mask], 2/N * np.abs(X_numpy[mask]), '--', label='NumPy FFT', alpha=0.7)
plt.xlabel('Freq (Hz)')
plt.ylabel('Amp')
plt.title('FFT (50Hz + 120Hz sin())')
plt.legend()
plt.grid(True)
plt.show()

递归版FFT

# 递归实现
deffft_recursive(x):
    N = len(x)
# 基例：N=1时，DFT就是自身
if N == 1:
return x

# 分解为奇偶序列
    even = fft_recursive(x[::2])  # 偶下标：0,2,4...
    odd = fft_recursive(x[1::2])   # 奇下标：1,3,5...

# 初始化结果
    X = np.zeros(N, dtype=np.complex128)
# 计算旋转因子 W_N^k = e^(-j*2πk/N)
for k in range(N//2):
        W = np.exp(-2j * np.pi * k / N)
        X[k] = even[k] + W * odd[k]
        X[k + N//2] = even[k] - W * odd[k]
return X