From f9d6d22fe2d3d77a1129d9fa2ee293417351354f Mon Sep 17 00:00:00 2001
From: Tatu Peltola <tatu.peltola@aalto.fi>
Date: Fri, 25 Mar 2016 18:00:20 +0200
Subject: [PATCH 1/3] Implement averaged FFT

---
 csdr.c    | 38 ++++++++++++++++++++++++++++++++++++++
 libcsdr.c | 13 +++++++++++++
 libcsdr.h |  2 ++
 3 files changed, 53 insertions(+)

diff --git a/csdr.c b/csdr.c
index 52de1ab..1ed5478 100644
--- a/csdr.c
+++ b/csdr.c
@@ -1299,6 +1299,44 @@ int main(int argc, char *argv[])
 		}
 	}
 
+	if(!strcmp(argv[1],"logaveragepower_cf"))
+	{
+		bigbufs=1;
+		if(argc<=4) return badsyntax("need required parameters (add_db, table_size, avgnumber)"); 
+		float add_db=0;
+		int avgnumber=0;
+		int fft_size=0;
+		
+		sscanf(argv[2],"%g",&add_db);
+		sscanf(argv[3],"%d",&fft_size);
+		sscanf(argv[4],"%d",&avgnumber);
+		
+		if(!getbufsize()) return -2; //dummy
+		if(!sendbufsize(initialize_buffers())) return -2;
+
+		if(fft_size != the_bufsize) return -2;
+		
+		//fprintf(stderr, "logaveragepower_cf %f %d=%d %d\n", add_db, fft_size, the_bufsize, avgnumber);
+		add_db -= 10*log10(avgnumber);
+		for(;;)
+		{
+			int i,n;
+			for(i = 0; i < the_bufsize; i++) {
+				output_buffer[i] = 0;
+			}
+			FEOF_CHECK;
+			for(n = 0; n < avgnumber; n++) {
+				FREAD_C;
+				//fprintf(stderr, "averaged %d\n", n);
+				accumulate_power_cf((complexf*)input_buffer, output_buffer, the_bufsize);	
+			}
+			log_ff(NULL, output_buffer, the_bufsize, add_db);
+			FWRITE_R;
+			TRY_YIELD;
+		}
+		return 0;
+	}
+
 	if(!strcmp(argv[1],"fft_exchange_sides_ff"))
 	{
 		if(argc<=2) return badsyntax("need required parameters (fft_size)");
diff --git a/libcsdr.c b/libcsdr.c
index e3b01bf..0be4cd0 100644
--- a/libcsdr.c
+++ b/libcsdr.c
@@ -949,6 +949,19 @@ void logpower_cf(complexf* input, float* output, int size, float add_db)
 	for(int i=0;i<size;i++) output[i]=10*output[i]+add_db; //@logpower_cf: pass 3
 }
 
