Merge pull request #21 from tejeez/master

Support for receivers providing real signal

Author: ha7ilm

csdr.c

@@ -836,6 +836,55 @@ int main(int argc, char *argv[])
 		return 0;
 	}
 
+	if(!strcmp(argv[1],"shift_addition_fc"))
+	{
+		bigbufs=1;
+
+		float starting_phase=0;
+		float rate;
+
+		int fd;
+		if(fd=init_fifo(argc,argv))
+		{
+			while(!read_fifo_ctl(fd,"%g\n",&rate)) usleep(10000);
+		}
+		else
+		{
+			if(argc<=2) return badsyntax("need required parameter (rate)");
+			sscanf(argv[2],"%g",&rate);
+		}
+
+		if(!sendbufsize(initialize_buffers())) return -2;
+		for(;;)
+		{
+			shift_addition_data_t data=shift_addition_init(rate);
+			fprintf(stderr,"shift_addition_fc: reinitialized to %g\n",rate);
+			int remain, current_size;
+			float* ibufptr;
+			float* obufptr;
+			for(;;)
+			{
+				FEOF_CHECK;
+				if(!FREAD_R) break;
+				remain=the_bufsize;
+				ibufptr=input_buffer;
+				obufptr=output_buffer;
+				while(remain)
+				{
+					current_size=(remain>1024)?1024:remain;
+					starting_phase=shift_addition_fc(ibufptr, (complexf*)obufptr, current_size, data, starting_phase);
+					ibufptr+=current_size;
+					obufptr+=current_size*2;
+					remain-=current_size;
+				}
+				FWRITE_C;
+				if(read_fifo_ctl(fd,"%g\n",&rate)) break;
+				TRY_YIELD;
+			}
+		}
+		return 0;
+	}
+
 	if(!strcmp(argv[1],"shift_addition_cc_test"))
 	{
 		if(argc<=2) return badsyntax("need required parameter (rate)");
@@ -1495,6 +1544,94 @@ int main(int argc, char *argv[])
 			TRY_YIELD;
 		}
 	}
+
+	if(!strcmp(argv[1],"fft_fc"))
+	{
+		/*
+		For real FFT, the parameter is the number of output complex bins
+		instead of the actual FFT size.
+		Number of input samples used for each FFT is twice the given parameter.
+		This makes it easier to replace fft_cc by fft_fc in some applications. */
+		if(argc<=3) return badsyntax("need required parameters (fft_out_size, out_of_every_n_samples)");
+		int fft_in_size=0, fft_out_size=0;
+		sscanf(argv[2],"%d",&fft_out_size);
+		if(log2n(fft_out_size)==-1) return badsyntax("fft_out_size should be power of 2");
+		fft_in_size = 2*fft_out_size;
+		int every_n_samples;
+		sscanf(argv[3],"%d",&every_n_samples);
+		int benchmark=0;
+		int octave=0;
+		window_t window = WINDOW_DEFAULT;
+		if(argc>=5)
+		{
+			window=firdes_get_window_from_string(argv[4]);
+		}
+		if(argc>=6)
+		{
+			benchmark|=!strcmp("--benchmark",argv[5]);
+			octave|=!strcmp("--octave",argv[5]);
+		}
+		if(argc>=7)
+		{
+			benchmark|=!strcmp("--benchmark",argv[6]);
+			octave|=!strcmp("--octave",argv[6]);
+		}
+
+		if(!initialize_buffers()) return -2;
+		sendbufsize(fft_out_size);
+
+		//make FFT plan
+		float* input=(float*)fft_malloc(sizeof(float)*fft_in_size);
+		float* windowed=(float*)fft_malloc(sizeof(float)*fft_in_size);
+		complexf* output=(complexf*)fft_malloc(sizeof(complexf)*fft_out_size);
+		if(benchmark) fprintf(stderr,"fft_cc: benchmarking...");
+		FFT_PLAN_T* plan=make_fft_r2c(fft_in_size, windowed, output, 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); // TODO
+		float *windowt;
+		windowt = precalculate_window(fft_in_size, window);
+		for(;;)
+		{
+			FEOF_CHECK;
+			if(every_n_samples>fft_in_size)
+			{
+				fread(input, sizeof(float), fft_in_size, stdin);
+				//skipping samples before next FFT (but fseek doesn't work for pipes)
+				for(int seek_remain=every_n_samples-fft_in_size;seek_remain>0;seek_remain-=the_bufsize)
+				{
+					fread(temp_f, sizeof(complexf), MIN_M(the_bufsize,seek_remain), stdin);
+				}
+			}
+			else
+			{
+				//overlapped FFT
+				for(int i=0;i<fft_in_size-every_n_samples;i++) input[i]=input[i+every_n_samples];
+				fread(input+fft_in_size-every_n_samples, sizeof(float), every_n_samples, stdin);
+			}
+			//apply_window_c(input,windowed,fft_size,window);
+			apply_precalculated_window_f(input,windowed,fft_in_size,windowt);
+			fft_execute(plan);
+			if(octave)
+			{
+#if 0
+			// TODO
+				printf("fftdata=[");
+				//we have to swap the two parts of the array to get a valid spectrum
+				for(int i=fft_size/2;i<fft_size;i++) printf("(%g)+(%g)*i ",iof(output,i),qof(output,i));
+				for(int i=0;i<fft_size/2;i++) printf("(%g)+(%g)*i ",iof(output,i),qof(output,i));
+				printf(
+					"];\n"
+					"y=abs(fftdata);\n"
+					"refreshdata;\n"
+				);
+#endif
+			}
+			else fwrite(output, sizeof(complexf), fft_out_size, stdout);
+			TRY_YIELD;
+		}
+	}
+
+
 	#define LOGPOWERCF_BUFSIZE 64
 	if(!strcmp(argv[1],"logpower_cf"))
 	{

fft_fftw.h

@@ -22,6 +22,7 @@ struct fft_plan_s
 #include "libcsdr.h"
 
 FFT_PLAN_T* make_fft_c2c(int size, complexf* input, complexf* output, int forward, int benchmark);
+FFT_PLAN_T* make_fft_r2c(int size, float* input, complexf* output, int benchmark);
 void fft_execute(FFT_PLAN_T* plan);
 void fft_destroy(FFT_PLAN_T* plan);
 

libcsdr.c

@@ -1246,6 +1246,13 @@ void apply_precalculated_window_c(complexf* input, complexf* output, int size, f
 	}
 }
 
