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