+void accumulate_power_cf(complexf* input, float* output, int size)
+{
+	for(int i=0;i<size;i++) output[i] += iof(input,i)*iof(input,i) + qof(input,i)*qof(input,i); //@logpower_cf: pass 1
+	
+}
+
+void log_ff(float* input, float* output, int size, float add_db) {
+	for(int i=0;i<size;i++) output[i]=log10(output[i]); //@logpower_cf: pass 2
+
+	for(int i=0;i<size;i++) output[i]=10*output[i]+add_db; //@logpower_cf: pass 3
+}
+
+
 /*
   _____        _                                            _
  |  __ \      | |                                          (_)
diff --git a/libcsdr.h b/libcsdr.h
index 8bcbc2a..68cfaa8 100644
--- a/libcsdr.h
+++ b/libcsdr.h
@@ -138,6 +138,8 @@ void rational_resampler_get_lowpass_f(float* output, int output_size, int interp
 void apply_window_c(complexf* input, complexf* output, int size, window_t window);
 void apply_window_f(float* input, float* output, int size, window_t window);
 void logpower_cf(complexf* input, float* output, int size, float add_db);
+void accumulate_power_cf(complexf* input, float* output, int size);
+void log_ff(float* input, float* output, int size, float add_db);
 
 typedef struct fractional_decimator_ff_s
 {

From d2a00992275e0fbeeb99865b515f62dac449a10f Mon Sep 17 00:00:00 2001
From: Tatu Peltola <tatu.peltola@aalto.fi>
Date: Sat, 26 Mar 2016 11:44:28 +0200
Subject: [PATCH 2/3] Precalculate FFT window

Before this calculating FFT window was taking much more CPU time than calculating the FFT itself.
---
 csdr.c    |  5 ++++-
 libcsdr.c | 23 +++++++++++++++++++++++
 libcsdr.h |  2 ++
 3 files changed, 29 insertions(+), 1 deletion(-)

diff --git a/csdr.c b/csdr.c
index 1ed5478..c8e8d90 100644
--- a/csdr.c
+++ b/csdr.c
@@ -1245,6 +1245,8 @@ int main(int argc, char *argv[])
 		FFT_PLAN_T* plan=make_fft_c2c(fft_size, windowed, output, 1, benchmark);
 		if(benchmark) fprintf(stderr," done\n");
 		if(octave) printf("setenv(\"GNUTERM\",\"X11 noraise\");y=zeros(1,%d);semilogy(y,\"ydatasource\",\"y\");\n",fft_size);
+		float *windowt;
+		windowt = precalculate_window(fft_size, window);
 		for(;;)
 		{
 			FEOF_CHECK;
@@ -1263,7 +1265,8 @@ int main(int argc, char *argv[])
 				for(int i=0;i<fft_size-every_n_samples;i++) input[i]=input[i+every_n_samples];
 				fread(input+fft_size-every_n_samples, sizeof(complexf), every_n_samples, stdin);
 			}
-			apply_window_c(input,windowed,fft_size,window);
+			//apply_window_c(input,windowed,fft_size,window);
+			apply_precalculated_window_c(input,windowed,fft_size,windowt);
 			fft_execute(plan);
 			if(octave)
 			{
diff --git a/libcsdr.c b/libcsdr.c
index 0be4cd0..5c9051a 100644
--- a/libcsdr.c
+++ b/libcsdr.c
@@ -930,6 +930,29 @@ void apply_window_c(complexf* input, complexf* output, int size, window_t window
 	}
 }
 
+float *precalculate_window(int size, window_t window)
+{
+	float (*window_function)(float)=firdes_get_window_kernel(window);
+	float *windowt;
+	windowt = malloc(sizeof(float) * size);
+	for(int i=0;i<size;i++) //@precalculate_window
+	{
+		float rate=(float)i/(size-1);
+		windowt[i] = window_function(2.0*rate+1.0);
+	}
+	return windowt;
+}
+
+void apply_precalculated_window_c(complexf* input, complexf* output, int size, float *windowt)
+{
+	for(int i=0;i<size;i++) //@apply_precalculated_window_c
+	{
+		iof(output,i)=iof(input,i)*windowt[i];
+		qof(output,i)=qof(input,i)*windowt[i];
+	}
+}
+
+
 void apply_window_f(float* input, float* output, int size, window_t window)
 {
 	float (*window_function)(float)=firdes_get_window_kernel(window);
diff --git a/libcsdr.h b/libcsdr.h
index 68cfaa8..89eb3be 100644
--- a/libcsdr.h
+++ b/libcsdr.h
@@ -135,7 +135,9 @@ typedef struct rational_resampler_ff_s
 rational_resampler_ff_t rational_resampler_ff(float *input, float *output, int input_size, int interpolation, int decimation, float *taps, int taps_length, int last_taps_delay);
 void rational_resampler_get_lowpass_f(float* output, int output_size, int interpolation, int decimation, window_t window);
 
+float *precalculate_window(int size, window_t window);
 void apply_window_c(complexf* input, complexf* output, int size, window_t window);
+void apply_precalculated_window_c(complexf* input, complexf* output, int size, float *windowt);
 void apply_window_f(float* input, float* output, int size, window_t window);
 void logpower_cf(complexf* input, float* output, int size, float add_db);
 void accumulate_power_cf(complexf* input, float* output, int size);

From 38d567d96ef41b57c05cea453c07eb46d89c6a7d Mon Sep 17 00:00:00 2001
From: Tatu Peltola <tatu.peltola@aalto.fi>
Date: Sat, 22 Oct 2016 22:16:47 +0300
Subject: [PATCH 3/3] Fix logaveragepower_cf for FFT sizes below 16384

---
 csdr.c    | 22 +++++++++-------------
 libcsdr.c |  2 +-
 2 files changed, 10 insertions(+), 14 deletions(-)

diff --git a/csdr.c b/csdr.c
index 1598735..18160d6 100644
--- a/csdr.c
+++ b/csdr.c
@@ -1344,27 +1344,23 @@ int main(int argc, char *argv[])
 		sscanf(argv[3],"%d",&fft_size);
 		sscanf(argv[4],"%d",&avgnumber);
 		
-		if(!getbufsize()) return -2; //dummy
-		if(!sendbufsize(initialize_buffers())) return -2;
+		float *input = malloc(sizeof(float)*2 * fft_size);
+		float *output = malloc(sizeof(float) * fft_size);
 
-		if(fft_size != the_bufsize) return -2;
-		
-		//fprintf(stderr, "logaveragepower_cf %f %d=%d %d\n", add_db, fft_size, the_bufsize, avgnumber);
-		add_db -= 10*log10(avgnumber);
+		add_db -= 10.0*log10(avgnumber);
 		for(;;)
 		{
 			int i,n;
-			for(i = 0; i < the_bufsize; i++) {
-				output_buffer[i] = 0;
+			for(i = 0; i < fft_size; i++) {
+				output[i] = 0;
 			}
 			FEOF_CHECK;
 			for(n = 0; n < avgnumber; n++) {
-				FREAD_C;
-				//fprintf(stderr, "averaged %d\n", n);
-				accumulate_power_cf((complexf*)input_buffer, output_buffer, the_bufsize);	
+				fread (input, sizeof(float)*2, fft_size, stdin);
+				accumulate_power_cf((complexf*)input, output, fft_size);
 			}
-			log_ff(NULL, output_buffer, the_bufsize, add_db);
-			FWRITE_R;
+			log_ff(output, output, fft_size, add_db);
+			fwrite (output, sizeof(float), fft_size, stdout);
 			TRY_YIELD;
 		}
 		return 0;
diff --git a/libcsdr.c b/libcsdr.c
index f76bc8e..e38e15e 100644
--- a/libcsdr.c
+++ b/libcsdr.c
@@ -979,7 +979,7 @@ void accumulate_power_cf(complexf* input, float* output, int size)
 }
 
 void log_ff(float* input, float* output, int size, float add_db) {
-	for(int i=0;i<size;i++) output[i]=log10(output[i]); //@logpower_cf: pass 2
+	for(int i=0;i<size;i++) output[i]=log10(input[i]); //@logpower_cf: pass 2
 
 	for(int i=0;i<size;i++) output[i]=10*output[i]+add_db; //@logpower_cf: pass 3
 }