Verifying correctness of FFT algorithm

Question

Today I wrote an algorithm to compute the Fast Fourier Transform from a given array of points representing a discrete function. Now I'm trying to test it to see if it is working. I've tried about a dozen different input sets, and they seem to match up with examples I've found online. However, for my final test, I gave it the input of cos(i / 2), with i from 0 to 31, and I've gotten 3 different results based on which solver I use. My solution seems to be the least accurate:

Does this indicate a problem with my algorithm, or is it simply a result of the relatively small data set?

My code is below, in case it helps:

/**
 * 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;
        }
    }
}

Basanite

The second graph is clearly wrong; a real input should result in a symmetric output. – Paleogeography

If you want to verify your implementation (floating-point numerical issues notwithstanding), just write a naive DFT (about five lines of code), and compare. – Paleogeography

Answer 1

1. Validation

   look here: slow DFT,iDFT at the end is mine slow implementation of DFT and iDFT they are tested and correct. I also used them for fast implementations validation in the past.

2. Your code

   stop recursion is wrong (you forget to set the return element) mine looks like this:

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

   so when your N==1 set the output element to Re=2.0*real[0], Im=2.0*imaginary[0] before return. Also I am a bit lost in your complex math (t,t1,t2) and to lazy to analyze.

Just to be sure here is mine fast implementation. It need too much things from class hierarchy so it will not be of another use for you then visual comparison to your code.

My Fast implementation (cc means complex output and input):

//---------------------------------------------------------------------------
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[] and src[] are not overlapping !!! so you cannot transform array to itself .
_cos and _sin are precomputed tables of cos and sin values (computed by init() function like this:

    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); }

N is power of 2 (zero padded size of data set) (last n from init(n) call)

Just to be complete here is mine complex to complex slow version:

//---------------------------------------------------------------------------
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;
        }
    }
//---------------------------------------------------------------------------

Philbin

Is the D in DFFTcc meaning discrete? Asking this since I thought FFT would already imply the discrete aspect of it. – Thrust

@EduardoReis yes D means discrete and no FFT does not imply discrete as it might be implemented by analog computer/circuit too... – Philbin

Answer 2

I'd use something authoritative like Wolfram Alpha to verify.

If I evalute cos(i/2) for 0 <= i < 32, I get this array:

[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]

If I give that as input to Wolfram Alpha's FFT function I get this result.

The plot that I get looks symmetric, which makes sense. The plot looks nothing like any of the ones you supplied.

Sherilyn