Forums

fftw and derivative. once again

Started by ntiy January 9, 2011
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.  



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
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
>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.
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
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
On 01/10/2011 07:33 AM, ntiy wrote:
>> 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