• API
• FAQ
• Tools
• Archive
daily pastebin goal
67%
SHARE
TWEET

# Untitled

a guest Oct 18th, 2018 67 Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
1. #include <stdio.h>
2. #include <stdlib.h>
3. #include <math.h>
4. #define FOREVER while(1){ (void)0; }
5. #define PI 3.141592654f
6. //this marco only usable by FFT_Compute() function assuming we already have var realComp, imagineComp defined
7. //four inputs: real is real component value, imagine is imaginary component value
8. //kn is muplication of time domain sample location n and frequency domain sample location k
9. // sum is total number of samples in current stage (stageNumber*2)
10. #define EXP_NEG_TWOPI_JAY_N_K_OVER_N(real, imagine, kn,sum) { realComp = real*cos(-2.0f*PI*kn/(float)sum)-imagine*sin(-2.0f*PI*kn/(float)sum); \
11.                                                               imagineComp = imagine*cos(-2.0f*PI*kn/(float)sum)+real*sin(-2.0f*PI*kn/(float)sum);}
12.
13. typedef unsigned int UINT;
14.
15. //given a M (M<32) bits number(0 - 2^M) revert bits
16. void reverseInPlace(UINT * input, UINT m)
17. {
18.     UINT out = 0, scrape = 0, k;
19.     for (k=m; k>=1; k--) {
20.         //right shift till see the leftmost bit
21.         scrape =(*input)>>(k-1) & 0x00000001;
22.         //out is assembled bit by bit from leftmost to rightmost with respect to original number
23.             out += ((scrape == 0) ?  0 : (1<<(m-k)));
24.     }
25.
26.     *input = out;
27. }
28.
29. //FFT_Compute accepts two float array of length 2^m each represent Complex Number array in time domain, real component in realArray
30. //imaginary component in imagineArray, Output Arrays in frequency domain are stored in place repectively
31. //computational complexity O(nlogn), n is the length of arrays
32. void FFT_Compute(float *realArray, float *imagineArray, int m)
33. {
34.     //two temp variables used in computation
35.     float realComp, imagineComp;
36.     //length of our array is 2^m
37.     UINT len = 0x00000001 << (m);
38.     UINT k=0, j=0, temp, magicMask;
39.     //scrambleIndex array only used in stage 0
40.     UINT *scrambleIndex;
41.     //temp variable used for storing intermediate results
42.     float *realTemp, *imgTemp;
43.
44.     if (m >= 32) {
45.         printf("We don't support bit length greater than 31");
46.         return;
47.     }
48.
49.     scrambleIndex = (UINT *)malloc(sizeof(UINT) * len);
50.     realTemp = (float *)malloc(sizeof(float) * len);
51.     imgTemp = (float *)malloc(sizeof(float) * len);
52.
53.     for (k=0; k<len; k++) {
54.         temp = k;
55.         reverseInPlace(&temp, m);
56.         scrambleIndex[k] =  temp;
57.     }
58.     //compute output for each stage, intermediate result stored in place
59.     //j is the stage number, k is the index to sample location
60.
61.     //stage zero is a bit special cannot be written in general loop
62.     //the following loop just for stage 0
63.     for (k=0; k<len; k++) {
64.         realTemp[k] = realArray[scrambleIndex[k]];
65.         imgTemp[k] = imagineArray[scrambleIndex[k]];
66.         //here we only have zero
67.         /* dont need following since Weight is Identity Matrix
68.         EXP_NEG_TWOPI_JAY_N_K_OVER_N(realTemp[k], imgTemp[k], 0, len)
69.
70.         if ((k&0x00000001) != 0) //odd number need to multiply weight
71.         {
72.             realTemp[k] = realComp;
73.             imgTemp[k] = imagineComp;
74.         } */
75.     } //copy to output array finish stage 0
76.     for (k=0; k<len; k++) {
77.         realArray[k] = realTemp[k];
78.         imagineArray[k] = imgTemp[k];
79.     }
80.     //now start from stage 1 til finish
81.     for (j=1; j<=m; j++) {
82.         //first sum then apply weight
83.         //magicMask to get butterfly location
85.         //compute each sample location in freq domain
86.         for (k=0; k<len; k++) {
87.             //adding counter part of sum instead of comparison (<) may also use & (eg. k&magicMask == 0? 1: -1)
88.             realTemp[k] = ( k < (k^magicMask) ? 1.f :-1.f) * realArray[k] + realArray[k^magicMask];
89.             imgTemp[k] = ( k < (k^magicMask) ? 1.f :-1.f) * imagineArray[k] + imagineArray[k^magicMask];
90.
91.             //applying weight only odd subsequence needs weight
92.             //donot apply weight in last stage
93.             if (j!=m && (k > (k^(magicMask<<1)))) {
94.                 //apply weight -- by left shift then right shift, we can zero out j leftmost bits
96.                 realTemp[k] = realComp;
97.                     imgTemp[k] = imagineComp;
98.             }
99.         }
100.         //copy result over finish this stage
101.         for (k=0; k<len; k++) {
102.             realArray[k] = realTemp[k];
103.             imagineArray[k] = imgTemp[k];
104.             }
105.     }
106.
107.     free(scrambleIndex);
108.     free(realTemp);
109.     free(imgTemp);
110. }
111.
112. /*
113.    Direct fourier transform for comparison purpose
114.    used for verify our Implementation is correct
115. */
116. int DFT(int dir,int m,float *x1,float *y1)
117. {
118.    long i,k;
119.    float arg;
120.    float cosarg,sinarg;
121.    float *x2,*y2;
122.
123.    x2 = (float *)malloc(m*sizeof(float));
124.    y2 = (float *)malloc(m*sizeof(float));
125.    if (x2 == NULL || y2 == NULL)
126.       return(-1);
127.
128.    for (i=0;i<m;i++) {
129.       x2[i] = 0;
130.       y2[i] = 0;
131.       arg = - dir * 2.0 * 3.141592654 * (float)i / (float)m;
132.       for (k=0;k<m;k++) {
133.          cosarg = cos(k * arg);
134.          sinarg = sin(k * arg);
135.          x2[i] += (x1[k] * cosarg - y1[k] * sinarg);
136.          y2[i] += (x1[k] * sinarg + y1[k] * cosarg);
137.       }
138.    }
139.
140.    /* Copy the data back */
141.    dir = 1;
142.    if (dir == 1) {
143.       for (i=0;i<m;i++) {
144.          x1[i] = x2[i];
145.          y1[i] = y2[i];
146.       }
147.    } else {
148.       for (i=0;i<m;i++) {
149.          x1[i] = x2[i];
150.          y1[i] = y2[i];
151.       }
152.    }
153.
154.    free(x2);
155.    free(y2);
156.    return(0);
157. }
158.
159. void main()
160. {
161.     float testReal[32], testImagine[32];
162.     int bitCount = 5, i, length = 32;
163.     for (i=0; i<length; i++) {
164.         testReal[i] = sin(2.0f*PI*(float)i/(float)length);
165.         testImagine[i] = 0;
166.     }
167.
168.     printf("usage: %d(input number) %d(total bit count)\n");
169.
170.     FFT_Compute(testReal, testImagine, bitCount);
171.     //DFT(1, length, testReal, testImagine);
172.     for (i=0; i<length; i++) {
173.         printf("Complex Number for sample %d: Real = %f Imaginary = %f\n", i, testReal[i], testImagine[i]);
174.     }
175.
176.     FOREVER
177. }
RAW Paste Data
We use cookies for various purposes including analytics. By continuing to use Pastebin, you agree to our use of cookies as described in the Cookies Policy.

Top