1 /*
2 * rtl-sdr, turns your Realtek RTL2832 based DVB dongle into a SDR receiver
3 * Copyright (C) 2012 by Steve Markgraf <steve@steve-m.de>
4 * Copyright (C) 2012 by Hoernchen <la@tfc-server.de>
5 * Copyright (C) 2012 by Kyle Keen <keenerd@gmail.com>
6 *
7 * This program is free software: you can redistribute it and/or modify
8 * it under the terms of the GNU General Public License as published by
9 * the Free Software Foundation, either version 2 of the License, or
10 * (at your option) any later version.
11 *
12 * This program is distributed in the hope that it will be useful,
13 * but WITHOUT ANY WARRANTY; without even the implied warranty of
14 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
15 * GNU General Public License for more details.
16 *
17 * You should have received a copy of the GNU General Public License
18 * along with this program. If not, see <http://www.gnu.org/licenses/>.
19 */
20
21
22 /*
23 * rtl_power: general purpose FFT integrator
24 * -f low_freq:high_freq:max_bin_size
25 * -i seconds
26 * outputs CSV
27 * time, low, high, step, db, db, db ...
28 * db optional? raw output might be better for noise correction
29 * todo:
30 * threading
31 * randomized hopping
32 * noise correction
33 * continuous IIR
34 * general astronomy usefulness
35 * multiple dongles
36 * multiple FFT workers
37 * check edge cropping for off-by-one and rounding errors
38 * 1.8MS/s for hiding xtal harmonics
39 */
40
41 #include <errno.h>
42 #include <signal.h>
43 #include <string.h>
44 #include <stdio.h>
45 #include <stdlib.h>
46 #include <time.h>
47
48 #ifndef _WIN32
49 #include <unistd.h>
50 #else
51 #include <windows.h>
52 #include <fcntl.h>
53 #include <io.h>
54 #include "getopt/getopt.h"
55 #define usleep(x) Sleep(x/1000)
56 #if defined(_MSC_VER) && (_MSC_VER < 1800)
57 #define round(x) (x > 0.0 ? floor(x + 0.5): ceil(x - 0.5))
58 #endif
59 #define _USE_MATH_DEFINES
60 #endif
61
62 #include <math.h>
63 #include <pthread.h>
64 #include <libusb.h>
65
66 #include "rtl-sdr.h"
67 #include "convenience/convenience.h"
68
69 #define MAX(x, y) (((x) > (y)) ? (x) : (y))
70
71 #define DEFAULT_BUF_LENGTH (1 * 16384)
72 #define AUTO_GAIN -100
73 #define BUFFER_DUMP (1<<12)
74
75 #define MAXIMUM_RATE 2800000
76 #define MINIMUM_RATE 1000000
77
78 static volatile int do_exit = 0;
79 static rtlsdr_dev_t *dev = NULL;
80 FILE *file;
81
82 int16_t* Sinewave;
83 double* power_table;
84 int N_WAVE, LOG2_N_WAVE;
85 int next_power;
86 int16_t *fft_buf;
87 int *window_coefs;
88
89 struct tuning_state
90 /* one per tuning range */
91 {
92 int freq;
93 int rate;
94 int bin_e;
95 long *avg; /* length == 2^bin_e */
96 int samples;
97 int downsample;
98 int downsample_passes; /* for the recursive filter */
99 double crop;
100 //pthread_rwlock_t avg_lock;
101 //pthread_mutex_t avg_mutex;
102 /* having the iq buffer here is wasteful, but will avoid contention */
103 uint8_t *buf8;
104 int buf_len;
105 //int *comp_fir;
106 //pthread_rwlock_t buf_lock;
107 //pthread_mutex_t buf_mutex;
108 };
109
110 /* 3000 is enough for 3GHz b/w worst case */
111 #define MAX_TUNES 3000
112 struct tuning_state tunes[MAX_TUNES];
113 int tune_count = 0;
114
115 int boxcar = 1;
116 int comp_fir_size = 0;
117 int peak_hold = 0;
118
usage(void)119 void usage(void)
120 {
121 fprintf(stderr,
122 "rtl_power, a simple FFT logger for RTL2832 based DVB-T receivers\n\n"
123 "Use:\trtl_power -f freq_range [-options] [filename]\n"
124 "\t-f lower:upper:bin_size [Hz]\n"
125 "\t (bin size is a maximum, smaller more convenient bins\n"
126 "\t will be used. valid range 1Hz - 2.8MHz)\n"
127 "\t[-i integration_interval (default: 10 seconds)]\n"
128 "\t (buggy if a full sweep takes longer than the interval)\n"
129 "\t[-1 enables single-shot mode (default: off)]\n"
130 "\t[-e exit_timer (default: off/0)]\n"
131 //"\t[-s avg/iir smoothing (default: avg)]\n"
132 //"\t[-t threads (default: 1)]\n"
133 "\t[-d device_index (default: 0)]\n"
134 "\t[-g tuner_gain (default: automatic)]\n"
135 "\t[-p ppm_error (default: 0)]\n"
136 "\t[-T enable bias-T on GPIO PIN 0 (works for rtl-sdr.com v3 dongles)]\n"
137 "\tfilename (a '-' dumps samples to stdout)\n"
138 "\t (omitting the filename also uses stdout)\n"
139 "\n"
140 "Experimental options:\n"
141 "\t[-w window (default: rectangle)]\n"
142 "\t (hamming, blackman, blackman-harris, hann-poisson, bartlett, youssef)\n"
143 // kaiser
144 "\t[-c crop_percent (default: 0%%, recommended: 20%%-50%%)]\n"
145 "\t (discards data at the edges, 100%% discards everything)\n"
146 "\t (has no effect for bins larger than 1MHz)\n"
147 "\t[-F fir_size (default: disabled)]\n"
148 "\t (enables low-leakage downsample filter,\n"
149 "\t fir_size can be 0 or 9. 0 has bad roll off,\n"
150 "\t try with '-c 50%%')\n"
151 "\t[-P enables peak hold (default: off)]\n"
152 "\t[-D enable direct sampling (default: off)]\n"
153 "\t[-O enable offset tuning (default: off)]\n"
154 "\n"
155 "CSV FFT output columns:\n"
156 "\tdate, time, Hz low, Hz high, Hz step, samples, dbm, dbm, ...\n\n"
157 "Examples:\n"
158 "\trtl_power -f 88M:108M:125k fm_stations.csv\n"
159 "\t (creates 160 bins across the FM band,\n"
160 "\t individual stations should be visible)\n"
161 "\trtl_power -f 100M:1G:1M -i 5m -1 survey.csv\n"
162 "\t (a five minute low res scan of nearly everything)\n"
163 "\trtl_power -f ... -i 15m -1 log.csv\n"
164 "\t (integrate for 15 minutes and exit afterwards)\n"
165 "\trtl_power -f ... -e 1h | gzip > log.csv.gz\n"
166 "\t (collect data for one hour and compress it on the fly)\n\n"
167 "Convert CSV to a waterfall graphic with:\n"
168 "\t http://kmkeen.com/tmp/heatmap.py.txt \n");
169 exit(1);
170 }
171
multi_bail(void)172 void multi_bail(void)
173 {
174 if (do_exit == 1)
175 {
176 fprintf(stderr, "Signal caught, finishing scan pass.\n");
177 }
178 if (do_exit >= 2)
179 {
180 fprintf(stderr, "Signal caught, aborting immediately.\n");
181 }
182 }
183
184 #ifdef _WIN32
185 BOOL WINAPI
sighandler(int signum)186 sighandler(int signum)
187 {
188 if (CTRL_C_EVENT == signum) {
189 do_exit++;
190 multi_bail();
191 return TRUE;
192 }
193 return FALSE;
194 }
195 #else
sighandler(int signum)196 static void sighandler(int signum)
197 {
198 do_exit++;
199 multi_bail();
200 }
201 #endif
202
203 /* more cond dumbness */
204 #define safe_cond_signal(n, m) pthread_mutex_lock(m); pthread_cond_signal(n); pthread_mutex_unlock(m)
205 #define safe_cond_wait(n, m) pthread_mutex_lock(m); pthread_cond_wait(n, m); pthread_mutex_unlock(m)
206
207 /* {length, coef, coef, coef} and scaled by 2^15
208 for now, only length 9, optimal way to get +85% bandwidth */
209 #define CIC_TABLE_MAX 10
210 int cic_9_tables[][10] = {
211 {0,},
212 {9, -156, -97, 2798, -15489, 61019, -15489, 2798, -97, -156},
213 {9, -128, -568, 5593, -24125, 74126, -24125, 5593, -568, -128},
214 {9, -129, -639, 6187, -26281, 77511, -26281, 6187, -639, -129},
215 {9, -122, -612, 6082, -26353, 77818, -26353, 6082, -612, -122},
216 {9, -120, -602, 6015, -26269, 77757, -26269, 6015, -602, -120},
217 {9, -120, -582, 5951, -26128, 77542, -26128, 5951, -582, -120},
218 {9, -119, -580, 5931, -26094, 77505, -26094, 5931, -580, -119},
219 {9, -119, -578, 5921, -26077, 77484, -26077, 5921, -578, -119},
220 {9, -119, -577, 5917, -26067, 77473, -26067, 5917, -577, -119},
221 {9, -199, -362, 5303, -25505, 77489, -25505, 5303, -362, -199},
222 };
223
224 #if defined(_MSC_VER) && (_MSC_VER < 1800)
log2(double n)225 double log2(double n)
226 {
227 return log(n) / log(2.0);
228 }
229 #endif
230
231 /* FFT based on fix_fft.c by Roberts, Slaney and Bouras
232 http://www.jjj.de/fft/fftpage.html
233 16 bit ints for everything
234 -32768..+32768 maps to -1.0..+1.0
235 */
236
sine_table(int size)237 void sine_table(int size)
238 {
239 int i;
240 double d;
241 LOG2_N_WAVE = size;
242 N_WAVE = 1 << LOG2_N_WAVE;
243 Sinewave = malloc(sizeof(int16_t) * N_WAVE*3/4);
244 power_table = malloc(sizeof(double) * N_WAVE);
245 for (i=0; i<N_WAVE*3/4; i++)
246 {
247 d = (double)i * 2.0 * M_PI / N_WAVE;
248 Sinewave[i] = (int)round(32767*sin(d));
249 //printf("%i\n", Sinewave[i]);
250 }
251 }
252
FIX_MPY(int16_t a,int16_t b)253 static inline int16_t FIX_MPY(int16_t a, int16_t b)
254 /* fixed point multiply and scale */
255 {
256 int c = ((int)a * (int)b) >> 14;
257 b = c & 0x01;
258 return (c >> 1) + b;
259 }
260
fix_fft(int16_t iq[],int m)261 int fix_fft(int16_t iq[], int m)
262 /* interleaved iq[], 0 <= n < 2**m, changes in place */
263 {
264 int mr, nn, i, j, l, k, istep, n, shift;
265 int16_t qr, qi, tr, ti, wr, wi;
266 n = 1 << m;
267 if (n > N_WAVE)
268 {return -1;}
269 mr = 0;
270 nn = n - 1;
271 /* decimation in time - re-order data */
272 for (m=1; m<=nn; ++m) {
273 l = n;
274 do
275 {l >>= 1;}
276 while (mr+l > nn);
277 mr = (mr & (l-1)) + l;
278 if (mr <= m)
279 {continue;}
280 // real = 2*m, imag = 2*m+1
281 tr = iq[2*m];
282 iq[2*m] = iq[2*mr];
283 iq[2*mr] = tr;
284 ti = iq[2*m+1];
285 iq[2*m+1] = iq[2*mr+1];
286 iq[2*mr+1] = ti;
287 }
288 l = 1;
289 k = LOG2_N_WAVE-1;
290 while (l < n) {
291 shift = 1;
292 istep = l << 1;
293 for (m=0; m<l; ++m) {
294 j = m << k;
295 wr = Sinewave[j+N_WAVE/4];
296 wi = -Sinewave[j];
297 if (shift) {
298 wr >>= 1; wi >>= 1;}
299 for (i=m; i<n; i+=istep) {
300 j = i + l;
301 tr = FIX_MPY(wr,iq[2*j]) - FIX_MPY(wi,iq[2*j+1]);
302 ti = FIX_MPY(wr,iq[2*j+1]) + FIX_MPY(wi,iq[2*j]);
303 qr = iq[2*i];
304 qi = iq[2*i+1];
305 if (shift) {
306 qr >>= 1; qi >>= 1;}
307 iq[2*j] = qr - tr;
308 iq[2*j+1] = qi - ti;
309 iq[2*i] = qr + tr;
310 iq[2*i+1] = qi + ti;
311 }
312 }
313 --k;
314 l = istep;
315 }
316 return 0;
317 }
318
rectangle(int i,int length)319 double rectangle(int i, int length)
320 {
321 return 1.0;
322 }
323
hamming(int i,int length)324 double hamming(int i, int length)
325 {
326 double a, b, w, N1;
327 a = 25.0/46.0;
328 b = 21.0/46.0;
329 N1 = (double)(length-1);
330 w = a - b*cos(2*i*M_PI/N1);
331 return w;
332 }
333
blackman(int i,int length)334 double blackman(int i, int length)
335 {
336 double a0, a1, a2, w, N1;
337 a0 = 7938.0/18608.0;
338 a1 = 9240.0/18608.0;
339 a2 = 1430.0/18608.0;
340 N1 = (double)(length-1);
341 w = a0 - a1*cos(2*i*M_PI/N1) + a2*cos(4*i*M_PI/N1);
342 return w;
343 }
344
blackman_harris(int i,int length)345 double blackman_harris(int i, int length)
346 {
347 double a0, a1, a2, a3, w, N1;
348 a0 = 0.35875;
349 a1 = 0.48829;
350 a2 = 0.14128;
351 a3 = 0.01168;
352 N1 = (double)(length-1);
353 w = a0 - a1*cos(2*i*M_PI/N1) + a2*cos(4*i*M_PI/N1) - a3*cos(6*i*M_PI/N1);
354 return w;
355 }
356
hann_poisson(int i,int length)357 double hann_poisson(int i, int length)
358 {
359 double a, N1, w;
360 a = 2.0;
361 N1 = (double)(length-1);
362 w = 0.5 * (1 - cos(2*M_PI*i/N1)) * \
363 pow(M_E, (-a*(double)abs((int)(N1-1-2*i)))/N1);
364 return w;
365 }
366
youssef(int i,int length)367 double youssef(int i, int length)
368 /* really a blackman-harris-poisson window, but that is a mouthful */
369 {
370 double a, a0, a1, a2, a3, w, N1;
371 a0 = 0.35875;
372 a1 = 0.48829;
373 a2 = 0.14128;
374 a3 = 0.01168;
375 N1 = (double)(length-1);
376 w = a0 - a1*cos(2*i*M_PI/N1) + a2*cos(4*i*M_PI/N1) - a3*cos(6*i*M_PI/N1);
377 a = 0.0025;
378 w *= pow(M_E, (-a*(double)abs((int)(N1-1-2*i)))/N1);
379 return w;
380 }
381
kaiser(int i,int length)382 double kaiser(int i, int length)
383 // todo, become more smart
384 {
385 return 1.0;
386 }
387
bartlett(int i,int length)388 double bartlett(int i, int length)
389 {
390 double N1, L, w;
391 L = (double)length;
392 N1 = L - 1;
393 w = (i - N1/2) / (L/2);
394 if (w < 0) {
395 w = -w;}
396 w = 1 - w;
397 return w;
398 }
399
rms_power(struct tuning_state * ts)400 void rms_power(struct tuning_state *ts)
401 /* for bins between 1MHz and 2MHz */
402 {
403 int i, s;
404 uint8_t *buf = ts->buf8;
405 int buf_len = ts->buf_len;
406 long p, t;
407 double dc, err;
408
409 p = t = 0L;
410 for (i=0; i<buf_len; i++) {
411 s = (int)buf[i] - 127;
412 t += (long)s;
413 p += (long)(s * s);
414 }
415 /* correct for dc offset in squares */
416 dc = (double)t / (double)buf_len;
417 err = t * 2 * dc - dc * dc * buf_len;
418 p -= (long)round(err);
419
420 if (!peak_hold) {
421 ts->avg[0] += p;
422 } else {
423 ts->avg[0] = MAX(ts->avg[0], p);
424 }
425 ts->samples += 1;
426 }
427
frequency_range(char * arg,double crop)428 void frequency_range(char *arg, double crop)
429 /* flesh out the tunes[] for scanning */
430 // do we want the fewest ranges (easy) or the fewest bins (harder)?
431 {
432 char *start, *stop, *step;
433 int i, j, upper, lower, max_size, bw_seen, bw_used, bin_e, buf_len;
434 int downsample, downsample_passes;
435 double bin_size;
436 struct tuning_state *ts;
437 /* hacky string parsing */
438 start = arg;
439 stop = strchr(start, ':') + 1;
440 stop[-1] = '\0';
441 step = strchr(stop, ':') + 1;
442 step[-1] = '\0';
443 lower = (int)atofs(start);
444 upper = (int)atofs(stop);
445 max_size = (int)atofs(step);
446 stop[-1] = ':';
447 step[-1] = ':';
448 downsample = 1;
449 downsample_passes = 0;
450 /* evenly sized ranges, as close to MAXIMUM_RATE as possible */
451 // todo, replace loop with algebra
452 for (i=1; i<1500; i++) {
453 bw_seen = (upper - lower) / i;
454 bw_used = (int)((double)(bw_seen) / (1.0 - crop));
455 if (bw_used > MAXIMUM_RATE) {
456 continue;}
457 tune_count = i;
458 break;
459 }
460 /* unless small bandwidth */
461 if (bw_used < MINIMUM_RATE) {
462 tune_count = 1;
463 downsample = MAXIMUM_RATE / bw_used;
464 bw_used = bw_used * downsample;
465 }
466 if (!boxcar && downsample > 1) {
467 downsample_passes = (int)log2(downsample);
468 downsample = 1 << downsample_passes;
469 bw_used = (int)((double)(bw_seen * downsample) / (1.0 - crop));
470 }
471 /* number of bins is power-of-two, bin size is under limit */
472 // todo, replace loop with log2
473 for (i=1; i<=21; i++) {
474 bin_e = i;
475 bin_size = (double)bw_used / (double)((1<<i) * downsample);
476 if (bin_size <= (double)max_size) {
477 break;}
478 }
479 /* unless giant bins */
480 if (max_size >= MINIMUM_RATE) {
481 bw_seen = max_size;
482 bw_used = max_size;
483 tune_count = (upper - lower) / bw_seen;
484 bin_e = 0;
485 crop = 0;
486 }
487 if (tune_count > MAX_TUNES) {
488 fprintf(stderr, "Error: bandwidth too wide.\n");
489 exit(1);
490 }
491 buf_len = 2 * (1<<bin_e) * downsample;
492 if (buf_len < DEFAULT_BUF_LENGTH) {
493 buf_len = DEFAULT_BUF_LENGTH;
494 }
495 /* build the array */
496 for (i=0; i<tune_count; i++) {
497 ts = &tunes[i];
498 ts->freq = lower + i*bw_seen + bw_seen/2;
499 ts->rate = bw_used;
500 ts->bin_e = bin_e;
501 ts->samples = 0;
502 ts->crop = crop;
503 ts->downsample = downsample;
504 ts->downsample_passes = downsample_passes;
505 ts->avg = (long*)malloc((1<<bin_e) * sizeof(long));
506 if (!ts->avg) {
507 fprintf(stderr, "Error: malloc.\n");
508 exit(1);
509 }
510 for (j=0; j<(1<<bin_e); j++) {
511 ts->avg[j] = 0L;
512 }
513 ts->buf8 = (uint8_t*)malloc(buf_len * sizeof(uint8_t));
514 if (!ts->buf8) {
515 fprintf(stderr, "Error: malloc.\n");
516 exit(1);
517 }
518 ts->buf_len = buf_len;
519 }
520 /* report */
521 fprintf(stderr, "Number of frequency hops: %i\n", tune_count);
522 fprintf(stderr, "Dongle bandwidth: %iHz\n", bw_used);
523 fprintf(stderr, "Downsampling by: %ix\n", downsample);
524 fprintf(stderr, "Cropping by: %0.2f%%\n", crop*100);
525 fprintf(stderr, "Total FFT bins: %i\n", tune_count * (1<<bin_e));
526 fprintf(stderr, "Logged FFT bins: %i\n", \
527 (int)((double)(tune_count * (1<<bin_e)) * (1.0-crop)));
528 fprintf(stderr, "FFT bin size: %0.2fHz\n", bin_size);
529 fprintf(stderr, "Buffer size: %i bytes (%0.2fms)\n", buf_len, 1000 * 0.5 * (float)buf_len / (float)bw_used);
530 }
531
retune(rtlsdr_dev_t * d,int freq)532 void retune(rtlsdr_dev_t *d, int freq)
533 {
534 uint8_t dump[BUFFER_DUMP];
535 int n_read;
536 rtlsdr_set_center_freq(d, (uint32_t)freq);
537 /* wait for settling and flush buffer */
538 usleep(5000);
539 rtlsdr_read_sync(d, &dump, BUFFER_DUMP, &n_read);
540 if (n_read != BUFFER_DUMP) {
541 fprintf(stderr, "Error: bad retune.\n");}
542 }
543
fifth_order(int16_t * data,int length)544 void fifth_order(int16_t *data, int length)
545 /* for half of interleaved data */
546 {
547 int i;
548 int a, b, c, d, e, f;
549 a = data[0];
550 b = data[2];
551 c = data[4];
552 d = data[6];
553 e = data[8];
554 f = data[10];
555 /* a downsample should improve resolution, so don't fully shift */
556 /* ease in instead of being stateful */
557 data[0] = ((a+b)*10 + (c+d)*5 + d + f) >> 4;
558 data[2] = ((b+c)*10 + (a+d)*5 + e + f) >> 4;
559 data[4] = (a + (b+e)*5 + (c+d)*10 + f) >> 4;
560 for (i=12; i<length; i+=4) {
561 a = c;
562 b = d;
563 c = e;
564 d = f;
565 e = data[i-2];
566 f = data[i];
567 data[i/2] = (a + (b+e)*5 + (c+d)*10 + f) >> 4;
568 }
569 }
570
remove_dc(int16_t * data,int length)571 void remove_dc(int16_t *data, int length)
572 /* works on interleaved data */
573 {
574 int i;
575 int16_t ave;
576 long sum = 0L;
577 for (i=0; i < length; i+=2) {
578 sum += data[i];
579 }
580 ave = (int16_t)(sum / (long)(length));
581 if (ave == 0) {
582 return;}
583 for (i=0; i < length; i+=2) {
584 data[i] -= ave;
585 }
586 }
587
generic_fir(int16_t * data,int length,int * fir)588 void generic_fir(int16_t *data, int length, int *fir)
589 /* Okay, not at all generic. Assumes length 9, fix that eventually. */
590 {
591 int d, temp, sum;
592 int hist[9] = {0,};
593 /* cheat on the beginning, let it go unfiltered */
594 for (d=0; d<18; d+=2) {
595 hist[d/2] = data[d];
596 }
597 for (d=18; d<length; d+=2) {
598 temp = data[d];
599 sum = 0;
600 sum += (hist[0] + hist[8]) * fir[1];
601 sum += (hist[1] + hist[7]) * fir[2];
602 sum += (hist[2] + hist[6]) * fir[3];
603 sum += (hist[3] + hist[5]) * fir[4];
604 sum += hist[4] * fir[5];
605 data[d] = (int16_t)(sum >> 15) ;
606 hist[0] = hist[1];
607 hist[1] = hist[2];
608 hist[2] = hist[3];
609 hist[3] = hist[4];
610 hist[4] = hist[5];
611 hist[5] = hist[6];
612 hist[6] = hist[7];
613 hist[7] = hist[8];
614 hist[8] = temp;
615 }
616 }
617
downsample_iq(int16_t * data,int length)618 void downsample_iq(int16_t *data, int length)
619 {
620 fifth_order(data, length);
621 //remove_dc(data, length);
622 fifth_order(data+1, length-1);
623 //remove_dc(data+1, length-1);
624 }
625
real_conj(int16_t real,int16_t imag)626 long real_conj(int16_t real, int16_t imag)
627 /* real(n * conj(n)) */
628 {
629 return ((long)real*(long)real + (long)imag*(long)imag);
630 }
631
scanner(void)632 void scanner(void)
633 {
634 int i, j, j2, f, n_read, offset, bin_e, bin_len, buf_len, ds, ds_p;
635 int32_t w;
636 struct tuning_state *ts;
637 bin_e = tunes[0].bin_e;
638 bin_len = 1 << bin_e;
639 buf_len = tunes[0].buf_len;
640 for (i=0; i<tune_count; i++) {
641 if (do_exit >= 2)
642 {return;}
643 ts = &tunes[i];
644 f = (int)rtlsdr_get_center_freq(dev);
645 if (f != ts->freq) {
646 retune(dev, ts->freq);}
647 rtlsdr_read_sync(dev, ts->buf8, buf_len, &n_read);
648 if (n_read != buf_len) {
649 fprintf(stderr, "Error: dropped samples.\n");}
650 /* rms */
651 if (bin_len == 1) {
652 rms_power(ts);
653 continue;
654 }
655 /* prep for fft */
656 for (j=0; j<buf_len; j++) {
657 fft_buf[j] = (int16_t)ts->buf8[j] - 127;
658 }
659 ds = ts->downsample;
660 ds_p = ts->downsample_passes;
661 if (boxcar && ds > 1) {
662 j=2, j2=0;
663 while (j < buf_len) {
664 fft_buf[j2] += fft_buf[j];
665 fft_buf[j2+1] += fft_buf[j+1];
666 fft_buf[j] = 0;
667 fft_buf[j+1] = 0;
668 j += 2;
669 if (j % (ds*2) == 0) {
670 j2 += 2;}
671 }
672 } else if (ds_p) { /* recursive */
673 for (j=0; j < ds_p; j++) {
674 downsample_iq(fft_buf, buf_len >> j);
675 }
676 /* droop compensation */
677 if (comp_fir_size == 9 && ds_p <= CIC_TABLE_MAX) {
678 generic_fir(fft_buf, buf_len >> j, cic_9_tables[ds_p]);
679 generic_fir(fft_buf+1, (buf_len >> j)-1, cic_9_tables[ds_p]);
680 }
681 }
682 remove_dc(fft_buf, buf_len / ds);
683 remove_dc(fft_buf+1, (buf_len / ds) - 1);
684 /* window function and fft */
685 for (offset=0; offset<(buf_len/ds); offset+=(2*bin_len)) {
686 // todo, let rect skip this
687 for (j=0; j<bin_len; j++) {
688 w = (int32_t)fft_buf[offset+j*2];
689 w *= (int32_t)(window_coefs[j]);
690 //w /= (int32_t)(ds);
691 fft_buf[offset+j*2] = (int16_t)w;
692 w = (int32_t)fft_buf[offset+j*2+1];
693 w *= (int32_t)(window_coefs[j]);
694 //w /= (int32_t)(ds);
695 fft_buf[offset+j*2+1] = (int16_t)w;
696 }
697 fix_fft(fft_buf+offset, bin_e);
698 if (!peak_hold) {
699 for (j=0; j<bin_len; j++) {
700 ts->avg[j] += real_conj(fft_buf[offset+j*2], fft_buf[offset+j*2+1]);
701 }
702 } else {
703 for (j=0; j<bin_len; j++) {
704 ts->avg[j] = MAX(real_conj(fft_buf[offset+j*2], fft_buf[offset+j*2+1]), ts->avg[j]);
705 }
706 }
707 ts->samples += ds;
708 }
709 }
710 }
711
csv_dbm(struct tuning_state * ts)712 void csv_dbm(struct tuning_state *ts)
713 {
714 int i, len, ds, i1, i2, bw2, bin_count;
715 long tmp;
716 double dbm;
717 len = 1 << ts->bin_e;
718 ds = ts->downsample;
719 /* fix FFT stuff quirks */
720 if (ts->bin_e > 0) {
721 /* nuke DC component (not effective for all windows) */
722 ts->avg[0] = ts->avg[1];
723 /* FFT is translated by 180 degrees */
724 for (i=0; i<len/2; i++) {
725 tmp = ts->avg[i];
726 ts->avg[i] = ts->avg[i+len/2];
727 ts->avg[i+len/2] = tmp;
728 }
729 }
730 /* Hz low, Hz high, Hz step, samples, dbm, dbm, ... */
731 bin_count = (int)((double)len * (1.0 - ts->crop));
732 bw2 = (int)(((double)ts->rate * (double)bin_count) / (len * 2 * ds));
733 fprintf(file, "%i, %i, %.2f, %i, ", ts->freq - bw2, ts->freq + bw2,
734 (double)ts->rate / (double)(len*ds), ts->samples);
735 // something seems off with the dbm math
736 i1 = 0 + (int)((double)len * ts->crop * 0.5);
737 i2 = (len-1) - (int)((double)len * ts->crop * 0.5);
738 for (i=i1; i<=i2; i++) {
739 dbm = (double)ts->avg[i];
740 dbm /= (double)ts->rate;
741 dbm /= (double)ts->samples;
742 dbm = 10 * log10(dbm);
743 fprintf(file, "%.2f, ", dbm);
744 }
745 dbm = (double)ts->avg[i2] / ((double)ts->rate * (double)ts->samples);
746 if (ts->bin_e == 0) {
747 dbm = ((double)ts->avg[0] / \
748 ((double)ts->rate * (double)ts->samples));}
749 dbm = 10 * log10(dbm);
750 fprintf(file, "%.2f\n", dbm);
751 for (i=0; i<len; i++) {
752 ts->avg[i] = 0L;
753 }
754 ts->samples = 0;
755 }
756
main(int argc,char ** argv)757 int main(int argc, char **argv)
758 {
759 #ifndef _WIN32
760 struct sigaction sigact;
761 #endif
762 char *filename = NULL;
763 int i, length, r, opt, wb_mode = 0;
764 int f_set = 0;
765 int gain = AUTO_GAIN; // tenths of a dB
766 int dev_index = 0;
767 int dev_given = 0;
768 int ppm_error = 0;
769 int interval = 10;
770 int fft_threads = 1;
771 int smoothing = 0;
772 int single = 0;
773 int direct_sampling = 0;
774 int offset_tuning = 0;
775 int enable_biastee = 0;
776 double crop = 0.0;
777 char *freq_optarg;
778 time_t next_tick;
779 time_t time_now;
780 time_t exit_time = 0;
781 char t_str[50];
782 struct tm *cal_time;
783 double (*window_fn)(int, int) = rectangle;
784 freq_optarg = "";
785
786 while ((opt = getopt(argc, argv, "f:i:s:t:d:g:p:e:w:c:F:1PDOhT")) != -1) {
787 switch (opt) {
788 case 'f': // lower:upper:bin_size
789 freq_optarg = strdup(optarg);
790 f_set = 1;
791 break;
792 case 'd':
793 dev_index = verbose_device_search(optarg);
794 dev_given = 1;
795 break;
796 case 'g':
797 gain = (int)(atof(optarg) * 10);
798 break;
799 case 'c':
800 crop = atofp(optarg);
801 break;
802 case 'i':
803 interval = (int)round(atoft(optarg));
804 break;
805 case 'e':
806 exit_time = (time_t)((int)round(atoft(optarg)));
807 break;
808 case 's':
809 if (strcmp("avg", optarg) == 0) {
810 smoothing = 0;}
811 if (strcmp("iir", optarg) == 0) {
812 smoothing = 1;}
813 break;
814 case 'w':
815 if (strcmp("rectangle", optarg) == 0) {
816 window_fn = rectangle;}
817 if (strcmp("hamming", optarg) == 0) {
818 window_fn = hamming;}
819 if (strcmp("blackman", optarg) == 0) {
820 window_fn = blackman;}
821 if (strcmp("blackman-harris", optarg) == 0) {
822 window_fn = blackman_harris;}
823 if (strcmp("hann-poisson", optarg) == 0) {
824 window_fn = hann_poisson;}
825 if (strcmp("youssef", optarg) == 0) {
826 window_fn = youssef;}
827 if (strcmp("kaiser", optarg) == 0) {
828 window_fn = kaiser;}
829 if (strcmp("bartlett", optarg) == 0) {
830 window_fn = bartlett;}
831 break;
832 case 't':
833 fft_threads = atoi(optarg);
834 break;
835 case 'p':
836 ppm_error = atoi(optarg);
837 break;
838 case '1':
839 single = 1;
840 break;
841 case 'P':
842 peak_hold = 1;
843 break;
844 case 'D':
845 direct_sampling = 1;
846 break;
847 case 'O':
848 offset_tuning = 1;
849 break;
850 case 'F':
851 boxcar = 0;
852 comp_fir_size = atoi(optarg);
853 break;
854 case 'T':
855 enable_biastee = 1;
856 break;
857 case 'h':
858 default:
859 usage();
860 break;
861 }
862 }
863
864 if (!f_set) {
865 fprintf(stderr, "No frequency range provided.\n");
866 exit(1);
867 }
868
869 if ((crop < 0.0) || (crop > 1.0)) {
870 fprintf(stderr, "Crop value outside of 0 to 1.\n");
871 exit(1);
872 }
873
874 frequency_range(freq_optarg, crop);
875
876 if (tune_count == 0) {
877 usage();}
878
879 if (argc <= optind) {
880 filename = "-";
881 } else {
882 filename = argv[optind];
883 }
884
885 if (interval < 1) {
886 interval = 1;}
887
888 fprintf(stderr, "Reporting every %i seconds\n", interval);
889
890 if (!dev_given) {
891 dev_index = verbose_device_search("0");
892 }
893
894 if (dev_index < 0) {
895 exit(1);
896 }
897
898 r = rtlsdr_open(&dev, (uint32_t)dev_index);
899 if (r < 0) {
900 fprintf(stderr, "Failed to open rtlsdr device #%d.\n", dev_index);
901 exit(1);
902 }
903 #ifndef _WIN32
904 sigact.sa_handler = sighandler;
905 sigemptyset(&sigact.sa_mask);
906 sigact.sa_flags = 0;
907 sigaction(SIGINT, &sigact, NULL);
908 sigaction(SIGTERM, &sigact, NULL);
909 sigaction(SIGQUIT, &sigact, NULL);
910 sigaction(SIGPIPE, &sigact, NULL);
911 #else
912 SetConsoleCtrlHandler( (PHANDLER_ROUTINE) sighandler, TRUE );
913 #endif
914
915 if (direct_sampling) {
916 verbose_direct_sampling(dev, 1);
917 }
918
919 if (offset_tuning) {
920 verbose_offset_tuning(dev);
921 }
922
923 /* Set the tuner gain */
924 if (gain == AUTO_GAIN) {
925 verbose_auto_gain(dev);
926 } else {
927 gain = nearest_gain(dev, gain);
928 verbose_gain_set(dev, gain);
929 }
930
931 verbose_ppm_set(dev, ppm_error);
932
933 rtlsdr_set_bias_tee(dev, enable_biastee);
934 if (enable_biastee)
935 fprintf(stderr, "activated bias-T on GPIO PIN 0\n");
936
937 if (strcmp(filename, "-") == 0) { /* Write log to stdout */
938 file = stdout;
939 #ifdef _WIN32
940 // Is this necessary? Output is ascii.
941 _setmode(_fileno(file), _O_BINARY);
942 #endif
943 } else {
944 file = fopen(filename, "wb");
945 if (!file) {
946 fprintf(stderr, "Failed to open %s\n", filename);
947 exit(1);
948 }
949 }
950
951 /* Reset endpoint before we start reading from it (mandatory) */
952 verbose_reset_buffer(dev);
953
954 /* actually do stuff */
955 rtlsdr_set_sample_rate(dev, (uint32_t)tunes[0].rate);
956 sine_table(tunes[0].bin_e);
957 next_tick = time(NULL) + interval;
958 if (exit_time) {
959 exit_time = time(NULL) + exit_time;}
960 fft_buf = malloc(tunes[0].buf_len * sizeof(int16_t));
961 length = 1 << tunes[0].bin_e;
962 window_coefs = malloc(length * sizeof(int));
963 for (i=0; i<length; i++) {
964 window_coefs[i] = (int)(256*window_fn(i, length));
965 }
966 while (!do_exit) {
967 scanner();
968 time_now = time(NULL);
969 if (time_now < next_tick) {
970 continue;}
971 // time, Hz low, Hz high, Hz step, samples, dbm, dbm, ...
972 cal_time = localtime(&time_now);
973 strftime(t_str, 50, "%Y-%m-%d, %H:%M:%S", cal_time);
974 for (i=0; i<tune_count; i++) {
975 fprintf(file, "%s, ", t_str);
976 csv_dbm(&tunes[i]);
977 }
978 fflush(file);
979 while (time(NULL) >= next_tick) {
980 next_tick += interval;}
981 if (single) {
982 do_exit = 1;}
983 if (exit_time && time(NULL) >= exit_time) {
984 do_exit = 1;}
985 }
986
987 /* clean up */
988
989 if (do_exit) {
990 fprintf(stderr, "\nUser cancel, exiting...\n");}
991 else {
992 fprintf(stderr, "\nLibrary error %d, exiting...\n", r);}
993
994 if (file != stdout) {
995 fclose(file);}
996
997 rtlsdr_close(dev);
998 free(fft_buf);
999 free(window_coefs);
1000 //for (i=0; i<tune_count; i++) {
1001 // free(tunes[i].avg);
1002 // free(tunes[i].buf8);
1003 //}
1004 return r >= 0 ? r : -r;
1005 }
1006
1007 // vim: tabstop=8:softtabstop=8:shiftwidth=8:noexpandtab
1008