> 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
Reply by Rune Allnor●December 11, 20062006-12-11
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
Reply by Hanno Rein●December 11, 20062006-12-11
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
Reply by Rune Allnor●December 11, 20062006-12-11
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:
>
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
Reply by Hanno Rein●December 11, 20062006-12-11
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