快速傅里叶变换（FFT）算法的程序原理解析与实现指南

F(ω) = ∫[-∞ to ∞] f(t)e^(-jωt)dt
f(t) = (1/2π) ∫[-∞ to ∞] F(ω)e^(jωt)dω

X[k] = Σ[n=0 to N-1] x[n] * e^(-j2πkn/N), k = 0, 1, ..., N-1

X[k] = Σ[n=0 to N-1] x[n] * W_N^(kn)
     = Σ[r=0 to N/2-1] x[2r] * W_N^(k(2r)) + Σ[r=0 to N/2-1] x[2r+1] * W_N^(k(2r+1))
     = Σ[r=0 to N/2-1] x[2r] * W_(N/2)^(kr) + W_N^k * Σ[r=0 to N/2-1] x[2r+1] * W_(N/2)^(kr)

import numpy as np
import cmath
 
def fft(x):
    """
    快速傅里叶变换实现
    :param x: 输入序列（长度必须为2的幂次）
    :return: 频域系数
    """
    n = len(x)
    
    # 基本情况：单点DFT
    if n <= 1:
        return x
    
    # 确保长度为偶数
    if n % 2 > 0:
        raise ValueError("输入长度必须是2的幂次")
    
    # 分治：分离偶数和奇数索引
    even = fft(x[0::2])
    odd = fft(x[1::2])
    
    # 组合结果
    result = [0] * n
    for k in range(n // 2):
        # 旋转因子
        twiddle = cmath.exp(-2j * cmath.pi * k / n)
        result[k] = even[k] + twiddle * odd[k]
        result[k + n // 2] = even[k] - twiddle * odd[k]
    
    return result
 
# 测试代码
def test_fft():
    # 生成测试信号：两个正弦波的叠加
    sample_rate = 1000  # 采样率
    duration = 1.0      # 持续时间
    t = np.linspace(0, duration, int(sample_rate * duration), endpoint=False)
    
    # 50Hz和120Hz的正弦波
    signal = np.sin(2 * np.pi * 50 * t) + 0.5 * np.sin(2 * np.pi * 120 * t)
    
    # 添加噪声
    signal += 0.2 * np.random.randn(len(t))
    
    # 执行FFT
    fft_result = fft(signal[:1024])  # 取前1024个点
    
    # 计算频率轴
    freqs = np.fft.fftfreq(1024, 1/sample_rate)
    
    # 绘制结果
    import matplotlib.pyplot as plt
    
    plt.figure(figsize=(12, 6))
    
    plt.subplot(2, 1, 1)
    plt.plot(t[:1024], signal[:1024])
    plt.title('时域信号')
    plt.xlabel('时间 (s)')
    plt.ylabel('幅度')
    
    plt.subplot(2, 1, 2)
    plt.plot(freqs[:512], np.abs(fft_result)[:512])
    plt.title('频域表示')
    plt.xlabel('频率 (Hz)')
    plt.ylabel('幅度')
    
    plt.tight_layout()
    plt.show()
 
if __name__ == "__main__":
    test_fft()

def fft_iterative(x):
    """
    迭代实现的FFT算法（Cooley-Tukey）
    """
    n = len(x)
    
    # 位反转重排
    def bit_reverse_copy(a):
        result = a.copy()
        j = 0
        for i in range(1, n):
            bit = n >> 1
            while j & bit:
                j ^= bit
                bit >>= 1
            j ^= bit
            if i < j:
                result[i], result[j] = result[j], result[i]
        return result
    
    # 位反转重排
    x = bit_reverse_copy(x)
    
    # 蝶形运算
    size = 2
    while size <= n:
        half_size = size // 2
        twiddle_step = cmath.exp(-2j * cmath.pi / size)
        
        for i in range(0, n, size):
            twiddle = 1
            for j in range(half_size):
                # 蝶形运算
                temp = twiddle * x[i + j + half_size]
                x[i + j + half_size] = x[i + j] - temp
                x[i + j] = x[i + j] + temp
                twiddle *= twiddle_step
        
        size *= 2
    
    return x

class FFT {
    constructor() {
        this.PI2 = Math.PI * 2;
    }
    
    // 复数乘法
    complexMultiply(a, b) {
        return {
            real: a.real * b.real - a.imag * b.imag,
            imag: a.real * b.imag + a.imag * b.real
        };
    }
    
    // 复数加法
    complexAdd(a, b) {
        return {
            real: a.real + b.real,
            imag: a.imag + b.imag
        };
    }
    
    // 复数减法
    complexSubtract(a, b) {
        return {
            real: a.real - b.real,
            imag: a.imag - b.imag
        };
    }
    
    // 快速傅里叶变换
    fft(signal) {
        const n = signal.length;
        
        if (n <= 1) return signal;
        
        if (n % 2 !== 0) {
            throw new Error("信号长度必须是2的幂次");
        }
        
        // 分离偶数和奇数索引
        const even = [];
        const odd = [];
        
        for (let i = 0; i < n; i += 2) {
            even.push(signal[i]);
            odd.push(signal[i + 1]);
        }
        
        // 递归计算
        const evenFFT = this.fft(even);
        const oddFFT = this.fft(odd);
        
        // 组合结果
        const result = new Array(n);
        
        for (let k = 0; k < n / 2; k++) {
            // 计算旋转因子
            const angle = -this.PI2 * k / n;
            const twiddle = {
                real: Math.cos(angle),
                imag: Math.sin(angle)
            };
            
            const twiddleOdd = this.complexMultiply(twiddle, oddFFT[k]);
            
            result[k] = this.complexAdd(evenFFT[k], twiddleOdd);
            result[k + n / 2] = this.complexSubtract(evenFFT[k], twiddleOdd);
        }
        
        return result;
    }
    
    // 计算幅度谱
    magnitudeSpectrum(fftResult) {
        return fftResult.map(complex => 
            Math.sqrt(complex.real * complex.real + complex.imag * complex.imag)
        );
    }
    
    // 计算功率谱
    powerSpectrum(fftResult) {
        return this.magnitudeSpectrum(fftResult).map(mag => mag * mag);
    }
}
 
// Web Audio API集成
class AudioFFTProcessor {
    constructor() {
        this.fft = new FFT();
        this.audioContext = null;
        this.analyser = null;
        this.dataArray = null;
    }
    
    async initialize() {
        this.audioContext = new (window.AudioContext || window.webkitAudioContext)();
        this.analyser = this.audioContext.createAnalyser();
        this.analyser.fftSize = 2048;
        
        const bufferLength = this.analyser.frequencyBinCount;
        this.dataArray = new Uint8Array(bufferLength);
        
        // 获取用户媒体
        const stream = await navigator.mediaDevices.getUserMedia({ audio: true });
        const source = this.audioContext.createMediaStreamSource(stream);
        source.connect(this.analyser);
    }
    
    processAudio() {
        this.analyser.getByteTimeDomainData(this.dataArray);
        
        // 转换为复数格式
        const complexSignal = Array.from(this.dataArray).map(value => ({
            real: (value - 128) / 128, // 归一化
            imag: 0
        }));
        
        // 执行FFT
        const fftResult = this.fft.fft(complexSignal);
        
        // 获取幅度谱
        const magnitude = this.fft.magnitudeSpectrum(fftResult);
        
        return {
            timeDomain: this.dataArray,
            frequencyDomain: magnitude
        };
    }
    
    visualize(canvas) {
        const ctx = canvas.getContext('2d');
        const width = canvas.width;
        const height = canvas.height;
        
        const draw = () => {
            const data = this.processAudio();
            
            ctx.fillStyle = 'rgb(0, 0, 0)';
            ctx.fillRect(0, 0, width, height);
            
            // 绘制频谱
            const barWidth = width / data.frequencyDomain.length;
            
            for (let i = 0; i < data.frequencyDomain.length; i++) {
                const barHeight = (data.frequencyDomain[i] / 255) * height;
                
                const r = barHeight + 25 * (i / data.frequencyDomain.length);
                const g = 250 * (i / data.frequencyDomain.length);
                const b = 50;
                
                ctx.fillStyle = `rgb(${r},${g},${b})`;
                ctx.fillRect(i * barWidth, height - barHeight, barWidth, barHeight);
            }
            
            requestAnimationFrame(draw);
        };
        
        draw();
    }
}
 
// 使用示例
const audioProcessor = new AudioFFTProcessor();
 
audioProcessor.initialize().then(() => {
    const canvas = document.getElementById('fftCanvas');
    audioProcessor.visualize(canvas);
}).catch(error => {
    console.error('初始化失败:', error);
});

class OptimizedFFT:
    def __init__(self, max_size):
        self.max_size = max_size
        self.precomputed_twiddles = {}
        self._precompute_twiddles()
    
    def _precompute_twiddles(self):
        """预计算所有可能的旋转因子"""
        for size in [2**i for i in range(1, int(np.log2(self.max_size)) + 1)]:
            twiddles = []
            for k in range(size // 2):
                angle = -2j * np.pi * k / size
                twiddles.append(np.exp(angle))
            self.precomputed_twiddles[size] = twiddles
    
    def fft(self, x):
        """使用预计算旋转因子的FFT"""
        n = len(x)
        if n <= 1:
            return x
        
        # 使用预计算的旋转因子
        even = self.fft(x[0::2])
        odd = self.fft(x[1::2])
        
        result = [0] * n
        twiddles = self.precomputed_twiddles.get(n, [])
        
        for k in range(n // 2):
            twiddle = twiddles[k] if k < len(twiddles) else np.exp(-2j * np.pi * k / n)
            result[k] = even[k] + twiddle * odd[k]
            result[k + n // 2] = even[k] - twiddle * odd[k]
        
        return result

def fft_in_place(x):
    """原地FFT，减少内存分配"""
    n = len(x)
    
    # 位反转重排（原地）
    def bit_reverse_swap(arr):
        j = 0
        for i in range(1, n):
            bit = n >> 1
            while j & bit:
                j ^= bit
                bit >>= 1
            j ^= bit
            if i < j:
                arr[i], arr[j] = arr[j], arr[i]
    
    bit_reverse_swap(x)
    
    # 蝶形运算（原地）
    size = 2
    while size <= n:
        half_size = size // 2
        twiddle_step = np.exp(-2j * np.pi / size)
        
        for i in range(0, n, size):
            twiddle = 1
            for j in range(half_size):
                # 原地蝶形运算
                temp = twiddle * x[i + j + half_size]
                x[i + j + half_size] = x[i + j] - temp
                x[i + j] = x[i + j] + temp
                twiddle *= twiddle_step
        
        size *= 2
    
    return x

import concurrent.futures
import numpy as np
 
def parallel_fft(x, num_threads=4):
    """并行FFT实现"""
    n = len(x)
    
    if n <= 1024:  # 小数组使用串行算法
        return fft(x)
    
    # 分治并行处理
    with concurrent.futures.ThreadPoolExecutor(max_workers=num_threads) as executor:
        # 并行计算偶数和奇数部分
        even_future = executor.submit(fft, x[0::2])
        odd_future = executor.submit(fft, x[1::2])
        
        even = even_future.result()
        odd = odd_future.result()
    
    # 组合结果
    result = [0] * n
    for k in range(n // 2):
        twiddle = np.exp(-2j * np.pi * k / n)
        result[k] = even[k] + twiddle * odd[k]
        result[k + n // 2] = even[k] - twiddle * odd[k]
    
    return result

class AudioProcessor:
    def __init__(self, sample_rate=44100):
        self.sample_rate = sample_rate
        self.fft_size = 2048
        self.hop_size = 512
    
    def analyze_pitch(self, audio_data):
        """音高检测"""
        # 分帧处理
        frames = self._frame_audio(audio_data)
        
        pitches = []
        for frame in frames:
            # 加窗处理
            windowed = frame * np.hanning(len(frame))
            
            # FFT变换
            spectrum = fft(windowed)
            
            # 寻找峰值频率
            magnitudes = np.abs(spectrum[:len(spectrum)//2])
            peak_idx = np.argmax(magnitudes)
            
            # 转换为频率
            frequency = peak_idx * self.sample_rate / len(frame)
            pitches.append(frequency)
        
        return pitches
    
    def remove_noise(self, audio_data, noise_threshold=0.1):
        """降噪处理"""
        # 转换为频域
        spectrum = fft(audio_data)
        
        # 阈值处理
        magnitude = np.abs(spectrum)
        phase = np.angle(spectrum)
        
        # 软阈值降噪
        magnitude_clean = np.maximum(0, magnitude - noise_threshold)
        
        # 重构信号
        clean_spectrum = magnitude_clean * np.exp(1j * phase)
        clean_audio = np.real(ifft(clean_spectrum))
        
        return clean_audio
    
    def _frame_audio(self, audio_data):
        """音频分帧"""
        frames = []
        for i in range(0, len(audio_data) - self.fft_size + 1, self.hop_size):
            frames.append(audio_data[i:i + self.fft_size])
        return frames

class ImageFFTProcessor:
    def __init__(self):
        pass
    
    def apply_low_pass_filter(self, image, cutoff_frequency=0.1):
        """低通滤波器"""
        # 2D FFT
        fft_result = np.fft.fft2(image)
        fft_shifted = np.fft.fftshift(fft_result)
        
        # 创建滤波器
        rows, cols = image.shape
        crow, ccol = rows // 2, cols // 2
        
        mask = np.zeros((rows, cols), np.uint8)
        radius = int(min(rows, cols) * cutoff_frequency / 2)
        cv2.circle(mask, (ccol, crow), radius, 1, -1)
        
        # 应用滤波器
        filtered_fft = fft_shifted * mask
        
        # 逆变换
        filtered_image = np.fft.ifft2(np.fft.ifftshift(filtered_fft))
        return np.real(filtered_image)
    
    def detect_edges(self, image):
        """边缘检测"""
        # 2D FFT
        fft_result = np.fft.fft2(image)
        
        # 高通滤波器
        rows, cols = image.shape
        center_row, center_col = rows // 2, cols // 2
        
        # 创建高通滤波器
        mask = np.ones((rows, cols))
        mask[center_row-10:center_row+10, center_col-10:center_col+10] = 0
        
        # 应用滤波器
        filtered_fft = fft_result * mask
        
        # 逆变换
        edge_image = np.fft.ifft2(filtered_fft)
        return np.abs(edge_image)
    
    def compress_image(self, image, compression_ratio=0.8):
        """图像压缩"""
        # 2D FFT
        fft_result = np.fft.fft2(image)
        
        # 获取幅度和相位
        magnitude = np.abs(fft_result)
        phase = np.angle(fft_result)
        
        # 保留主要频率成分
        total_coeffs = magnitude.size
        retain_coeffs = int(total_coeffs * (1 - compression_ratio))
        
        # 找到最重要的系数
        flat_magnitude = magnitude.flatten()
        threshold_idx = np.argsort(flat_magnitude)[-retain_coeffs]
        threshold = flat_magnitude[threshold_idx]
        
        # 阈值处理
        mask = magnitude >= threshold
        compressed_fft = fft_result * mask
        
        # 逆变换
        compressed_image = np.fft.ifft2(compressed_fft)
        return np.real(compressed_image)

class OFDMModulator:
    def __init__(self, num_subcarriers=64, cp_length=16):
        self.num_subcarriers = num_subcarriers
        self.cp_length = cp_length
        self.subcarrier_spacing = 1.0  # 子载波间隔
    
    def modulate(self, data_symbols):
        """OFDM调制"""
        # QAM调制（简化版）
        modulated = self._qam_modulate(data_symbols)
        
        # IFFT转换为时域信号
        time_domain = np.fft.ifft(modulated, self.num_subcarriers)
        
        # 添加循环前缀
        cp = time_domain[-self.cp_length:]
        ofdm_symbol = np.concatenate([cp, time_domain])
        
        return ofdm_symbol
    
    def demodulate(self, received_signal):
        """OFDM解调"""
        # 移除循环前缀
        useful_signal = received_signal[self.cp_length:]
        
        # FFT转换为频域
        frequency_domain = np.fft.fft(useful_signal)
        
        # QAM解调
        demodulated = self._qam_demodulate(frequency_domain)
        
        return demodulated
    
    def _qam_modulate(self, data_bits):
        """QAM调制（16-QAM简化版）"""
        # 将比特映射到复数符号
        symbols = []
        for i in range(0, len(data_bits), 4):
            if i + 3 < len(data_bits):
                # 16-QAM映射
                real_part = (data_bits[i] * 2 + data_bits[i+1]) - 1.5
                imag_part = (data_bits[i+2] * 2 + data_bits[i+3]) - 1.5
                symbols.append(complex(real_part, imag_part))
        
        # 填充到子载波数量
        while len(symbols) < self.num_subcarriers:
            symbols.append(0)
        
        return np.array(symbols[:self.num_subcarriers])
    
    def _qam_demodulate(self, received_symbols):
        """QAM解调"""
        demodulated_bits = []
        
        for symbol in received_symbols:
            # 16-QAM解映射
            real_part = 1 if symbol.real > 0 else 0
            imag_part = 1 if symbol.imag > 0 else 0
            
            demodulated_bits.extend([real_part, imag_part])
        
        return demodulated_bits

def debug_fft(x, debug_level=1):
    """带调试信息的FFT实现"""
    n = len(x)
    
    if debug_level >= 1:
        print(f"FFT输入: 长度={n}, 数据={x[:5]}...")  # TRAE断点位置
    
    if n <= 1:
        return x
    
    # 分治处理
    even = debug_fft(x[0::2], debug_level + 1)
    odd = debug_fft(x[1::2], debug_level + 1)
    
    if debug_level >= 2:
        print(f"第{debug_level}层 - 偶数部分结果: {even[:3]}...")
        print(f"第{debug_level}层 - 奇数部分结果: {odd[:3]}...")
    
    result = [0] * n
    for k in range(n // 2):
        twiddle = cmath.exp(-2j * cmath.pi * k / n)
        result[k] = even[k] + twiddle * odd[k]
        result[k + n // 2] = even[k] - twiddle * odd[k]
        
        if debug_level >= 3 and k < 2:
            print(f"k={k}: twiddle={twiddle}, result[{k}]={result[k]}")
    
    return result

import time
import memory_profiler
 
@memory_profiler.profile
def benchmark_fft implementations():
    """性能基准测试"""
    sizes = [128, 256, 512, 1024, 2048, 4096]
    
    results = {}
    
    for size in sizes:
        # 生成测试数据
        test_data = np.random.random(size) + 1j * np.random.random(size)
        
        # 测试自定义FFT
        start_time = time.time()
        custom_result = fft(test_data)
        custom_time = time.time() - start_time
        
        # 测试NumPy FFT
        start_time = time.time()
        numpy_result = np.fft.fft(test_data)
        numpy_time = time.time() - start_time
        
        # 计算误差
        error = np.mean(np.abs(custom_result - numpy_result))
        
        results[size] = {
            'custom_time': custom_time,
            'numpy_time': numpy_time,
            'speedup': numpy_time / custom_time,
            'error': error
        }
        
        print(f"数组大小 {size}:")
        print(f"  自定义FFT: {custom_time:.6f}s")
        print(f"  NumPy FFT: {numpy_time:.6f}s")
        print(f"  加速比: {numpy_time/custom_time:.2f}x")
        print(f"  平均误差: {error:.2e}")
        print()
    
    return results
 
# 运行基准测试
if __name__ == "__main__":
    benchmark_results = benchmark_fft_implementations()