usrp: support DCT in addition to complex FFT

- usrp/evscan (usage, main): treat all arguments but the last as options
  and pass them to "fft"
- usrp/fft.c (fft_complex, do_fft): moved actual FFT calculation to
  separate function
- usrp/fft.c (fft_real, do_fft): added real-valued FFT (DCT) as alternative
- usrp/fft.c (main): undocumented option -a to set the FFT algorithm

usrp/evscan

@@ -2,17 +2,22 @@
 
 usage()
 {
-    echo "usage: $0 base-dir" 1>&2
+    echo "usage: $0 [options] base-dir" 1>&2
     exit 1
 }
 
 
+[ "$1" ] || usage
+opts=
+while [ "$2" ]; do
+    opts="$opts $1"
+    shift
+done
 [ -d "$1" ] || usage
-[ "$2" ] && usage
 
 for f in "$1"/24*; do
     echo -n `basename "$f"` ''
     for n in "$f"/data*; do
-	./fft -s 100 0 20 50 <"$n" || echo "fft failed for $n" 1>&1
+	./fft $opts -s 100 0 20 50 <"$n" || echo "fft failed for $n" 1>&1
     done | ./range -v 2
 done

usrp/fft.c

@@ -9,13 +9,61 @@
 
 #define  DEFAULT_THRESHOLD	100
 
+static int alg = 0;
+
+
+static void fft_complex(int n, const float *re, const float *im, double *res)
+{
+	fftw_complex *in, *out;
+	fftw_plan plan;
+	int i;
+	double a;
+
+	in = fftw_malloc(sizeof(fftw_complex)*n);
+	out = fftw_malloc(sizeof(fftw_complex)*n);
+
+	for (i = 0; i != n; i++) {
+		in[i][0] = re[i];
+		in[i][1] = im[i];
+	}
+
+	plan = fftw_plan_dft_1d(n, in, out, FFTW_FORWARD, FFTW_ESTIMATE);
+	fftw_execute(plan);
+
+	for (i = 0; i != n; i++) {
+		a = hypot(out[i][0], out[i][1])/n;
+		a = a*a;
+		res[i] = a;
+	}
+}
+
+
+static void fft_real(int n, const float *re, const float *im, double *res)
+{
+	double *in;
+	fftw_plan plan;
+	int i;
+	double a ;
+
+	in = fftw_malloc(sizeof(double)*n);
+	(void) a;
+
+	for (i = 0; i != n; i++)
+		in[i] = hypot(re[i], im[i]);
+
+	plan =  fftw_plan_r2r_1d(n, in, res, FFTW_REDFT10, FFTW_ESTIMATE);
+	fftw_execute(plan);
+
+	for (i = 0; i != n; i++)
+		res[i] /= n;
+}
+
 
 static void do_fft(int skip, int dump, int low, int high, double threshold)
 {
 	float c[2];
 	float *re = NULL, *im = NULL;
-	fftw_complex *in, *out;
+	double *res;
-	fftw_plan plan;
 	int e = 0, n = 0;
 	int i;
 	double a;
@@ -58,23 +106,26 @@ static void do_fft(int skip, int dump, int low, int high, double threshold)
 	im += skip;
 	n -= skip;
 
-	in = fftw_malloc(sizeof(fftw_complex)*n);
+	res = malloc(n*sizeof(double));
-	out = fftw_malloc(sizeof(fftw_complex)*n);
+	if (!res) {
-
+		perror("malloc");
-	for (i = 0; i != n; i++) {
+		exit(1);
-		in[i][0] = re[i];
-		in[i][1] = im[i];
 	}
 
-	plan = fftw_plan_dft_1d(n, in, out, FFTW_FORWARD, FFTW_ESTIMATE);
+	switch (alg) {
-	fftw_execute(plan);
+	case 0:
+		fft_complex(n, re, im, res);
+		break;
+	case 1:
+		fft_real(n, re, im, res);
+		break;
+	default:
+		abort();
+	}
 
 	if (dump) {
-		for (i = 0; i < n; i += 1) {
+		for (i = 0; i < n; i += 1)
-			a = hypot(out[i][0], out[i][1])/n;
+			printf("%d %g\n", i, 10*log(res[i])/log(10));
-//			a = a*a;
-			printf("%d %g\n", i, 10*log(a)/log(10));
-		}
 	} else {
 		double s = 0;
 		double db;
@@ -91,8 +142,7 @@ static void do_fft(int skip, int dump, int low, int high, double threshold)
 		if (low == high)
 			low--;
 		for (i = low; i != high; i++) {
-			a = hypot(out[i][0], out[i][1])/n;
+			a = res[i];
-//			a = a*a;
 			db = 10*log(a)/log(10);
 			if (db >= threshold)
 				s += a;
@@ -124,8 +174,11 @@ int main(int argc, char **argv)
 	double threshold = DEFAULT_THRESHOLD;
 	int c;
 
-	while ((c = getopt(argc, argv, "ds:")) != EOF)
+	while ((c = getopt(argc, argv, "a:ds:")) != EOF)
 		switch (c) {
+		case 'a':
+			alg = atoi(optarg);
+			break;
 		case 'd':
 			dump = 1;
 			break;