comparison src/btrack_plus/accFFT.cpp @ 0:3dcbd77efc94

added files for OF project

+//
+//  accFFT.cpp
+//  AccelerateFFTtool
+//
+//  Created by Adam Stark on 17/07/2012.
+//  Copyright (c) 2012 __MyCompanyName__. All rights reserved.
+//
+
+#include "accFFT.h"
+
+
+accFFT :: accFFT(int fft_size,int type)
+{
+    fft_type = type;
+    
+    fftSize = fft_size;           
+    fftSizeOver2 = fftSize/2;
+    log2n = log2f(fftSize);  
+    log2nOver2 = log2n/2;
+    
+    if (fft_type == 0)
+    {
+        split.realp = (float *) malloc(fftSize * sizeof(float));
+        split.imagp = (float *) malloc(fftSize * sizeof(float));
+        
+        // allocate the fft object once
+        fftSetup = vDSP_create_fftsetup(log2n, FFT_RADIX2);
+        if (fftSetup == NULL) {
+            //printf("FFT Setup failed\n");
+        }
+    }
+    else if (fft_type == 1) 
+    {
+        d_split.realp = (double *) malloc(fftSize * sizeof(double));
+        d_split.imagp = (double *) malloc(fftSize * sizeof(double));
+        
+        // allocate the fft object once
+        fftSetupD = vDSP_create_fftsetupD(log2n, FFT_RADIX2);
+        if (fftSetupD == NULL) {
+            //printf("FFT Setup failed\n");
+        }
+    }
+        
+    
+    
+}
+
+accFFT :: ~accFFT()
+{
+    if (fft_type == 0)
+    {
+        free(split.realp);
+        free(split.imagp);
+        vDSP_destroy_fftsetup(fftSetup);
+    }
+    else if (fft_type == 1)
+    {
+        free(d_split.realp);
+        free(d_split.imagp);
+        vDSP_destroy_fftsetupD(fftSetupD);
+    }
+    
+    
+}
+
+void accFFT :: forward_FFT_f(float *buffer,float *real,float *imag)
+{        
+    //convert to split complex format with evens in real and odds in imag
+    vDSP_ctoz((COMPLEX *) buffer, 2, &split, 1, fftSizeOver2);
+    
+    //calc fft
+    vDSP_fft_zrip(fftSetup, &split, 1, log2n, FFT_FORWARD);
+    
+    // set Nyquist component to imaginary of 0 component
+    split.realp[fftSizeOver2] = split.imagp[0];
+    split.imagp[fftSizeOver2] = 0.0;
+    
+    // set 0 component to zero
+    split.imagp[0] = 0.0;
+    
+    // multiply by 0.5 to get correct output (to do with Apple's FFT implementation)
+    for (i = 0; i <= fftSizeOver2; i++)
+    {
+        split.realp[i] *= 0.5;
+        split.imagp[i] *= 0.5;
+    }
+    
+    // set values above N/2+1 which are complex conjugate mirror image of those below
+    for (i = fftSizeOver2 - 1;i > 0;--i)
+    {
+        split.realp[2*fftSizeOver2 - i] = split.realp[i];
+        split.imagp[2*fftSizeOver2 - i] = -1*split.imagp[i];
+        
+        //cout << split_data.realp[2*fftSizeOver2 - i] << "   " << split_data.imagp[2*fftSizeOver2 - i] << "i" << endl;
+    }
+    
+    for (i = 0;i < fftSize;i++)
+    {
+        real[i] = split.realp[i];
+        imag[i] = split.imagp[i];
+    }
+}
+
+
+
+void accFFT :: forward_FFT_d(double *buffer,fft_complex *out)
+{        
+    //convert to split complex format with evens in real and odds in imag
+    vDSP_ctozD((DOUBLE_COMPLEX *) buffer, 2, &d_split, 1, fftSizeOver2);
+    
+    //calc fft
+    vDSP_fft_zripD(fftSetupD, &d_split, 1, log2n, FFT_FORWARD);
+    
+    // set Nyquist component to imaginary of 0 component
+    d_split.realp[fftSizeOver2] = d_split.imagp[0];
+    d_split.imagp[fftSizeOver2] = 0.0;
+    
+    // set 0 component to zero
+    d_split.imagp[0] = 0.0;
+    
+    // multiply by 0.5 to get correct output (to do with Apple's FFT implementation)
+    for (i = 0; i <= fftSizeOver2; i++)
+    {
+        d_split.realp[i] *= 0.5;
+        d_split.imagp[i] *= 0.5;
+    }
+    
+    // set values above N/2+1 which are complex conjugate mirror image of those below
+    for (i = fftSizeOver2 - 1;i > 0;--i)
+    {
+        d_split.realp[2*fftSizeOver2 - i] = d_split.realp[i];
+        d_split.imagp[2*fftSizeOver2 - i] = -1*d_split.imagp[i];
+    }
+    
+    for (i = 0;i < fftSize;i++)
+    {
+        out[i][0] = d_split.realp[i]; 
+        out[i][1] = d_split.imagp[i];
+    }
+}