daily pastebin goal
78%
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
  84.         magicMask = (0x00000001<<(j-1));
  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
  95.                 EXP_NEG_TWOPI_JAY_N_K_OVER_N(realTemp[k], imgTemp[k], ((k<<(32-j))>>(32-j)), (magicMask<<2));
  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. OK, I Understand
 
Top