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; }
[FFTW] Cyclic results with FFTW?
Started by ●March 17, 2005
Reply by ●March 17, 20052005-03-17
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... > > FredFred: 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
Reply by ●March 17, 20052005-03-17
(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.
Reply by ●March 17, 20052005-03-17
(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.)
Reply by ●March 17, 20052005-03-17
"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