Subversion Repositories tcp_socket

/* fft.c */
/* Jonathan P Dawson */
/* 2013-12-23 */
 
#include "math.h"
 
/*globals*/
const int n = 1024;
const int m = 10;
float twiddle_step_real[m];
float twiddle_step_imaginary[m];
 
 
/*calculate twiddle factors and store them*/
void calculate_twiddles(){
    unsigned stage, subdft_size, span;
    for(stage=1; stage<=m; stage++){
        subdft_size = 1 << stage;
        span = subdft_size >> 1;
        twiddle_step_real[stage] = cos(M_PI/span);
        twiddle_step_imaginary[stage] = -sin(M_PI/span);
    }
}
 
/*bit reverse*/
unsigned bit_reverse(unsigned forward){
    unsigned reversed=0;
    unsigned i;
    for(i=0; i<m; i++){
        reversed <<= 1;
        reversed |= forward & 1;
        forward >>= 1;
    }
    return reversed;
}
 
/*calculate fft*/
void fft(float reals[], float imaginaries[]){
 
    int stage, subdft_size, span, i, ip, j;
    float sr, si, temp_real, temp_imaginary, imaginary_twiddle, real_twiddle;
 
    //read data into array
    for(i=0; i<n; i++){
        ip = bit_reverse(i);
        if(i < ip){
            temp_real = reals[i];
            temp_imaginary = imaginaries[i];
            reals[i] = reals[ip];
            imaginaries[i] = imaginaries[ip];
            reals[ip] = temp_real;
            imaginaries[ip] = temp_imaginary;
        }
    }
 
    //butterfly multiplies
    for(stage=1; stage<=m; stage++){
        report(stage);
        subdft_size = 1 << stage;
        span = subdft_size >> 1;
 
        //initialize trigonometric recurrence
 
        real_twiddle=1.0;
        imaginary_twiddle=0.0;
 
        sr = twiddle_step_real[stage];
        si = twiddle_step_imaginary[stage];
 
        for(j=0; j<span; j++){
            for(i=j; i<n; i+=subdft_size){
                ip=i+span;
 
                temp_real      = reals[ip]*real_twiddle      - imaginaries[ip]*imaginary_twiddle;
                temp_imaginary = reals[ip]*imaginary_twiddle + imaginaries[ip]*real_twiddle;
 
                reals[ip]       = reals[i]-temp_real;
                imaginaries[ip] = imaginaries[i]-temp_imaginary;
 
                reals[i]       = reals[i]+temp_real;
                imaginaries[i] = imaginaries[i]+temp_imaginary;
            }
            //trigonometric recreal_twiddlerence
            temp_real=real_twiddle;
            real_twiddle      = temp_real*sr - imaginary_twiddle*si;
            imaginary_twiddle = temp_real*si + imaginary_twiddle*sr;
        }
    }
}
 
void main(){
    float reals[n];
    float imaginaries[n];
    unsigned i;
 
    /* pre-calculate sine and cosine*/
    calculate_twiddles();
 
    /* generate a 64 sample cos wave */
    for(i=0; i<n; i++){
        reals[i] = 0.0;
        imaginaries[i] = 0.0;
    }
    for(i=0; i<=64; i++){
        reals[i] = sin(2.0 * M_PI * (i/64.0));
    }
 
    /* output time domain signal to a file */
    for(i=0; i<n; i++){
        file_write(reals[i], "x_re");
        file_write(imaginaries[i], "x_im");
    }
 
    /* transform into frequency domain */
    fft(reals, imaginaries);
 
    /* output frequency domain signal to a file */
    for(i=0; i<n; i++){
        file_write(reals[i], "fft_x_re");
        file_write(imaginaries[i], "fft_x_im");
    }
}