验证FFT算法的正确性

3
今天我编写了一个算法,用于计算给定表示离散函数的点阵的快速傅里叶变换。现在我正在尝试测试它是否正常工作。我尝试了大约十几组不同的输入集,并且它们似乎与我在网上找到的示例相匹配。但是,在我的最后一个测试中,我将cos(i/2)作为输入,其中i从0到31,根据使用的求解器,我得到了3个不同的结果。我的解决方案似乎是最不准确的:
以下是我的代码,以防有帮助:
/**
 * Slices the original array, starting with start, grabbing every stride elements.
 * For example, slice(A, 3, 4, 5) would return elements 3, 8, 13, and 18 from array A.
 * @param array     The array to be sliced
 * @param start     The starting index
 * @param newLength The length of the final array
 * @param stride    The spacing between elements to be selected
 * @return          A sliced copy of the input array
 */
public double[] slice(double[] array, int start, int newLength, int stride) {
    double[] newArray = new double[newLength];
    int count = 0;
    for (int i = start; count < newLength && i < array.length; i += stride) {
        newArray[count++] = array[i];
    }
    return newArray;
}

/**
 * Calculates the fast fourier transform of the given function.  The parameters are updated with the calculated values
 * To ignore all imaginary output, leave imaginary null
 * @param real An array representing the real part of a discrete-time function
 * @param imaginary An array representing the imaginary part of a discrete-time function
 * Pre: If imaginary is not null, the two arrays must be the same length, which must be a power of 2
 */
public void fft(double[] real, double[] imaginary) throws IllegalArgumentException {
    if (real == null) {
        throw new NullPointerException("Real array cannot be null");
    }

    int N = real.length;

    // Make sure the length is a power of 2
    if ((Math.log(N) / Math.log(2)) % 1 != 0) {
        throw new IllegalArgumentException("The array length must be a power of 2");
    }

    if (imaginary != null && imaginary.length != N) {
        throw new IllegalArgumentException("The two arrays must be the same length");
    }

    if (N == 1) {
        return;
    }

    double[] even_re = slice(real, 0, N/2, 2);
    double[] odd_re = slice(real, 1, N/2, 2);

    double[] even_im = null;
    double[] odd_im = null;
    if (imaginary != null) {
        even_im = slice(imaginary, 0, N/2, 2);
        odd_im = slice(imaginary, 1, N/2, 2);
    }

    fft(even_re, even_im);
    fft(odd_re, odd_im);

    // F[k] = real[k] + imaginary[k]

    //              even   odd
    //       F[k] = E[k] + O[k] * e^(-i*2*pi*k/N)
    // F[k + N/2] = E[k] - O[k] * e^(-i*2*pi*k/N)

    // Split complex arrays into component arrays:
    // E[k] = er[k] + i*ei[k]
    // O[k] = or[k] + i*oi[k]

    // e^ix = cos(x) + i*sin(x)

    // Let x = -2*pi*k/N
    // F[k] = er[k] + i*ei[k] + (or[k] + i*oi[k])(cos(x) + i*sin(x))
    //      = er[k] + i*ei[k] + or[k]cos(x) + i*or[k]sin(x) + i*oi[k]cos(x) - oi[k]sin(x)
    //      = (er[k] + or[k]cos(x) - oi[k]sin(x)) + i*(ei[k] + or[k]sin(x) + oi[k]cos(x))
    //        {               real              }     {            imaginary            }

    // F[k + N/2] = (er[k] - or[k]cos(x) + oi[k]sin(x)) + i*(ei[k] - or[k]sin(x) - oi[k]cos(x))
    //              {               real              }     {            imaginary            }

    // Ignoring all imaginary parts (oi = 0):
    //       F[k] = er[k] + or[k]cos(x)
    // F[k + N/2] = er[k] - or[k]cos(x)

    for (int k = 0; k < N/2; ++k) {
        double t = odd_re[k] * Math.cos(-2 * Math.PI * k/N);
        real[k]       = even_re[k] + t;
        real[k + N/2] = even_re[k] - t;

        if (imaginary != null) {
            t = odd_im[k] * Math.sin(-2 * Math.PI * k/N);
            real[k]       -= t;
            real[k + N/2] += t;

            double t1 = odd_re[k] * Math.sin(-2 * Math.PI * k/N);
            double t2 = odd_im[k] * Math.cos(-2 * Math.PI * k/N);
            imaginary[k]       = even_im[k] + t1 + t2;
            imaginary[k + N/2] = even_im[k] - t1 - t2;
        }
    }
}

3
第二张图明显有误;实际输入应该产生对称的输出。 - Oliver Charlesworth
4
如果你想验证你的实现(忽略浮点数的数值问题),只需编写一个简单的DFT(大约五行代码)并进行比较即可。 - Oliver Charlesworth
2个回答

