1 //
2 // fft_example.c
3 //
4 // This example demonstrates the interface to the fast discrete Fourier
5 // transform (FFT).
6 // SEE ALSO: mdct_example.c
7 // fct_example.c
8 //
9
10 #include <stdio.h>
11 #include <stdlib.h>
12 #include <math.h>
13 #include <getopt.h>
14 #include "liquid.h"
15
16 // print usage/help message
usage()17 void usage()
18 {
19 printf("fft_example [options]\n");
20 printf(" h : print help\n");
21 printf(" v/q : verbose/quiet\n");
22 printf(" n : fft size, default: 16\n");
23 }
24
main(int argc,char * argv[])25 int main(int argc, char*argv[])
26 {
27 // options
28 unsigned int nfft = 16; // transform size
29 int method = 0; // fft method (ignored)
30 int verbose = 0; // verbose output?
31
32 int dopt;
33 while ((dopt = getopt(argc,argv,"hvqn:")) != EOF) {
34 switch (dopt) {
35 case 'h': usage(); return 0;
36 case 'v': verbose = 1; break;
37 case 'q': verbose = 0; break;
38 case 'n': nfft = atoi(optarg); break;
39 default:
40 exit(1);
41 }
42 }
43
44 // allocate memory arrays
45 float complex * x = (float complex*) malloc(nfft*sizeof(float complex));
46 float complex * y = (float complex*) malloc(nfft*sizeof(float complex));
47 float complex * z = (float complex*) malloc(nfft*sizeof(float complex));
48
49 // initialize input
50 unsigned int i;
51 for (i=0; i<nfft; i++)
52 x[i] = (float)i - _Complex_I*(float)i;
53
54 // create fft plans
55 fftplan pf = fft_create_plan(nfft, x, y, LIQUID_FFT_FORWARD, method);
56 fftplan pr = fft_create_plan(nfft, y, z, LIQUID_FFT_BACKWARD, method);
57
58 // print fft plans
59 fft_print_plan(pf);
60 //fft_print_plan(pr);
61
62 // execute fft plans
63 fft_execute(pf);
64 fft_execute(pr);
65
66 // destroy fft plans
67 fft_destroy_plan(pf);
68 fft_destroy_plan(pr);
69
70 // normalize inverse
71 for (i=0; i<nfft; i++)
72 z[i] /= (float) nfft;
73
74 if (verbose) {
75 // print results
76 printf("original signal, x[n]:\n");
77 for (i=0; i<nfft; i++)
78 printf(" x[%3u] = %8.4f + j%8.4f\n", i, crealf(x[i]), cimagf(x[i]));
79 printf("y[n] = fft( x[n] ):\n");
80 for (i=0; i<nfft; i++)
81 printf(" y[%3u] = %8.4f + j%8.4f\n", i, crealf(y[i]), cimagf(y[i]));
82 printf("z[n] = ifft( y[n] ):\n");
83 for (i=0; i<nfft; i++)
84 printf(" z[%3u] = %8.4f + j%8.4f\n", i, crealf(z[i]), cimagf(z[i]));
85 }
86
87 // compute RMSE between original and result
88 float rmse = 0.0f;
89 for (i=0; i<nfft; i++) {
90 float complex d = x[i] - z[i];
91 rmse += crealf(d * conjf(d));
92 }
93 rmse = sqrtf( rmse / (float)nfft );
94 printf("rmse = %12.4e\n", rmse);
95
96 // free allocated memory
97 free(x);
98 free(y);
99 free(z);
100
101 printf("done.\n");
102 return 0;
103 }
104
105