0
$\begingroup$

I am using FFTW3 to transform the autocorrelation function of a particular function. I expect the transform to be real since the autocorrelation function is symmetric with respect to negative and positive lags and its Fourier transform give the power spectrum. However, I get nonzero complex part and negative real parts after the transformation. What I may be doing wrong?

I am providing the code below:

/* Compilation and execution command:
 * gcc q1part2.c nrutil.c qtrap.c trapzd.c -lfftw3 -lm -o q1part2 && ./q1part2 
 */

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <complex.h>
#include <fftw3.h>
#include "nr.h"
#include "nrutil.h"
double h = 0.01;
double a = 0.0, b = 10.0, p;

double func(double t)
{
    return  cos(10 * t) + sin( 20 * t);
}
double intgrd(double t)
{
    return func(t) * func(t + p);
}

int main()
{
    int  N = round(((b - a) / h) + 1);
    double y[2*N - 1], t[2*N - 1]; /* t array hold time delays */
    int i;
    fftw_complex *in, *out;
    fftw_plan MYPLAN;
    in = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * (2 * N - 1));
    out = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * (2 * N - 1));
    MYPLAN = fftw_plan_dft_1d(2*N - 1, in, out, FFTW_FORWARD, FFTW_ESTIMATE);
    for(i = 0; i < 2*N - 1; i++)
    {
        /* The definition of autocorrelation from the book "A First Course on
        * Computational Physics" by Paul DeVries was used.
        */
        p = i * h - 10.0;
        t[i] = p;
        if(p >= 0)
        {
            y[i] = sqrt(1 / (2 * M_PI)) * qtrap(intgrd, a, b - p);
            in[i] = y[i];
        }
        else
        {
            y[i] = sqrt(1 / (2 * M_PI)) * qtrap(intgrd, -p, b);
            in[i] = y[i];
        }
        //printf("%f \t %f \n", t[i], y[i]);
    }
    fftw_execute(MYPLAN);
    for(i = 0; i < N; i++)
    {
        printf("%f \t %f  %f \n", i * 2*  M_PI / (2 * (b - a)), creal(out[i]),cimag(out[i])); 
    }
    fftw_destroy_plan(MYPLAN);
    fftw_free(in);
    fftw_free(out);
    return 0;
}

The function qtrap is a numerical integration algorithm provided in the code from the second edition (plain C) of the book Numerical Recipes. I am really confused about the usage of the library and would like to get pointed in the right direction.

Thanks in advance

$\endgroup$
2
$\begingroup$

If the imaginary part is very small you are probably not doing anything wrong. This is a common occurrence due to floating point quantization noise. As long as you know the result should be real and the imaginary part is small, just drop the imaginary part.

| improve this answer | |
$\endgroup$
  • $\begingroup$ Well thanks for the answer, my mistake was to include one last point which messed things up pretty bad. Apart from that this transform is the FFT of the autocorrelation function of a finite waveform namely $$ cos(10t) + sin(20t) $$ Hence its FFT should give the power spectrum I guess but although the values are real I get negative values in between. My aim is to compare the power spectrum obtained from the autocorrelation function with the power spectrum obtained from applying FFT directly on the time series of the function in the time interval from 0 to 10. $\endgroup$ – Vesnog May 18 '14 at 1:49
1
$\begingroup$

For your information, I just took one additional data point that I should not have taken. Replacing every 2*N-1 with 2*N-2 resolved the issue. My original (problematic sequence) was as follows: $$ [0,1,2,3,3,2,1,0]$$ The problem was resolved when the series became like by removing the last term. $$ [0,1,2,3,3,2,1]$$

| improve this answer | |
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.