DSPRelated.com
Forums

[FFTW] Cyclic results with FFTW?

Started by Fred March 17, 2005
Dear all, 

I implemented some test FFT's from various libraries, having
decomposed and reconstructed a simple signal (sine), my code worked
for all of them however, the signal I got using fftw was not
reconstructed correctly, indeed, the transform is not cyclic.

How to get cyclic output from fftw? 

Below is a code in 1D I used. sorry to ask such a trivial question and
thanks for your answer...

Fred



--------- C CODE ----------

#include <stdio.h>
#include <math.h>

#include <fftw3.h>

unsigned flags = FFTW_ESTIMATE  | FFTW_PRESERVE_INPUT; 

/*    volontary not equal to 2^n (n E ZZ+*)   */
    F = fopen("FULLT", "w");


    for(i=0; i<ARRAY_SIZE; i++) {
        A[i] = sin(i/10.);
        B[i] = 0;
    }

    for(i=0; i<ARRAY_SIZE; i++) {
        fprintf(o, "%u\t%f\n", i, A[i]);
    }

    /* this just interleaves A & B into C */
    interleave(C, A, B, ARRAY_SIZE, 0);

    for(i=0; i<ARRAY_SIZE; i++) {
        fprintf(R, "%u\t%f\n", i, A[i]);
    }

    p = fftw_plan_dft_r2c_1d(ARRAY_SIZE, A, C, flags);
    fftw_execute(p);

    /* this just splits C into A & B */
    split(A, B, (const fftw_complex * const) C, ARRAY_SIZE, 0);

    for(i=0; i<ARRAY_SIZE; i++) {
        fprintf(R, "%u\t%f\n", i, A[i]);
        fprintf(I, "%u\t%f\n", i, B[i]);
        fprintf(F, "%u\t%f\n", i, sqrt(A[i]*A[i] + B[i]*B[i]));
    }

    interleave(A, B, ARRAY_SIZE, 0);

    invp = fftw_plan_dft_c2r_1d(ARRAY_SIZE, C, A, flags);
    fftw_execute(invp);

    split B, (const fftw_complex * const) C, ARRAY_SIZE, 0);

    for(i=0; i<ARRAY_SIZE; i++) {
        fprintf(O, "%u\t%f\n", i, A[i]/ARRAY_SIZE);
    }

    fftw_destroy_plan(p);
    fftw_destroy_plan(invp);
    fftw_cleanup();

    fclose(I);
    fclose(R);
    fclose(O);
    fclose(o);
    fclose(F);

    return EXIT_SUCCESS;
}
Fred wrote:
> Dear all, > > I implemented some test FFT's from various libraries, having > decomposed and reconstructed a simple signal (sine), my code worked > for all of them however, the signal I got using fftw was not > reconstructed correctly, indeed, the transform is not cyclic. > > How to get cyclic output from fftw? > > Below is a code in 1D I used. sorry to ask such a trivial question and > thanks for your answer... > > Fred
Fred: The FFTW output is not in [-f,0,+f] order. FFTW stores the data starting at DC, then positive frequencies, and lastly the negative frequencies beginning with the most negative frequency. Good luck, OUP
(What do you mean by "not cyclic"?)

Rest assured that FFTW computes the DFT correctly, so any error is
almost certainly in your own code.

Most likely you are simply misunderstanding the output format of FFTW's
r2c transforms (which are DFTs specialized for real inputs). For N
inputs, the output is a complex array, but it is *not* of the same
length N.  Rather, it is of length N/2+1.  The reason is that the
outputs of a real-input DFT are redundant by a conjugate symmetry, so
only ~half of them are outputted.

This is all described in detail by the FFTW manual.

Cordially,
Steven G. Johnson

PS. Or it could be that you have bugs in your "split" routines, etc.
These sorts of transformations aren't needed for FFTW.  e.g. it's not
clear to me why you would call split after the c2r transform, whose
output is already a real array. It's impossible to evaluate, of course,
because you don't post working code...not that you should expect others
to debug your code, anyway.

(I'm not aware of *any* widely used FFT code that outputs in [-f, 0,
+f] order, so I doubt that is the poster's problem.)

"Fred" <protected.ggroups.5.nomades@spamgourmet.com> wrote in message 
news:aaa048d2.0503170512.116504d2@posting.google.com...
> Dear all, > > I implemented some test FFT's from various libraries, having > decomposed and reconstructed a simple signal (sine), my code worked > for all of them however, the signal I got using fftw was not > reconstructed correctly, indeed, the transform is not cyclic. > > How to get cyclic output from fftw?
Fred, First, do you mean cyclic in time or in frequency? You won't get "cyclic" output from fftw or any other fft (i.e. forward transform time > frequency) unless you do some particular things: By definition, a normal fft addresses only a single "cycle" of the spectrum. You can repeat it if you want to make it cyclic. That's equivalent to increasing the sample rate - zero stuffing between the time samples - without adding non-zero sample values (as would be the case in interpolation - which is a bit different or more involved). If the temporal record is cyclic then a "round trip" FFT / IFFT should result in the same sequence as was input -neglecting any numerical roundoff type noise that's added. Fred