aesa/src/FFT_Radix2_DIF.c

   1 #include "FFT_Radix2_DIF.h"
   2
   3 /************************************************************************************/
   4
   5 void FFT_Radix2_DIF_init(DOLProcess * p)
   6 {
   7         // int i,j;
   8
   9         // for (i=0; i<NFFT; i=i+1)
  10         // {
  11         //      p->local->Samples_Inp[i].real = 0.0;    // Buffer A of NFFT size
  12         //      p->local->Samples_Inp[i].imag = 0.0;
  13         //      p->local->Samples_Out[i].real = 0.0;    // Buffer B of NFFT size
  14         //      p->local->Samples_Out[i].imag = 0.0;
  15         // }
  16
  17         p->local->pow_2[0]=1;
  18         int i;
  19         for (i=1; i < MAX_POW; i++)
  20         {
  21                 p->local->pow_2[i] = p->local->pow_2[i-1]*2;
  22         }
  23
  24         // p->local->sample_counter = 0;
  25         p->local->n_ffts = 0;
  26         p->local->n_fft_batches = 0;
  27         p->local->n_iterations = 0;
  28         p->local->excount = 0;
  29
  30 }
  31
  32 /************************************************************************************/
  33
  34 int FFT_Radix2_DIF_fire(DOLProcess * p)
  35 {
  36         // ComplexNumber next_sample;
  37         ComplexNumber samples_in[NFFT];
  38         // ComplexNumber samples_out[NFFT];
  39         p->local->excount++;
  40
  41
  42     int N1,N2, k1, c;
  43
  44         DOL_read((void*)PORT_DATA_IN, &samples_in, NFFT*sizeof(ComplexNumber), p);
  45
  46         // p->local->Samples_Inp[p->local->sample_counter] = next_sample;
  47         // next_sample = p->local->Samples_Out[p->local->sample_counter];
  48
  49         // The value is calculated, the Dol write is placed as the bottom of the fire function
  50
  51         // p->local->sample_counter = p->local->sample_counter + 1;
  52         // if (p->local->sample_counter == NFFT)
  53         // {
  54     N1 = 2;
  55     N2 = NFFT/2;
  56
  57     // FFT_Radix2(&(p->local->Samples_Inp[0]), NFFT);  // Compute first stage
  58     FFT_Radix2(&(samples_in[0]), NFFT);  // Compute first stage
  59
  60
  61     // !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  62     // I don't really understand what's happening here! \/
  63         while (N2 != 1)  // Compute log NFFT stages
  64         {
  65                   for (k1 = 0; k1 < N1; k1++)
  66                 {
  67             FFT_Radix2(&(samples_in[N2*k1]), N2);
  68                 }
  69       N2 = N2/2;
  70         }
  71
  72           // Bit_Reverse_Reorder(&(p->local->Samples_Inp[0]), &(p->local->Samples_Out[0]), NFFT, p);  // Reorder data for output
  73
  74                 // p->local->sample_counter = 0;
  75                 p->local->n_ffts = p->local->n_ffts + 1;
  76                 if (p->local->n_ffts == NOF_FFT_PER_RANGE_BIN)
  77                 {
  78                         p->local->n_ffts = 0;
  79                         p->local->n_fft_batches = p->local->n_fft_batches + 1;
  80                         if (p->local->n_fft_batches == NOF_FFT_BATCHES)
  81                         {
  82                                 p->local->n_fft_batches = 0;
  83                                 p->local->n_iterations = p->local->n_iterations + 1;
  84                         }
  85                 }
  86         // }
  87                 // printf("\n FFT batches : %d\n",p->local->n_fft_batches);
  88
  89         DOL_write((void*)PORT_DATA_OUT, &samples_in, NFFT*sizeof(ComplexNumber), p);
  90
  91         // what i think is happening: is that the output is placed in sample in, but in reverse reorder,
  92         // it looks at how many inputs we have got, and it switches them. not necessary anymore!
  93
  94         if (p->local->n_iterations == NUMBER_OF_ITERATIONS)
  95         {
  96                 printf("\n\n:: FFT_Radix2_DIF process finished ::\n\n");
  97     // printf("\ttotal number of executions:  %d",p->local->excount);
  98
  99                 DOL_detach(p);
 100         }
 101
 102         return 0;
 103 }
 104
 105 /************************************************************************************/
 106
 107 void FFT_Radix2(ComplexNumber * data, int N)
 108 {
 109         int n2, k1, N1, N2;
 110         ComplexNumber W, butterfly[2];
 111
 112         N1=2;
 113         N2=N/2;
 114
 115         /** Do 2 Point DFTs */
 116         for (n2 = 0; n2 < N2; n2++)
 117         {
 118           W.real=cos(n2*2.0*PI/(double)N);
 119           W.imag=-sin(n2*2.0*PI/(double)N);
 120
 121                 /** Compute one butterfly **/
 122                 butterfly[0].real = (data[n2].real + data[N2 + n2].real);
 123                 butterfly[0].imag = (data[n2].imag + data[N2 + n2].imag);
 124                 butterfly[1].real = (data[n2].real - data[N2 + n2].real) * W.real - ((data[n2].imag - data[N2 + n2].imag) * W.imag);
 125                 butterfly[1].imag = (data[n2].imag - data[N2 + n2].imag) * W.real + ((data[n2].real - data[N2 + n2].real) * W.imag);
 126
 127                 /** In-place results */
 128                 for (k1 = 0; k1 < N1; k1++)
 129                 {
 130             data[n2 + N2*k1].real = butterfly[k1].real;
 131             data[n2 + N2*k1].imag = butterfly[k1].imag;
 132                 }
 133         }
 134 }
 135
 136 /************************************************************************************/
 137
 138 void Bit_Reverse_Reorder(ComplexNumber * input, ComplexNumber * output, int N, DOLProcess * p)
 139 {
 140         int bits, i, j, k;
 141         bits = 0;
 142
 143         for (i=0; i<MAX_POW; i++)
 144         {
 145                 if (p->local->pow_2[i] == N)
 146                 {
 147                         bits=i;
 148                 }
 149         }
 150
 151         for (i=0; i<N; i++)
 152         {
 153                 j=0;
 154                 for (k=0; k<bits; k++)
 155                 {
 156             if (i & p->local->pow_2[k])
 157                         {
 158                                 j += p->local->pow_2[bits-k-1];
 159                         }
 160                 }
 161                 if (j>=i)
 162                 {
 163       output[i].real = input[j].real;
 164       output[i].imag = input[j].imag;
 165       output[j].real = input[i].real;
 166       output[j].imag = input[i].imag;
 167                 }
 168         }
 169 }