The Code-Bin
Links
Home
Add your code!
All Listings
About
Latest Entry
Featured Scripts
Author's Website
Latest Entries
FFMPEG Thumbnail Scr...
PHP, 0.8KB
Jul. 29, 10:24pm
John
Z80 Assembler, 190 bytes
Feb. 17, 3:36am
John
Z80 Assembler, 176 bytes
Sep. 13, 2:19am
John
Z80 Assembler, 77 bytes
Sep. 13, 2:18am
John
Z80 Assembler, 209 bytes
Sep. 13, 2:17am
untitled C Code
Posted by: M | September 19, 2010 @ 9:58pm
C Code
[
Download
]
#include <stdio.h> #include <stdlib.h> #include <math.h> #define logb(v, b) (log(v)/log(b)) #define SAMPLES 8 #define PI 3.1415926535897932 // complex datatype typedef struct complex_t { double re; // Real part double im; // Imag part } complex; // function headers complex c_from_polar(double r, double radians); complex c_add(complex left, complex right); complex c_sub(complex left, complex right); complex c_mult(complex left, complex right); double c_magnitude(complex c); int* twiddleIndexGen(int* twindex, int log2N); complex* twiddleGen(complex* twiddles, int* twindex, int log2N); complex* DFT_naive(complex* x, int N); // definitions // convert polar value to complex complex c_from_polar(double r, double radians){ complex result; result.re = r * cos(radians); result.im = r * sin(radians); return result; } // add two complex numbers complex c_add(complex left, complex right){ complex result; result.re = left.re + right.re; result.im = left.re + right.re; return result; } // subtract two complex numbers complex c_sub(complex left, complex right){ complex result; result.re = left.re - right.re; result.im = left.re - right.re; return result; } // multiply two complex numbers complex c_mult(complex left, complex right){ complex result; result.re = left.re * right.re - left.im * right.im; result.im = left.re * right.im - left.im * right.re; return result; } // compute magnitude of complex number double c_magnitude(complex c){ return sqrt(c.re * c.re + c.im * c.im); } // compute DFT for input x complex* DFT_naive(complex* x, int N){ complex* X = (complex*) malloc(sizeof(struct complex_t) * N); complex* w = (complex*) malloc(sizeof(struct complex_t) * N); int k, n; for (k = 0; k < N; k++){ w[k] = c_from_polar(1, -2*PI*k/N); } for (k = 0; k < N; k++){ X[k].re = X[k].im = 0.0; for (n = 0; n < N; n++){ //X[k] = c_add(X[k[, c_mult(x[n], w[(n*k)%N])); } } free(w); return X; } /* twiddleIndexGen: * Precalculates the final index order * of the twiddle factors based on * the number of samples. * Uses the bit reversal method */ int* twiddleIndexGen(int* twindex, int log2N){ int i, j; for (j = 0; j < (1<<log2N); j++){ // Index value int this_num = j; for (i = 0; i < log2N; i++){ // Bit twindex[j] <<= 1; twindex[j] += this_num & 1; // is this bit 1? this_num >>= 1; } printf("%d ", twindex[j]); } printf("\n"); return twindex; } complex* twiddleGen(complex* twiddles, int* twindex, int log2N){ twindex = twiddleIndexGen(twindex, log2N); int i, j; int N = (1<<log2N); for (i = 0; i < N; i++){ j = twindex[i]; twiddles[j] = c_from_polar(1, -2*PI*j/N); printf("%f, %f\n", twiddles[j].re, twiddles[j].im); } return twiddles; } /* RFFT: * Recursive algorithm. Traverses through the array * of sample values in a predefined order before * traversing back and computing FFTs using twiddle * factors. */ complex* RFFT(complex* top, complex* out, complex** left, complex** right, int l){ int i; printf("\nL= %i\n", l); for (i = 0; i < SAMPLES/(1<<l); i++){ left[l][i] = top[2*i]; // Fill LHS printf("%.0f ", left[l][i].re); right[l][i] = top[2*i+1]; // Fill RHS printf("%.0f\n", right[l][i].re); } if (1<<l>=SAMPLES){ } else { RFFT(left[l], out, left, right, l+1); // Traverse left RFFT(right[l], out, left, right, l+1); // Traverse right } return top; } int main(int argc, char **argv) { int i; int ls = logb(SAMPLES, 2); // # levels complex *main = (complex *)malloc(sizeof(complex)*SAMPLES); // main array complex *out = (complex *)malloc(sizeof(complex)*SAMPLES); // output array complex **left = (complex **)malloc(sizeof(complex)*SAMPLES); // array of left traverse arrays complex **right = (complex **)malloc(sizeof(complex)*SAMPLES); // array of right traverse arrays /* (only need one for each per level since only working with a single pair at one time) * each level down the tree requires only half the number of values as the previous * so array size is halved. Total array size is roughly 4x sample size */ // Populate MAIN array for (i = 0; i < SAMPLES; i++) { main[i] = (complex){i,0}; //printf("%.0f,%.0f\n\n", main[i].re, main[i].im); } // Generate twiddle indexes and twiddle factors int *twindex = (int*)malloc(sizeof(int)*SAMPLES); complex *twiddles = (complex*)malloc(sizeof(complex)*SAMPLES); twiddles = twiddleGen(twiddles, twindex, ls); // Define LEFT/RIGHT array sizes, each is 1/2 // as long as the previous array for (i = ls; i >= 0; i--){ left[ls-i+1] = (complex*)malloc(sizeof(complex)*(1<<i)); right[ls-i+1] = (complex*)malloc(sizeof(complex)*(1<<i)); } //Begin Recursion int level = 1; RFFT(main, out, left, right, level); // Cleanup for (i = ls; i >= 0; i--){ free(left[ls-i+1]); free(right[ls-i+1]); } free(left); free(right); free(twiddles); free(twindex); free(out); free(main); return 0; }
Syntax Highlighting
[
Open in new window
]
Author Comments
none
Rating
4.47 / 8
53 Votes
http://codebin.yi.org/845
page generated in 0.00 seconds