为何FFTW3的DST Type-1函数运行缓慢？如何用FFT实现快速DST？

核心思路 #

DST Type-1（对应FFTW的FFTW_RODFT00）可通过构造反对称扩展序列 + 复数FFT实现，借助FFT的O(n logn)时间复杂度提升运算速度，避开原生RODFT00的低效问题。

具体实现步骤 #

假设输入序列input[0..N-1]长度为N，要得到对应的DST Type-1变换结果output[0..N-1]，步骤如下：

1. 构造反对称扩展序列
   创建长度为M = 2*(N+1)的复数序列y，满足反对称规则：

   • y[0]的实部、虚部均为0
   • 对1 ≤ n ≤ N，y[n]的实部为input[n-1]，虚部为0
   • y[N+1]的实部、虚部均为0
   • 对1 ≤ n ≤ N，y[M - n]的实部为-y[n]的实部，虚部为-y[n]的虚部

2. 执行复数FFT
   用FFTW的复数FFT计算序列y的傅里叶变换，得到复数序列Y。

3. 提取结果并匹配输出
   DST Type-1的结果对应Y的虚部，直接提取即可匹配FFTW_RODFT00的无归一化输出：

   • 对0 ≤ k ≤ N-1，output[k] = Y[k+1]的虚部

FFTW代码示例 #

#include <fftw3.h>

void fast_dst_type1(int N, double *input, double *output) {
    int M = 2 * (N + 1);
    fftw_complex *y = fftw_malloc(sizeof(fftw_complex) * M);
    fftw_complex *Y = fftw_malloc(sizeof(fftw_complex) * M);
    
    // 初始化反对称扩展序列
    y[0][0] = 0.0; y[0][1] = 0.0;
    for (int n = 1; n <= N; ++n) {
        y[n][0] = input[n-1];
        y[n][1] = 0.0;
    }
    y[N+1][0] = 0.0; y[N+1][1] = 0.0;
    for (int n = 1; n <= N; ++n) {
        y[M - n][0] = -y[n][0];
        y[M - n][1] = -y[n][1];
    }
    
    // 创建FFT计划并执行
    fftw_plan plan = fftw_plan_dft_1d(M, y, Y, FFTW_FORWARD, FFTW_ESTIMATE);
    fftw_execute(plan);
    
    // 提取DST结果
    for (int k = 0; k < N; ++k) {
        output[k] = Y[k+1][1];
    }
    
    // 清理资源
    fftw_destroy_plan(plan);
    fftw_free(y);
    fftw_free(Y);
}

注意事项 #

• 扩展序列长度M=2(N+1)是为了保证反对称性完整，让FFT能精准提取正弦分量。
• 若需要正交归一化的DST结果，可在最后对output序列除以2*(N+1)。
• 该方法复用了FFTW高度优化的复数FFT实现，速度远优于O(n²)的手动实现和原生RODFT00。

内容的提问来源于stack exchange，提问作者Hansov