Take the 2-minute tour ×
Stack Overflow is a question and answer site for professional and enthusiast programmers. It's 100% free, no registration required.

I have problems calculating the autocorrelation of my time series using the FFT. I know that

 Corr(g,h)_j <-> G_k x H_k*

where G_k and H_k are the discrete Fourier transform of g_j and h_j and H_k* is the complex conjugation of H_k. My Code is:

void fft_corr(float *vin1, float *vin2, float *vout, int n)
{
    static complex *A1; //working array
    static complex *A2; //working array
    static complex *B; //working array
    int i,k;

    A1  = (complex*) alloc1(n,sizeof(complex));
    A2  = (complex*) alloc1(n,sizeof(complex));
    B  = (complex*) alloc1(n,sizeof(complex));

    //copy pressure array to complex array
    for(i=0;i<n;i++){
        A1[i].r = vin1[i];
        A2[i].r = vin2[i];
        A1[i].i = 0.;
        A2[i].i = 0.;
    }

    // perform  fft inverse flag=1
    pfacc(1,n,A1);
    pfacc(1,n,A2);

    for(k=0;k<n;k++){
        B[k].r = A1[k].r * A2[k].r + A1[k].i * A2[k].i;
        B[k].i = - A1[k].r * A2[k].i + A1[k].i * A2[k].r;
    }

    pfacc(-1,n,B);

    for(i=0;i<n;i++) vout[i]=B[i].r/n; //write into output and normalize

}//#

What am I missing? The result does not look like an autocorrelation if A1 == A2.

Thanks a lot result

share|improve this question
    
anybody...........? –  MichaelScott Jul 16 '13 at 15:47
    
The code looks OK at first glance, although using a cross-correlation routine to perform autocorrelation is rather inefficient (you probably realise this). When you say that the autocorrelation looks wrong, what exactly is the problem ? –  Paul R Jul 16 '13 at 15:53
    
Hey Paul, thanks a lot for your reply. I have added a plot of the result from a n=315 long input array. –  MichaelScott Jul 16 '13 at 16:20
    
I think it's probably OK - what might be confusing you is that you should only be plotting the first N/2 points, as the FFT is conjugate-symetric about N/2 for real-valued inputs. Note also that you're doing circular correlation (assumes periodic input) - if that's not what you want then you'll need to zero-pad the input data. –  Paul R Jul 16 '13 at 16:32
1  
My best guess is that your FFT routine only support N = 2^m internally, so it's using a 512 point FFT for your 315 point data set. That would probably make your scaling factor wrong. Note that 315/512 = 0.615234375, which is close to your peak value. You might want to take a look at the documentation for pfacc. –  Paul R Jul 16 '13 at 18:14

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Browse other questions tagged or ask your own question.