Using FFTW to calculate derivatives

Hello!

I want to use FFTW to calculate derivatives. But it doesn't work. I
think I make a big mistake, but I couldn't find it. I calculate the
fourier transform, multiply the coefficents with -i*k, and transform
back. My code is:

void TestFFTW() {
		 fstream o1,o2;
		 o1.open("out1.txt", ios::out);
		 o2.open("out2.txt", ios::out);
		 int N=1000;

		 std::complex<double>*  c=new std::complex<double> [N];

		 fftw_plan p  =
fftw_plan_dft_1d(N,reinterpret_cast<fftw_complex*>(c),

reinterpret_cast<fftw_complex*>(c),
                                          FFTW_FORWARD, FFTW_ESTIMATE);
		 // Generate input
		 for (int i=0;i<N;i++){
			c[i]=sin(2.*M_PI/(double)N*(double)i*2.)+cos(2.*M_PI/(double)N*(double)i);

			o1 << (double)i/(double)N << " " << c[i].real() << endl;
		 }

		 fftw_execute(p);

		 complex<double> I(0,1.);
		 for (int i=0;i<N;i++){
			// Calculate derivative
			c[i]*=-I*(double)i/(double)N;
	  	 }

		 fftw_execute(p);

		 for (int i=0;i<N;i++)
				// renormalize
				c[i]/=(double)N;

		 // output
		 for (int i=0;i<N;i++)
		 o2 << (double)i/(double)N << " "
                          << c[i].real() << " "
                          << c[i].imag() << std::endl;

		 fftw_destroy_plan(p);
		 delete [] c;
	
		 o1.close();
		 o2.close();
}


Thanks for any help!

 Hanno

Hanno Rein skrev:
> Hello!
>
> I want to use FFTW to calculate derivatives. But it doesn't work. I
> think I make a big mistake, but I couldn't find it. I calculate the
> fourier transform, multiply the coefficents with -i*k, and transform
> back. My code is:
>

> 		 complex<double> I(0,1.);
> 		 for (int i=0;i<N;i++){
> 			// Calculate derivative
> 			c[i]*=-I*(double)i/(double)N;
> 	  	 }

You probably need to negative frequencies into account. This loop
seems correct for 0 <= i < N/2, but from N/2 onwards you may want
to multiply by something like -(N/2-i). Remember, the spectrum
needs to preserve symmetry.

Rune

Rune Allnor schrieb:

> You probably need to negative frequencies into account. This loop
> seems correct for 0 <= i < N/2, but from N/2 onwards you may want
> to multiply by something like -(N/2-i). Remember, the spectrum
> needs to preserve symmetry.

I think you are right. However, It is still not working. But I tried
the real to complex transform where you only get N/2+1 complex numbers
back. It works fine with that if I multiply every coefficent by -I*i.

Hanno

Hanno Rein skrev:
> Rune Allnor schrieb:
>
> > You probably need to negative frequencies into account. This loop
> > seems correct for 0 <= i < N/2, but from N/2 onwards you may want
> > to multiply by something like -(N/2-i). Remember, the spectrum
> > needs to preserve symmetry.
>
> I think you are right. However, It is still not working. But I tried
> the real to complex transform where you only get N/2+1 complex numbers
> back. It works fine with that if I multiply every coefficent by -I*i.

I never get the details about the symmetries right the first time
around; the correct factor might be something like -((N-1)/2+i) or
something like that. My point was merely that you can't scale
naively with spectrum index.

Rune

Rune Allnor wrote:
> You probably need to negative frequencies into account. This loop
> seems correct for 0 <= i < N/2, but from N/2 onwards you may want
> to multiply by something like -(N/2-i). Remember, the spectrum
> needs to preserve symmetry.

Its (i-N) for i > N/2, actually.  Remember that the point of aliasing
is that you can add any *integer* multiple of N to the frequency index
i  (although I'd prefer to use "i" for sqrt(-1)).  To get a
minimal-slope interpolation, you want to add the integer multiple that
minimizes the frequency, which is 0*N for i < N/2 and -1*N for i > N/2.
 And the N/2 component (for even N) should be set to zero.

But, as the original poster discovered for himself, these difficulties
are largely avoided by using the r2c (specialized real-input FFT)
interface, which doesn't compute the redundant half of the spectrum in
the first place.

Steven G. Johnson