1
  1. 验证

    看这里:slow DFT,iDFT 最后是我慢速实现的 DFTiDFT,它们经过测试并且正确。我过去也用它们来进行快速实现的验证。

  2. 你的代码

    停止递归的方式是错误的(你忘记设置返回元素了),我的代码如下:

    if (n<=1) { if (n==1) { dst[0]=src[0]*2.0; dst[1]=src[1]*2.0; } return; }
    

    所以当你的 N==1 时,在 return 前将输出元素设置为 Re=2.0*real[0], Im=2.0*imaginary[0]。此外,我在你的复杂数学中有些迷失,并且懒得分析。

为了确保,这是我的快速实现。它需要从类层次结构中获取太多的东西,因此除了与您的代码进行视觉比较之外,对您没有其他用途。

我的快速实现(cc表示复杂的输出和输入):

//---------------------------------------------------------------------------
void transform::DFFTcc(double *dst,double *src,int n)
    {
    if (n>N) init(n);
    if (n<=1) { if (n==1) { dst[0]=src[0]*2.0; dst[1]=src[1]*2.0; } return; }
    int i,j,n2=n>>1,q,dq=+N/n,mq=N-1;
    // reorder even,odd (buterfly)
    for (j=0,i=0;i<n+n;) { dst[j]=src[i]; i++; j++; dst[j]=src[i]; i+=3; j++; }
    for (    i=2;i<n+n;) { dst[j]=src[i]; i++; j++; dst[j]=src[i]; i+=3; j++; }
    // recursion
    DFFTcc(src  ,dst  ,n2); // even
    DFFTcc(src+n,dst+n,n2); // odd
    // reorder and weight back (buterfly)
    double a0,a1,b0,b1,a,b;
    for (q=0,i=0,j=n;i<n;i+=2,j+=2,q=(q+dq)&mq)
        {
        a0=src[j  ]; a1=+_cos[q];
        b0=src[j+1]; b1=+_sin[q];
        a=(a0*a1)-(b0*b1);
        b=(a0*b1)+(a1*b0);
        a0=src[i  ]; a1=a;
        b0=src[i+1]; b1=b;
        dst[i  ]=(a0+a1)*0.5;
        dst[i+1]=(b0+b1)*0.5;
        dst[j  ]=(a0-a1)*0.5;
        dst[j+1]=(b0-b1)*0.5;
        }
    }
//---------------------------------------------------------------------------


dst[]src[]不重叠!!! 因此,您不能将数组转换为自身。
_cos_sin是预计算的cossin值的表格(由init()函数计算如下:

    double a,da; int i;
    da=2.0*M_PI/double(N);
    for (a=0.0,i=0;i<N;i++,a+=da) { _cos[i]=cos(a); _sin[i]=sin(a); }


N2的幂次方(数据集的零填充大小)(从init(n)调用中的最后一个n

为了完整起见,这是我复杂到复杂的慢版本:

//---------------------------------------------------------------------------
void transform::DFTcc(double *dst,double *src,int n)
    {
    int i,j;
    double a,b,a0,a1,_n,b0,b1,q,qq,dq;
    dq=+2.0*M_PI/double(n); _n=2.0/double(n);
    for (q=0.0,j=0;j<n;j++,q+=dq)
        {
        a=0.0; b=0.0;
        for (qq=0.0,i=0;i<n;i++,qq+=q)
            {
            a0=src[i+i  ]; a1=+cos(qq);
            b0=src[i+i+1]; b1=+sin(qq);
            a+=(a0*a1)-(b0*b1);
            b+=(a0*b1)+(a1*b0);
            }
        dst[j+j  ]=a*_n;
        dst[j+j+1]=b*_n;
        }
    }
//---------------------------------------------------------------------------

DFFTcc 中的 D 是否表示离散?我问这个问题是因为我认为 FFT 已经暗示了它的离散性质。 - Eduardo Reis
1
@EduardoReis 是的,D代表离散,而FFT并不意味着离散,因为它也可以通过模拟计算机/电路来实现... - Spektre

0

我会使用像 Wolfram Alpha 这样权威的工具进行验证。

如果我对 0 <= i < 32 进行 cos(i/2) 的计算,我会得到以下数组:

[1,0.878,0.540,0.071,-0.416,-0.801,-0.990,-0.936,-0.654,-0.211,0.284,0.709,0.960,0.977,0.754,0.347,-0.146,-0.602,-0.911,-0.997,-0.839,-0.476,0.004,0.483,0.844,0.998,0.907,0.595,0.137,-0.355,-0.760,-0.978]

如果我将其作为输入提供给 Wolfram Alpha 的 FFT 函数,我会得到 这个结果
我得到的图形看起来是对称的,这很有意义。这个图形看起来与你提供的任何一个都不一样。

网页内容由stack overflow 提供, 点击上面的
可以查看英文原文,
原文链接