>> On 01/09/2011 05:11 AM, ntiy wrote:
>>> Hi DSP guys.
>>>
>>> I know there were a few threads with this question already. I looked
>>> through them and I am still stuck. Can't figure out what could be wrong
> and
>>> my eyes and head hurt. So, I thought may be someone helps...
>>>
>>> I am trying to compute derivative with fftw. Here are key points of my
> C
>>> code (well, that's actually the whole code):
(code snipped)
>>> I do not get right result.
>>> Code is trivial, right? I take function, convert it to F-space,
> multiply
>>> F-components by -om^2 (properly aliased) and then transform the result
> back
>>> to direct space. But it does not seem to give right result. I see 3
>>> possible reasons:
>>> 1) I don't use fftw in a proper way (wrong omegas for example?)
>>> 2) there is non-related to fftw bug
>>> 3) I am stupid
>>>
>>> If it is 1) or 2) please help. I am ready to smash my head against the
>>> wall. If it's 3) then, of course nothing helps...
>>
>> You haven't told us what the right result is, and what you're actually
>> getting.
>>
>> Multiplying by the square of the frequency should get you the double
>> derivative, not the 1st derivative -- is that what you want?
>>
>> You're not using any frequency limiting -- are you sure you're not just
>> amplifying the heck out of your noise?
>>
>> In Scilab:
>>
>> -->N = 4096; n = 0:N-1;
>>
>> -->nn = pmodulo(n + N/2, N) - N/2;
>>
>> -->x = exp(-nn.^2 / 100000); xd = ifft(fft(x) .* (-nn .^ 2));
>>
>> -->clf; plot2d(n, [x' xd'])
>>
>> gives me what looks like a very clean double derivative of the input
>> Gaussian.
>>
>> --
>>
>> Tim Wescott
>> Wescott Design Services
>> http://www.wescottdesign.com
>>
>> Do you need to implement control loops in software?
>> "Applied Control Theory for Embedded Systems" was written for you.
>> See details at http://www.wescottdesign.com/actfes/actfes.html
>>
>
> Thanks for reply.
> 1) sorry for confusion. You are right I am looking for second derivative:
>
> a) f_x = 1/N \sum_{\omega}f_{\omega}e^{i\omega x} -- definition of
> fourier transform
> now I do differentiation:
> b) d^2f(x)/dx^2 = 1/N \sum_{\omega}f_{\omega}e^{i\omega x}*(i\omega)^2
>
> So far that is precise expression.
So far that is _not_ a precise expression! Your definition in both step
a and step b is of the discrete Fourier transform. But the definition
(d^2/dt^2) f(x) = F{-w^2 * X(w)}
is only true for the continuous Fourier transform, which implies a
continuous x. Your x is a sampled version of a Gaussian, which is not
defined over continuous time.
> The only error that get's into play that
> I see is that FFT is an approximation and not precise. But
> errors should be of the order 10^-16 and smaller.
An approximation of what? The FFT is an exact equivalent to a DFT, just
rearranged so that when you do the computation it'll be faster. An FFT
_implementation_ is an approximation, yes, with errors on the order of
10^-16 for some computer systems (Scilab is often a bit worse).
Sampling a continuous-time signal gives an approximation of the original
signal that can range anywhere from complete (if the signal is
bandlimited per Nyquist, Shannon, etc.) to atrocious. Even then it's
not "exact" -- I said "complete" because under certain narrow
circumstances you can capture all of the _information_ in the original
signal, but the sampled-time and continuous-time domains are just two
different things: comparing a continuous-time signal to a sampled-time
signal isn't comparing apples and oranges, it's comparing apples to
rocks with pictures of apples carved in them.
> Am I wrong here? What is "frequency
> limitation"? I seem to sum over all spectra...
So sampling can introduce some HUGE approximations, depending on what
the signal is and what you're planning on doing with the samples. One
of the places where the approximation really, truly breaks down is at
high frequencies -- specifically, in this case, at frequencies close to
the Nyquist rate. If your highest gain in the system is at the Nyquist
rate, then chances are high that you've got problems.
The three errors that stick out the most for me are:
* You are failing to make the distinction between a sampled-time signal
and a continuous-time one.
* You are failing to make the distinction between a discrete Fourier
transform and a continuous Fourier transform.
* You are assuming that differentiation exists in sampled time --
strictly speaking, it doesn't.
I hope this helps.
--
Tim Wescott
Wescott Design Services
http://www.wescottdesign.com
Do you need to implement control loops in software?
"Applied Control Theory for Embedded Systems" was written for you.
See details at http://www.wescottdesign.com/actfes/actfes.html
Reply by Fred Marshall●January 10, 20112011-01-10
On 1/10/2011 10:26 AM, Fred Marshall wrote:
>
> Once more, in order to do the circular convolution you're trying to do,
> you need to append zeros to give the convolution room: So the unit
> sample response would be N samples followed by N zeros for 2N samples.
> The same thing goes for the signal being operated on.....
>
> Fred
Also, you need to give time for the transient response of the "filter"
you're using. There will be a leading and lagging edge.
Assuming the signal is of N samples, appended to 2N samples with zeros,
and,
assuming the time kernel of the "filter" is also N samples, appended to
2N samples with zeros,
and
assuming that they both start at t=0,
then
the "filter" introduces a delay of N/2
and
the output of the filter is transient for N/2 samples
it is in steady state more or less for N samples
is in transient again for N/2 samples.
So, if there are 2N samples overall, you will have N samples that are
more or less the result you want.... and, if time is aligned as above
they will be in the middle from N/2 to 3N/2.
Fred
Reply by Fred Marshall●January 10, 20112011-01-10
On 1/10/2011 7:33 AM, ntiy wrote:
>
> So far that is precise expression. The only error that get's into play that
> I see is that FFT is an approximation and not precise. But errors should be
> of the order 10^-16 and smaller. Am I wrong here? What is "frequency
> limitation"? I seem to sum over all spectra...
>
> Thanks,
> Ntiy.
Well, a couple of things:
- the FFT isn't so much an approximation or imprecise unless you get
down into the quantization errors related to the word length you're
using. Normally that isn't a huge deal so I wouldn't start out by
"going there".
- the expression you give isn't precise because you are *also* using
finite length sequences which have ramifications that are normally not
ignored.
You are computing a frequency domain multiplication of two sequences
which amounts to a time domain circular convolution of their inverse
transforms .. the time domain sequences that you're starting with.
Multiply in frequency .... circular convolution in time. They are
directly equivalent.
I might suggest that you do the time domain convolution so that you can
see what's going on. But then you're going to have to have a unit
sample response of a "differentiator". That's why I suggested using a
filter design program or algorithm to get one that's likely to work.
There is really no way around this.
If it bothers you, try this:
Compute the inverse transform of just the om^2 term. What do you get?
Is it what you expect? Likely not.
Here is a picture of the frequency magnitude response of a
differentiator that works:
o o
oo oo
| o o o o |
| o o o o |
| o o o o |
| o |o o| o |
| o o o o |
| o o o o |
| o o o o |
| o o o o |
| o o o o |
| o |o o| o |
| o o o o |
| o o o o |
| o o o o |
| o o o o |
| o o o o |
| o |o o| o |
| o o o o |
| o o o o |
|o o o o|
+-----------------------------------------------------+
f=0 fs/2 fs
It goes through zero at fs/2 because the sign flips there. And, it's
zero around fs/2 because there needs to be a practical limit in its
frequency content. That's the "frequency limitation" mentioned.
Here is a cartoon of the time domain unit sample response:
o
o
| oo |
| oo |
| o o |
| o o |
| o o |
| o o o |
| o o o |
| o o o o o o |
|oo o o o o o o o o o o o o |
o--o-o-o--oo-o--oo---o-----o------o---oo--o-oo--o-o-o--
| o o o o o o o o o o o o oo
| o o o o o o |
| o o o |
| o o o |
| o o |
| o o |
| o o |
| o o |
| o o |
oo
t=0 oo t=(N-1)*T
o
o
t>
Note that the unit sample response of a differentiator (well, since
we're dealing with sampled data, it's really a "differencer" would be:
[0 0 0 0 0 0 0 ........... 1 -1 0 0 0 0 0 0 0.......]
And, since it has to be bandlimited, we effectively pass it through a
lowpass filter which gets the zero-valued region around fs/2 and gets
those wavy parts in the unit sample response in time.
Once more, in order to do the circular convolution you're trying to do,
you need to append zeros to give the convolution room: So the unit
sample response would be N samples followed by N zeros for 2N samples.
The same thing goes for the signal being operated on.....
Fred
Reply by ntiy●January 10, 20112011-01-10
>On 01/09/2011 05:11 AM, ntiy wrote:
>> Hi DSP guys.
>>
>> I know there were a few threads with this question already. I looked
>> through them and I am still stuck. Can't figure out what could be wrong
and
>> my eyes and head hurt. So, I thought may be someone helps...
>>
>> I am trying to compute derivative with fftw. Here are key points of my
C
>> code (well, that's actually the whole code):
>>
>> main()
>> {
>> ...
>> ...
>>
>> N = 128;
>> a0 = 1;
>> alpha = 1;
>> L = 80;
>> dx = L/N;
>> dom = 2*PI/L;
>> ommax = dom*(N/2);
>> size_t size = sizeof( fftw_complex )*N;
>>
>> fftw_complex *c, *c_xx, *fc, *fc_xx;
>> fftw_plan c_to_fc, fc_xx_to_c_xx;
>>
>> om = ( double *) malloc( sizeof(double)*N );
>> c = ( fftw_complex* ) fftw_malloc( size ); // c in Re
>> c_xx = ( fftw_complex* ) fftw_malloc( size ); // c_xx in Re
>> fc = ( fftw_complex* ) fftw_malloc( size ); // c in Fr
>> fc_xx = ( fftw_complex* ) fftw_malloc( size ); // laplacian
>>
>> c_to_fc = fftw_plan_dft_1d(N, c, fc,
>> FFTW_FORWARD, FFTW_ESTIMATE);
>> fc_xx_to_c_xx = fftw_plan_dft_1d(N, fc_xx, c_xx,
>> FFTW_BACKWARD, FFTW_ESTIMATE);
>>
>>
>> // initialize input for FFT procedure.
>> // I take simple gaussian as an input
>> for(i=0; i<N; i++)
>> {
>> x = -L/2 + i*dx;
>> c[i][0] = exp( -a0*x*x );
>> c[i][1] = 0;
>> }
>>
>> // init matrix of omega's. this piece seems to be crucial --
>> // how to rearrange frequencies so that they are in a right order
>> for (i=0; i<N; i++)
>> {
>> if(i<N/2) { om[i] = dom*i; }
>> else { om[i] = dom*(-N + i); }
>> }
>> om[N/2]=0;
>> // this gives me freqs in the order:
>> // dom*[0, 1, 2, ... N/2-1,0,-N/2+1,-N/2+2,...,-1]
>>
>>
>> fftw_execute( c_to_fc ); // compute fourier of my function
>> for(i=0; i<N; i++)
>> {
>> tom2 = om[i]*om[i];
>> fc_xx[i][0] = -alpha*tom2*fc[i][0]; // compute fourier of
d2c/dx2
>> fc_xx[i][1] = -alpha*tom2*fc[i][1]; // f[d2c/dx2] = -om*om f[c]
>> }
>> fftw_execute( fc_xx_to_c_xx ); // go back to direct space
>>
>> for(i=0; i<N; i++) // print result
>> {
>> x = -L/2 + i*dx;
>> printf("%f\t%f\t%e\n", x, c_xx[i][0]/N, c_xx[i][1]/N);
>> }
>>
>> // free memory after that
>> ...
>> ..
>> }
>>
>> I do not get right result.
>> Code is trivial, right? I take function, convert it to F-space,
multiply
>> F-components by -om^2 (properly aliased) and then transform the result
back
>> to direct space. But it does not seem to give right result. I see 3
>> possible reasons:
>> 1) I don't use fftw in a proper way (wrong omegas for example?)
>> 2) there is non-related to fftw bug
>> 3) I am stupid
>>
>> If it is 1) or 2) please help. I am ready to smash my head against the
>> wall. If it's 3) then, of course nothing helps...
>
>You haven't told us what the right result is, and what you're actually
>getting.
>
>Multiplying by the square of the frequency should get you the double
>derivative, not the 1st derivative -- is that what you want?
>
>You're not using any frequency limiting -- are you sure you're not just
>amplifying the heck out of your noise?
>
>In Scilab:
>
>-->N = 4096; n = 0:N-1;
>
>-->nn = pmodulo(n + N/2, N) - N/2;
>
>-->x = exp(-nn.^2 / 100000); xd = ifft(fft(x) .* (-nn .^ 2));
>
>-->clf; plot2d(n, [x' xd'])
>
>gives me what looks like a very clean double derivative of the input
>Gaussian.
>
>--
>
>Tim Wescott
>Wescott Design Services
>http://www.wescottdesign.com
>
>Do you need to implement control loops in software?
>"Applied Control Theory for Embedded Systems" was written for you.
>See details at http://www.wescottdesign.com/actfes/actfes.html
>
Thanks for reply.
1) sorry for confusion. You are right I am looking for second derivative:
a) f_x = 1/N \sum_{\omega}f_{\omega}e^{i\omega x} -- definition of
fourier transform
now I do differentiation:
b) d^2f(x)/dx^2 = 1/N \sum_{\omega}f_{\omega}e^{i\omega x}*(i\omega)^2
So far that is precise expression. The only error that get's into play that
I see is that FFT is an approximation and not precise. But errors should be
of the order 10^-16 and smaller. Am I wrong here? What is "frequency
limitation"? I seem to sum over all spectra...
Thanks,
Ntiy.
Reply by Tim Wescott●January 9, 20112011-01-09
On 01/09/2011 05:11 AM, ntiy wrote:
> Hi DSP guys.
>
> I know there were a few threads with this question already. I looked
> through them and I am still stuck. Can't figure out what could be wrong and
> my eyes and head hurt. So, I thought may be someone helps...
>
> I am trying to compute derivative with fftw. Here are key points of my C
> code (well, that's actually the whole code):
>
> main()
> {
> ...
> ...
>
> N = 128;
> a0 = 1;
> alpha = 1;
> L = 80;
> dx = L/N;
> dom = 2*PI/L;
> ommax = dom*(N/2);
> size_t size = sizeof( fftw_complex )*N;
>
> fftw_complex *c, *c_xx, *fc, *fc_xx;
> fftw_plan c_to_fc, fc_xx_to_c_xx;
>
> om = ( double *) malloc( sizeof(double)*N );
> c = ( fftw_complex* ) fftw_malloc( size ); // c in Re
> c_xx = ( fftw_complex* ) fftw_malloc( size ); // c_xx in Re
> fc = ( fftw_complex* ) fftw_malloc( size ); // c in Fr
> fc_xx = ( fftw_complex* ) fftw_malloc( size ); // laplacian
>
> c_to_fc = fftw_plan_dft_1d(N, c, fc,
> FFTW_FORWARD, FFTW_ESTIMATE);
> fc_xx_to_c_xx = fftw_plan_dft_1d(N, fc_xx, c_xx,
> FFTW_BACKWARD, FFTW_ESTIMATE);
>
>
> // initialize input for FFT procedure.
> // I take simple gaussian as an input
> for(i=0; i<N; i++)
> {
> x = -L/2 + i*dx;
> c[i][0] = exp( -a0*x*x );
> c[i][1] = 0;
> }
>
> // init matrix of omega's. this piece seems to be crucial --
> // how to rearrange frequencies so that they are in a right order
> for (i=0; i<N; i++)
> {
> if(i<N/2) { om[i] = dom*i; }
> else { om[i] = dom*(-N + i); }
> }
> om[N/2]=0;
> // this gives me freqs in the order:
> // dom*[0, 1, 2, ... N/2-1,0,-N/2+1,-N/2+2,...,-1]
>
>
> fftw_execute( c_to_fc ); // compute fourier of my function
> for(i=0; i<N; i++)
> {
> tom2 = om[i]*om[i];
> fc_xx[i][0] = -alpha*tom2*fc[i][0]; // compute fourier of d2c/dx2
> fc_xx[i][1] = -alpha*tom2*fc[i][1]; // f[d2c/dx2] = -om*om f[c]
> }
> fftw_execute( fc_xx_to_c_xx ); // go back to direct space
>
> for(i=0; i<N; i++) // print result
> {
> x = -L/2 + i*dx;
> printf("%f\t%f\t%e\n", x, c_xx[i][0]/N, c_xx[i][1]/N);
> }
>
> // free memory after that
> ...
> ..
> }
>
> I do not get right result.
> Code is trivial, right? I take function, convert it to F-space, multiply
> F-components by -om^2 (properly aliased) and then transform the result back
> to direct space. But it does not seem to give right result. I see 3
> possible reasons:
> 1) I don't use fftw in a proper way (wrong omegas for example?)
> 2) there is non-related to fftw bug
> 3) I am stupid
>
> If it is 1) or 2) please help. I am ready to smash my head against the
> wall. If it's 3) then, of course nothing helps...
You haven't told us what the right result is, and what you're actually
getting.
Multiplying by the square of the frequency should get you the double
derivative, not the 1st derivative -- is that what you want?
You're not using any frequency limiting -- are you sure you're not just
amplifying the heck out of your noise?
In Scilab:
-->N = 4096; n = 0:N-1;
-->nn = pmodulo(n + N/2, N) - N/2;
-->x = exp(-nn.^2 / 100000); xd = ifft(fft(x) .* (-nn .^ 2));
-->clf; plot2d(n, [x' xd'])
gives me what looks like a very clean double derivative of the input
Gaussian.
--
Tim Wescott
Wescott Design Services
http://www.wescottdesign.com
Do you need to implement control loops in software?
"Applied Control Theory for Embedded Systems" was written for you.
See details at http://www.wescottdesign.com/actfes/actfes.html
Reply by Fred Marshall●January 9, 20112011-01-09
On 1/9/2011 5:11 AM, ntiy wrote:
> Hi DSP guys.
>
> I know there were a few threads with this question already. I looked
> through them and I am still stuck. Can't figure out what could be wrong and
> my eyes and head hurt. So, I thought may be someone helps...
>
> I am trying to compute derivative with fftw. Here are key points of my C
> code (well, that's actually the whole code):
>
> main()
> {
> ...
> ...
>
> N = 128;
> a0 = 1;
> alpha = 1;
> L = 80;
> dx = L/N;
> dom = 2*PI/L;
> ommax = dom*(N/2);
> size_t size = sizeof( fftw_complex )*N;
>
> fftw_complex *c, *c_xx, *fc, *fc_xx;
> fftw_plan c_to_fc, fc_xx_to_c_xx;
>
> om = ( double *) malloc( sizeof(double)*N );
> c = ( fftw_complex* ) fftw_malloc( size ); // c in Re
> c_xx = ( fftw_complex* ) fftw_malloc( size ); // c_xx in Re
> fc = ( fftw_complex* ) fftw_malloc( size ); // c in Fr
> fc_xx = ( fftw_complex* ) fftw_malloc( size ); // laplacian
>
> c_to_fc = fftw_plan_dft_1d(N, c, fc,
> FFTW_FORWARD, FFTW_ESTIMATE);
> fc_xx_to_c_xx = fftw_plan_dft_1d(N, fc_xx, c_xx,
> FFTW_BACKWARD, FFTW_ESTIMATE);
>
>
> // initialize input for FFT procedure.
> // I take simple gaussian as an input
> for(i=0; i<N; i++)
> {
> x = -L/2 + i*dx;
> c[i][0] = exp( -a0*x*x );
> c[i][1] = 0;
> }
>
> // init matrix of omega's. this piece seems to be crucial --
> // how to rearrange frequencies so that they are in a right order
> for (i=0; i<N; i++)
> {
> if(i<N/2) { om[i] = dom*i; }
> else { om[i] = dom*(-N + i); }
> }
> om[N/2]=0;
> // this gives me freqs in the order:
> // dom*[0, 1, 2, ... N/2-1,0,-N/2+1,-N/2+2,...,-1]
>
>
> fftw_execute( c_to_fc ); // compute fourier of my function
> for(i=0; i<N; i++)
> {
> tom2 = om[i]*om[i];
> fc_xx[i][0] = -alpha*tom2*fc[i][0]; // compute fourier of d2c/dx2
> fc_xx[i][1] = -alpha*tom2*fc[i][1]; // f[d2c/dx2] = -om*om f[c]
> }
> fftw_execute( fc_xx_to_c_xx ); // go back to direct space
>
> for(i=0; i<N; i++) // print result
> {
> x = -L/2 + i*dx;
> printf("%f\t%f\t%e\n", x, c_xx[i][0]/N, c_xx[i][1]/N);
> }
>
> // free memory after that
> ...
> ..
> }
>
> I do not get right result.
> Code is trivial, right? I take function, convert it to F-space, multiply
> F-components by -om^2 (properly aliased) and then transform the result back
> to direct space. But it does not seem to give right result. I see 3
> possible reasons:
> 1) I don't use fftw in a proper way (wrong omegas for example?)
> 2) there is non-related to fftw bug
> 3) I am stupid
>
> If it is 1) or 2) please help. I am ready to smash my head against the
> wall. If it's 3) then, of course nothing helps...
>
> Thanks.
>
>
>
There's much I don't understand about what you're doing and I don't care
to decipher the code really. Just what is the algorithm that you're
trying to implement .. in simple terms?
I'll guess:
You have f(t).
You compute F(w).
You compute jw*F(w)except at N/2 you've set jw to zero?
You inverse transform jw*F(w)
If everything were continuous and infinite then that follows the math.
But, that's not the case here.
Might one suggest that you design a differentiator using Parks-McClellan
or some other suitable FIR design program .. of the proper length (N) ..
and use the FT of that filter for convolution by multiplication? That's
less risky.
Also, when you're done with the frequency domain multiplication, are
there enough samples to avoid time domain aliasing/overlap of the
results? That is, if the original time sequence is of length N then the
FFT is of length N. And, if the differentiator is of length N (as it
must be) then the temporal convolution will be of length 2N-1. Have you
accounted for that?
Normally one might zero-pad from N to 2N in both cases so that there's
room for the convolution results.
I didn't try to follow all the index stuff you're doing. Seems like it
should be more straightforward. But if you need to orient the
differentiator to be symmetrical in time around zero then one could
understand doing something like that. That is, 1/2 of the nonzero
elements would be from zero to N/2 and 1/2 of the nonzero elements would
be from N/2 to N-1 and the zero padding would be around N/2. And,
perhaps the same for the signal if you want to center it at zero.
Fred
Reply by ntiy●January 9, 20112011-01-09
Hi DSP guys.
I know there were a few threads with this question already. I looked
through them and I am still stuck. Can't figure out what could be wrong and
my eyes and head hurt. So, I thought may be someone helps...
I am trying to compute derivative with fftw. Here are key points of my C
code (well, that's actually the whole code):
main()
{
...
...
N = 128;
a0 = 1;
alpha = 1;
L = 80;
dx = L/N;
dom = 2*PI/L;
ommax = dom*(N/2);
size_t size = sizeof( fftw_complex )*N;
fftw_complex *c, *c_xx, *fc, *fc_xx;
fftw_plan c_to_fc, fc_xx_to_c_xx;
om = ( double *) malloc( sizeof(double)*N );
c = ( fftw_complex* ) fftw_malloc( size ); // c in Re
c_xx = ( fftw_complex* ) fftw_malloc( size ); // c_xx in Re
fc = ( fftw_complex* ) fftw_malloc( size ); // c in Fr
fc_xx = ( fftw_complex* ) fftw_malloc( size ); // laplacian
c_to_fc = fftw_plan_dft_1d(N, c, fc,
FFTW_FORWARD, FFTW_ESTIMATE);
fc_xx_to_c_xx = fftw_plan_dft_1d(N, fc_xx, c_xx,
FFTW_BACKWARD, FFTW_ESTIMATE);
// initialize input for FFT procedure.
// I take simple gaussian as an input
for(i=0; i<N; i++)
{
x = -L/2 + i*dx;
c[i][0] = exp( -a0*x*x );
c[i][1] = 0;
}
// init matrix of omega's. this piece seems to be crucial --
// how to rearrange frequencies so that they are in a right order
for (i=0; i<N; i++)
{
if(i<N/2) { om[i] = dom*i; }
else { om[i] = dom*(-N + i); }
}
om[N/2]=0;
// this gives me freqs in the order:
// dom*[0, 1, 2, ... N/2-1,0,-N/2+1,-N/2+2,...,-1]
fftw_execute( c_to_fc ); // compute fourier of my function
for(i=0; i<N; i++)
{
tom2 = om[i]*om[i];
fc_xx[i][0] = -alpha*tom2*fc[i][0]; // compute fourier of d2c/dx2
fc_xx[i][1] = -alpha*tom2*fc[i][1]; // f[d2c/dx2] = -om*om f[c]
}
fftw_execute( fc_xx_to_c_xx ); // go back to direct space
for(i=0; i<N; i++) // print result
{
x = -L/2 + i*dx;
printf("%f\t%f\t%e\n", x, c_xx[i][0]/N, c_xx[i][1]/N);
}
// free memory after that
...
..
}
I do not get right result.
Code is trivial, right? I take function, convert it to F-space, multiply
F-components by -om^2 (properly aliased) and then transform the result back
to direct space. But it does not seem to give right result. I see 3
possible reasons:
1) I don't use fftw in a proper way (wrong omegas for example?)
2) there is non-related to fftw bug
3) I am stupid
If it is 1) or 2) please help. I am ready to smash my head against the
wall. If it's 3) then, of course nothing helps...
Thanks.