diff options
Diffstat (limited to 'sem5/sig/mm13/fft.c')
| -rw-r--r-- | sem5/sig/mm13/fft.c | 182 |
1 files changed, 182 insertions, 0 deletions
diff --git a/sem5/sig/mm13/fft.c b/sem5/sig/mm13/fft.c new file mode 100644 index 0000000..dd93399 --- /dev/null +++ b/sem5/sig/mm13/fft.c @@ -0,0 +1,182 @@ +#include <stdio.h> +#include <stdlib.h> +#include <stdint.h> +#include <complex.h> +#include <math.h> + +// Needed for option parsing +#include <unistd.h> +#include <string.h> + +#define MAX_FILE_SIZE 10000 + +typedef uint16_t index_t; +typedef complex double val_t; + +// Make a better version of this +index_t reverse_bits(index_t i, unsigned fft_size) +{ + unsigned n = fft_size; + index_t reverse = 0; + + while (n--) { + reverse <<= 1; + reverse |= i & 1; + i >>= 1; + } + + return reverse; +} + +val_t twiddle(index_t N, int k) +{ + val_t in = (2 * M_PI)/(double)N; + val_t res = cexp(-in * k * I); + return res; +} + +void butterfly_single(val_t *buf, index_t a, index_t b, val_t twid) +{ + val_t common = twid * buf[b]; + + buf[b] = buf[a] - common; + buf[a] += common; +} + +void butterfly(val_t *buf, unsigned depth, index_t n, index_t N) +{ + if (n == 2) { + butterfly_single(buf, 0, 1, twiddle(N, 0)); + return; + } + + n = n/2; + index_t step = pow(2, depth); + + for (index_t i = 0; i < n; i++) { + butterfly_single(buf, i, i + n, twiddle(N, i * step)); + } +} + +void fft_inplace_recur(val_t *buf, unsigned depth, index_t n, index_t N) +{ + if (n > 2) { + index_t split = n/2; + fft_inplace_recur(buf, depth+1, split, N); + fft_inplace_recur(buf+split, depth+1, split, N); + } + + butterfly(buf, depth, n, N); +} + +val_t *fft(val_t *buf, size_t input_size, size_t N) +{ + // Check if power of two + if ((N & (N - 1)) != 0) { + fprintf(stderr, "FFT size %lu must be power of two\n", N); + return NULL; + } + + val_t *tmpbuf = (val_t *) malloc(sizeof(val_t) * N); + + // Reverse bits + int least = input_size < N ? input_size : N; + for (int i = 0; i < least; i++) { + index_t reverse = reverse_bits(i, log2(N)); + tmpbuf[i] = buf[reverse]; + } + + // Padding + for (int i = 0; i < (int)N - least; i++) { + tmpbuf[i + least ] = 0; + } + + fft_inplace_recur(tmpbuf, 0, N, N); + return tmpbuf; +} + +void printhelp(char *prog) { + printf("Usage: %s [options]\n", prog); + puts("Options:"); + puts(" -h This help screen"); + puts(" -o <FILE> Write output to file"); + puts(" -i <FILE> Read from file"); + puts(" -f <SIZE> FFT size, must be power of two"); + puts(""); + printf("FILE should not have more than %d lines, with line seperated values.\n", MAX_FILE_SIZE); + puts("FILE will default to stdin/stdout"); +} + +int main(int argc, char **argv) +{ + int opt; + int optfft = -1; + char *input = NULL; + char *output = NULL; + + while ((opt = getopt(argc, argv, "hi:o:f:")) != -1) { + switch(opt) { + case 'o': + output = strdup(optarg); + break; + case 'i': + input = strdup(optarg); + break; + case 'f': + optfft = strtol(optarg, NULL, 10); + break; + case 'h': + default: + printhelp(argv[0]); + return 0; + } + } + + if (optfft == -1) { + fprintf(stderr, "FFT not specified\n"); + printhelp(argv[0]); + return 1; + } + + // Allocate storage + val_t *buf = (val_t *) malloc(sizeof(val_t) * MAX_FILE_SIZE); + if (buf == NULL) { + printf("Could not allocate mega array\n"); + return 1; + } + + FILE *f = stdin; + if (input != NULL) { + f = fopen(input, "r"); + if (f == NULL) { + printf("Could not open file %s\n", input); + return 1; + } + } + + double val; + int length; + for (length = 0; fscanf(f, "%lf\n", &val) != EOF; length++) { + buf[length] = (val_t) val; + //printf("Val: %lf + j%lf\n", creal(buf[ii]), cimag(buf[ii])); + } + + val_t *res = fft(buf, length, optfft); + if (res == NULL) { + fprintf(stderr, "Error calculating fft\n"); + return 1; + } + + f = stdout; + if (output != NULL) { + f = fopen(output, "w"); + if (f == NULL) { + printf("Could not open file %s\n", input); + return 1; + } + } + + for (int ii = 0; ii < optfft; ii++) { + fprintf(f, "%lf%+lfj\n", creal(res[ii]), cimag(res[ii])); + } +} |