Root/usrp/fft.c

Source at commit fff1e1ed2f2c1654fa550107946c7858acd4484d created 8 years 10 months ago.
By Werner Almesberger, usrp/fft.c: added Hamming, make dump easier to interpret
1#include <stdlib.h>
2#include <stdio.h>
3#include <unistd.h>
4#include <math.h>
5#include <sys/types.h>
6
7#include <fftw3.h>
8
9
10#define DEFAULT_THRESHOLD 100
11
12static int alg = 0;
13
14
15static double window(int i, int n)
16{
17    return 0.54-0.46*cos(M_PI*2*i/(n-1));
18}
19
20
21static void fft_complex(int n, const float *re, const float *im, double *res)
22{
23    fftw_complex *in, *out;
24    fftw_plan plan;
25    int i;
26    double a;
27
28    in = fftw_malloc(sizeof(fftw_complex)*n);
29    out = fftw_malloc(sizeof(fftw_complex)*n);
30
31    for (i = 0; i != n; i++) {
32        double w = window(i, n);
33        in[i][0] = re[i]*w;
34        in[i][1] = im[i]*w;
35    }
36
37    plan = fftw_plan_dft_1d(n, in, out, FFTW_FORWARD, FFTW_ESTIMATE);
38    fftw_execute(plan);
39
40    for (i = 0; i != n; i++) {
41        a = hypot(out[i][0], out[i][1]); // /n;
42        a = a*a/n;
43        res[i] = a;
44    }
45}
46
47
48static void fft_real(int n, const float *re, const float *im, double *res)
49{
50    double *in;
51    fftw_plan plan;
52    int i;
53    double a ;
54
55    in = fftw_malloc(sizeof(double)*n);
56
57    for (i = 0; i != n; i++) {
58        a = hypot(re[i], im[i]);
59        in[i] = a*a;
60    }
61
62    plan = fftw_plan_r2r_1d(n, in, res, FFTW_REDFT10, FFTW_ESTIMATE);
63    fftw_execute(plan);
64
65    /* @@@ not sure at all about the scaling */
66    for (i = 0; i != n; i++)
67        res[i] = res[i]/sqrt(n);
68}
69
70
71static void do_fft(int skip, int dump, int low, int high, double threshold)
72{
73    float c[2];
74    float *re = NULL, *im = NULL;
75    double *res;
76    int e = 0, n = 0;
77    int i;
78    double a;
79
80    while (1) {
81        size_t s;
82
83        s = fread(c, sizeof(c), 1, stdin);
84        if (!s) {
85            if (!ferror(stdin))
86                break;
87            if (s < 0) {
88                perror("read");
89                exit(1);
90            }
91        }
92        if (skip) {
93            skip--;
94            continue;
95        }
96        if (e <= n) {
97            e = e ? e*2 : 10000;
98            re = realloc(re, e*sizeof(float));
99            im = realloc(im, e*sizeof(float));
100            if (!re || !im) {
101                perror("realloc");
102                exit(1);
103            }
104        }
105        re[n] = c[0];
106        im[n] = c[1];
107        n++;
108    }
109
110    if (skip >= n) {
111        fprintf(stderr, "cannot skip %d of %d entries\n", skip, n);
112        exit(1);
113    }
114    re += skip;
115    im += skip;
116    n -= skip;
117
118    res = malloc(n*sizeof(double));
119    if (!res) {
120        perror("malloc");
121        exit(1);
122    }
123
124    switch (alg) {
125    case 0:
126        fft_complex(n, re, im, res);
127        break;
128    case 1:
129        fft_real(n, re, im, res);
130        break;
131    default:
132        abort();
133    }
134
135    if (dump) {
136        for (i = 0; i < n; i += 1)
137            printf("%d %g\n", i,
138                10*log(res[(i+n/2) % n])/log(10));
139    } else {
140        double s = 0;
141        double db;
142
143        if (high >= n+skip) {
144            fprintf(stderr, "end %d > number of samples %d\n",
145                high, n+skip);
146            exit(1);
147        }
148        low = low*(double) n/(n+skip);
149        high = high*(double) n/(n+skip);
150        if (high < n)
151            high++;
152        if (low == high)
153            low--;
154        for (i = low; i != high; i++) {
155            a = res[i];
156            db = 10*log(a)/log(10);
157            if (db >= threshold)
158                s += a;
159        }
160        printf("%f\n", 10*log(s)/log(10));
161    }
162}
163
164
165static void usage(const char *name)
166{
167    fprintf(stderr,
168"usage: %s [-s skip] low high [threshold]\n"
169" %s [-s skip] -d\n\n"
170" threshold only use frequency bins with at least this power, in - dB.\n"
171" E.g., a threshold value of 60 would be -60 dB. (default: %d\n"
172" dB)\n"
173" -d dump frequency-domain \n"
174" -s skip skip this number of samples from the beginning (default: 0)\n"
175    , name, name, -DEFAULT_THRESHOLD);
176    exit(1);
177}
178
179
180int main(int argc, char **argv)
181{
182    int dump = 0, skip = 0;
183    int low, high;
184    double threshold = DEFAULT_THRESHOLD;
185    int c;
186
187    while ((c = getopt(argc, argv, "a:ds:")) != EOF)
188        switch (c) {
189        case 'a':
190            alg = atoi(optarg);
191            break;
192        case 'd':
193            dump = 1;
194            break;
195        case 's':
196            skip = atoi(optarg);
197            break;
198        default:
199            usage(*argv);
200        }
201
202    switch (argc-optind) {
203    case 0:
204        if (!dump)
205            usage(*argv);
206        do_fft(skip, 1, 0, 0, 0);
207        break;
208    case 3:
209        threshold = -atof(argv[optind+2]);
210        /* fall through */
211    case 2:
212        if (dump)
213            usage(*argv);
214        low = atoi(argv[optind]);
215        high = atoi(argv[optind+1]);
216        if (low > high)
217            usage(*argv);
218        do_fft(skip, 0, low, high, threshold);
219        break;
220    default:
221        usage(*argv);
222    }
223
224    return 0;
225}
226

Archive Download this file



interactive