Forums

fftw: real/even inverse fft

Started by tania July 4, 2005
I wrote the following program to compute the inverse fft of a real/even
list of data in the frequency domain.
I expected to obtain a real/even list of data in the time domain but
this is not the case.
Please, what am I doing wrong??

#include <stdio.h>
#include <string.h>
#include <fftw3.h>

#define N 512

main()
{

 int i;
 float t,f, xt, yt, xf, yf;
 fftw_complex *time, *freq;
 fftw_plan p;

 FILE *pftime, *pffreq;

 time = fftw_malloc(sizeof(fftw_complex)*N);
 freq = fftw_malloc(sizeof(fftw_complex)*N);
 p=fftw_plan_dft_1d(N, freq, time, FFTW_BACKWARD, FFTW_ESTIMATE);
 pftime = fopen ("result.dat","w");
 pffreq = fopen ("data.dat","r");

 // read the frequency 0
 fscanf(pffreq,"%f%f\n",&f,&xf);
 freq[0][0]=xf;
 freq[0][1]=0.;

 // read the positive frequencies

 for (i = 1; i < N/2 +1; i ++)
  {
   fscanf(pffreq,"%f%f\n",&f,&xf);
   freq[i][0]=freq[N-i][0]=xf.;
   freq[i][1]=freq[N-i][1]=0;
  }


 fftw_execute(p);
 for (i = 0; i < N; i ++)
  {
   xt=time[i][0];
   yt=time[i][1];
   fprintf(pftime,"%e %e\n",xt/N,yt/N);
  }

 fftw_destroy_plan(p);
 fftw_free(time);fftw_free(freq);
 fclose(pffreq);fclose(pftime);
 
}

"tania" <tania.re@wanadoo.fr> wrote in message 
news:1120514787.979440.44980@z14g2000cwz.googlegroups.com...
>I wrote the following program to compute the inverse fft of a real/even > list of data in the frequency domain. > I expected to obtain a real/even list of data in the time domain but > this is not the case. > Please, what am I doing wrong??
I almost never read code. Here is what could be going wrong: "Even" means "symmetrical around zero" Generally when we do an fft or ifft, the zero time or zero frequency sample is the first sample. An example of an even sequence of N=10 elements would be like this: A) [1 2 3 0 0 0 0 0 3 2] If you rotate this sequence it would look like this: B) [0 0 0 3 2 1 2 3 0 0] which looks "even" but entails a delay or a frequency shift - depending on which domain these elements belong to. The assumption is that the 1st term repeats as the N+1th=11th sample and is not part of the sequence. Unless you start with an even sequence represented in the form of (A) above, then the sequence may not appear to be even and you will get a complex result after transforming. Fred
Fred Marshall wrote:
..
> "Even" means "symmetrical around zero" > Generally when we do an fft or ifft, the zero time or zero frequency sample > is the first sample. > > An example of an even sequence of N=10 elements would be like this: > > A) [1 2 3 0 0 0 0 0 3 2] > > If you rotate this sequence it would look like this: > > B) [0 0 0 3 2 1 2 3 0 0] > which looks "even" but entails a delay or a frequency shift - depending on > which domain these elements belong to.
Fred, I would say both your sequences A and B are even. An even sequence of length 10 has the form (x[0], x[1], x[2], x[3], x[4], x[5], x[4], x[3], x[2], x[1]), which is satisfied by both A and B. The definition of even is: x[n] = x[N-n], n=1,2,...,N/2-1 (using zero-based indexing). For a DFT pair X = DFT(x), these two equivalences hold: i) x real and even <=> X real and even ii) x real and odd <=> X imaginary and odd. Agreed? Regards, Andor
Thanks to both of you!
I obtain a real part which is actually even with the convention
described by Andor.
The values of the imaginary part are not zero but are very small so
that I suspect they are due to precision errors.

Tania

"Andor" <an2or@mailcircuit.com> wrote in message 
news:1120635117.442542.260780@z14g2000cwz.googlegroups.com...
> Fred Marshall wrote: > .. >> "Even" means "symmetrical around zero" >> Generally when we do an fft or ifft, the zero time or zero frequency >> sample >> is the first sample. >> >> An example of an even sequence of N=10 elements would be like this: >> >> A) [1 2 3 0 0 0 0 0 3 2] >> >> If you rotate this sequence it would look like this: >> >> B) [0 0 0 3 2 1 2 3 0 0] >> which looks "even" but entails a delay or a frequency shift - depending >> on >> which domain these elements belong to. > > Fred, I would say both your sequences A and B are even. An even > sequence of length 10 has the form > > (x[0], x[1], x[2], x[3], x[4], x[5], x[4], x[3], x[2], x[1]), > > which is satisfied by both A and B. The definition of even is: > > x[n] = x[N-n], n=1,2,...,N/2-1 (using zero-based indexing). > > For a DFT pair X = DFT(x), these two equivalences hold: > > i) x real and even <=> X real and even > ii) x real and odd <=> X imaginary and odd. > > Agreed?
Andor, Yes. My example wasn't the best. It may have been better to use this one: B) [3 2 1 2 3 0 0 0 0 0 ....] .... with as many zeros at the end as you like. The nonzero part is "even" around the 3rd sample but the sequence is not even overall. It appears we accomplished our objective! Fred
Fred Marshall wrote:
...
> It appears we accomplished our objective!
Another happy customer! :-)