+void apply_precalculated_window_f(float* input, float* output, int size, float *windowt)
+{
+	for(int i=0;i<size;i++) //@apply_precalculated_window_f
+	{
+		output[i] = input[i] * windowt[i];
+	}
+}
 
 void apply_window_f(float* input, float* output, int size, window_t window)
 {

libcsdr.h

@@ -140,6 +140,7 @@ void rational_resampler_get_lowpass_f(float* output, int output_size, int interp
 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_precalculated_window_f(float* input, float* 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);

libcsdr_gpl.c

@@ -51,6 +51,33 @@ float shift_addition_cc(complexf *input, complexf* output, int input_size, shift
 	return starting_phase;
 }
 
+float shift_addition_fc(float *input, complexf* output, int input_size, shift_addition_data_t d, float starting_phase)
+{
+	//The original idea was taken from wdsp:
+	//http://svn.tapr.org/repos_sdr_hpsdr/trunk/W5WC/PowerSDR_HPSDR_mRX_PS/Source/wdsp/shift.c
+
+	//However, this method introduces noise (from floating point rounding errors), which increases until the end of the buffer.
+	//fprintf(stderr, "cosd=%g sind=%g\n", d.cosdelta, d.sindelta);
+	float cosphi=cos(starting_phase);
+	float sinphi=sin(starting_phase);
+	float cosphi_last, sinphi_last;
+	for(int i=0;i<input_size;i++) //@shift_addition_cc: work
+	{
+		iof(output,i)=cosphi*input[i];
+		qof(output,i)=sinphi*input[i];
+		//using the trigonometric addition formulas
+		//cos(phi+delta)=cos(phi)cos(delta)-sin(phi)*sin(delta)
+		cosphi_last=cosphi;
+		sinphi_last=sinphi;
+		cosphi=cosphi_last*d.cosdelta-sinphi_last*d.sindelta;
+		sinphi=sinphi_last*d.cosdelta+cosphi_last*d.sindelta;
+	}
+	starting_phase+=d.rate*PI*input_size;
+	while(starting_phase>PI) starting_phase-=2*PI; //@shift_addition_cc: normalize starting_phase
+	while(starting_phase<-PI) starting_phase+=2*PI;
+	return starting_phase;
+}
+
 shift_addition_data_t shift_addition_init(float rate)
 {
 	rate*=2;

libcsdr_gpl.h

@@ -31,6 +31,7 @@ typedef struct shift_addition_data_s
 } shift_addition_data_t;
 shift_addition_data_t shift_addition_init(float rate);
 float shift_addition_cc(complexf *input, complexf* output, int input_size, shift_addition_data_t d, float starting_phase);
+float shift_addition_fc(float *input, complexf* output, int input_size, shift_addition_data_t d, float starting_phase);
 void shift_addition_cc_test(shift_addition_data_t d);
 
 float agc_ff(float* input, float* output, int input_size, float reference, float attack_rate, float decay_rate, float max_gain, short hang_time, short attack_wait_time, float gain_filter_alpha, float last_gain);