aesa: make "make clean" also remove output.txt
[jump.git] / 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 }