Incremental update for the front-end calibration utility

git-svn-id: https://svn.code.sf.net/p/gnss-sdr/code/trunk@383 64b25241-fba3-4117-9849-534c7e92360d

conf/front-end-cal.conf

@@ -6,24 +6,24 @@
 [GNSS-SDR]
 
 ;######### INITIAL RECEIVER POSITIION ######
-;GNSS-SDR.init_latitude_deg=40.74846557442795
+GNSS-SDR.init_latitude_deg=40.74846557442795
-;GNSS-SDR.init_longitude_deg=-73.98593961814200
+GNSS-SDR.init_longitude_deg=-73.98593961814200
-;GNSS-SDR.init_altitude_m=329.11968943169342
+GNSS-SDR.init_altitude_m=329.11968943169342
 
-GNSS-SDR.init_latitude_deg=41.27481478485936
+;GNSS-SDR.init_latitude_deg=41.27481478485936
-GNSS-SDR.init_longitude_deg=1.98753271588628
+;GNSS-SDR.init_longitude_deg=1.98753271588628
-GNSS-SDR.init_altitude_m=25
+;GNSS-SDR.init_altitude_m=25
 
 ;######### GLOBAL OPTIONS ##################
 ;internal_fs_hz: Internal signal sampling frequency after the signal conditioning stage [Hz].
-GNSS-SDR.internal_fs_hz=2000000
+GNSS-SDR.internal_fs_hz=2000020
 
 ;######### CONTROL_THREAD CONFIG ############
 ControlThread.wait_for_flowgraph=false
 
 ;######### SUPL RRLP GPS assistance configuration #####
 GNSS-SDR.SUPL_gps_enabled=true
-GNSS-SDR.SUPL_read_gps_assistance_xml=false
+GNSS-SDR.SUPL_read_gps_assistance_xml=true
 GNSS-SDR.SUPL_gps_ephemeris_server=supl.nokia.com
 GNSS-SDR.SUPL_gps_ephemeris_port=7275
 GNSS-SDR.SUPL_gps_acquisition_server=supl.google.com
@@ -40,13 +40,13 @@ SignalSource.implementation=File_Signal_Source
 SignalSource.AGC_enabled=false
 
 ;#filename: path to file with the captured GNSS signal samples to be processed
-SignalSource.filename=/media/DATALOGGER_/signals/Agilent GPS Generator/New York/2msps.dat
+SignalSource.filename=/media/DATALOGGER_/signals/RTL-SDR/cap_-90dBm_IF15_RF40_usb_peq.dat
 
 ;#item_type: Type and resolution for each of the signal samples. Use only gr_complex in this version.
 SignalSource.item_type=gr_complex
 
 ;#sampling_frequency: Original Signal sampling frequency in [Hz] 
-SignalSource.sampling_frequency=2000000
+SignalSource.sampling_frequency=2000020
 
 ;#freq: RF front-end center frequency in [Hz] 
 SignalSource.freq=1575420000
@@ -155,8 +155,8 @@ InputFilter.grid_density=16
 ;#The following options are used only in Freq_Xlating_Fir_Filter implementation.
 ;#InputFilter.IF is the intermediate frequency (in Hz) shifted down to zero Hz
 
-InputFilter.sampling_frequency=2000000
+InputFilter.sampling_frequency=2000020
-InputFilter.IF=0
+InputFilter.IF=-16242
 
 InputFilter.decimation_factor=1
 
@@ -181,11 +181,11 @@ Acquisition.sampled_ms=1
 ;#implementation: Acquisition algorithm selection for this channel: [GPS_L1_CA_PCPS_Acquisition] or [Galileo_E1_PCPS_Ambiguous_Acquisition]
 Acquisition.implementation=GPS_L1_CA_PCPS_Acquisition
 ;#threshold: Acquisition threshold
-Acquisition.threshold=40
+Acquisition.threshold=60
 ;#doppler_max: Maximum expected Doppler shift [Hz]
-Acquisition.doppler_max=10000
+Acquisition.doppler_max=100000
 ;#doppler_max: Maximum expected Doppler shift [Hz]
-Acquisition.doppler_min=-10000
+Acquisition.doppler_min=-100000
 ;#doppler_step Doppler step in the grid search [Hz]
 Acquisition.doppler_step=250
 ;#maximum dwells

src/algorithms/acquisition/gnuradio_blocks/pcps_acquisition_fine_doppler_cc.cc

@@ -331,70 +331,112 @@ int pcps_acquisition_fine_doppler_cc::compute_and_accumulate_grid(gr_vector_cons
 	return d_fft_size;
 }
 
+inline int pow2roundup (int x)
+{
+    if (x < 0)
+        return 0;
+    --x;
+    x |= x >> 1;
+    x |= x >> 2;
+    x |= x >> 4;
+    x |= x >> 8;
+    x |= x >> 16;
+    return x+1;
+}
+
 int pcps_acquisition_fine_doppler_cc::estimate_Doppler(gr_vector_const_void_star &input_items, int available_samples)
 {
 
+	// Direct FFT
+	int zero_padding_factor=8;
+	int fft_size_extended=d_fft_size*zero_padding_factor;
+
+	gr::fft::fft_complex *fft_operator=new gr::fft::fft_complex(fft_size_extended,true);
+	//zero padding the entire vector
+	memset(fft_operator->get_inbuf(),0,fft_size_extended*sizeof(gr_complex));
+
 	//1. generate local code aligned with the acquisition code phase estimation
 	gr_complex *code_replica;
 	if (posix_memalign((void**)&code_replica, 16, d_fft_size * sizeof(gr_complex)) == 0){};
 	gps_l1_ca_code_gen_complex_sampled(code_replica, d_gnss_synchro->PRN, d_fs_in, 0);
 
 	int shift_index=(int)d_gnss_synchro->Acq_delay_samples;
-	std::cout<<"shift_index="<<shift_index<<std::endl;
+	//std::cout<<"shift_index="<<shift_index<<std::endl;
-	//std::rotate(code_replica,code_replica+shift_index,code_replica+available_samples);
+	// Rotate to align the local code replica using acquisition time delay estimation
+	std::rotate(code_replica,code_replica+(d_fft_size-shift_index),code_replica+d_fft_size-1);
 
 	//2. Perform code wipe-off
 	const gr_complex *in = (const gr_complex *)input_items[0]; //Get the input samples pointer
 
-	volk_32fc_x2_multiply_32fc_u(d_fft_if->get_inbuf(), in, code_replica, d_fft_size);
+	volk_32fc_x2_multiply_32fc_u(fft_operator->get_inbuf(), in, code_replica, d_fft_size);
 
-	// 3. Perform the FFT
+	// 3. Perform the FFT (zero padded!)
-	d_fft_if->execute();
+	fft_operator->execute();
 
 	// 4. Compute the magnitude and find the maximum
 	float* p_tmp_vector;
-	if (posix_memalign((void**)&p_tmp_vector, 16, d_fft_size * sizeof(float)) == 0){};
+	if (posix_memalign((void**)&p_tmp_vector, 16, fft_size_extended * sizeof(float)) == 0){};
-	volk_32fc_magnitude_squared_32f_a(p_tmp_vector, d_fft_if->get_outbuf(), d_fft_size);
+	volk_32fc_magnitude_squared_32f_a(p_tmp_vector, fft_operator->get_outbuf(), fft_size_extended);
 
 	unsigned int tmp_index_freq;
-	volk_32f_index_max_16u_a(&tmp_index_freq,p_tmp_vector,d_fft_size);
+	volk_32f_index_max_16u_a(&tmp_index_freq,p_tmp_vector,fft_size_extended);
 
-	float fftFreqBins[d_fft_size];
+	//std::cout<<"FFT maximum index present at "<<tmp_index_freq<<std::endl;
-	float curr_normal_freq=-1.0;
+
-	for (unsigned int i=0;i<d_fft_size;i++)
+	//case even
+	int counter=0;
+
+	float fftFreqBins[fft_size_extended];
+
+	for (int k=0;k<(fft_size_extended/2);k++)
 	{
-		fftFreqBins[i]=(d_fs_in/2.0)*curr_normal_freq;
+		fftFreqBins[counter]=((d_fs_in/2)*k)/(fft_size_extended/2);
-		curr_normal_freq+=2/d_fft_size;
+		counter++;
 	}
 
+	for (int k=fft_size_extended/2;k>0;k--)
+	{
+		fftFreqBins[counter]=((-d_fs_in/2)*k)/(fft_size_extended/2);
+		counter++;
+	}
 
-	std::stringstream filename;
+	// 5. Update the Doppler estimation in Hz
-	std::streamsize n = sizeof(gr_complex) * (d_fft_size);
+	if (abs(fftFreqBins[tmp_index_freq]-d_gnss_synchro->Acq_doppler_hz)<1000)
+	{
+		d_gnss_synchro->Acq_doppler_hz=(double)fftFreqBins[tmp_index_freq];
+		//std::cout<<"FFT maximum present at "<<fftFreqBins[tmp_index_freq]<<" [Hz]"<<std::endl;
+	}else{
+		//std::cout<<"Error estimating fine frequency Doppler"<<std::endl;
+		//debug log
+//
+//		std::cout<<"FFT maximum present at "<<fftFreqBins[tmp_index_freq]<<" [Hz]"<<std::endl;
+//		std::stringstream filename;
+//		std::streamsize n = sizeof(gr_complex) * (d_fft_size);
+//
+//		filename.str("");
+//		filename << "../data/code_prn_" << d_gnss_synchro->PRN << ".dat";
+//		d_dump_file.open(filename.str().c_str(), std::ios::out
+//				| std::ios::binary);
+//		d_dump_file.write((char*)code_replica, n); //write directly |abs(x)|^2 in this Doppler bin?
+//		d_dump_file.close();
+//
+//		filename.str("");
+//		filename << "../data/signal_prn_" << d_gnss_synchro->PRN << ".dat";
+//		d_dump_file.open(filename.str().c_str(), std::ios::out
+//				| std::ios::binary);
+//		d_dump_file.write((char*)in, n); //write directly |abs(x)|^2 in this Doppler bin?
+//		d_dump_file.close();
+//
+//
+//		n = sizeof(float) * (fft_size_extended);
+//		filename.str("");
+//		filename << "../data/fft_prn_" << d_gnss_synchro->PRN << ".dat";
+//		d_dump_file.open(filename.str().c_str(), std::ios::out
+//				| std::ios::binary);
+//		d_dump_file.write((char*)p_tmp_vector, n); //write directly |abs(x)|^2 in this Doppler bin?
+//		d_dump_file.close();
+	}
 
-	filename.str("");
-	filename << "../data/code_prn_" << d_gnss_synchro->PRN << ".dat";
-	d_dump_file.open(filename.str().c_str(), std::ios::out
-			| std::ios::binary);
-	d_dump_file.write((char*)code_replica, n); //write directly |abs(x)|^2 in this Doppler bin?
-	d_dump_file.close();
-
-	filename.str("");
-	filename << "../data/signal_prn_" << d_gnss_synchro->PRN << ".dat";
-	d_dump_file.open(filename.str().c_str(), std::ios::out
-			| std::ios::binary);
-	d_dump_file.write((char*)in, n); //write directly |abs(x)|^2 in this Doppler bin?
-	d_dump_file.close();
-
-
-	n = sizeof(float) * (d_fft_size);
-	filename.str("");
-	filename << "../data/fft_prn_" << d_gnss_synchro->PRN << ".dat";
-	d_dump_file.open(filename.str().c_str(), std::ios::out
-			| std::ios::binary);
-	d_dump_file.write((char*)p_tmp_vector, n); //write directly |abs(x)|^2 in this Doppler bin?
-	d_dump_file.close();
-
-	std::cout<<"FFT maximum present at "<<fftFreqBins[tmp_index_freq]<<" [Hz]"<<std::endl;
 
 	return d_fft_size;
 }
@@ -444,7 +486,8 @@ int pcps_acquisition_fine_doppler_cc::general_work(int noutput_items,
 			d_state=3;
 		}
 		d_sample_counter+=consumed_samples;
-		consume_each(consumed_samples);
+		//consume_each(consumed_samples);
+		consume_each(0);
 		break;
 	case 3: // Compute test statistics and decide
 		d_input_power=estimate_input_power(input_items);
@@ -462,8 +505,9 @@ int pcps_acquisition_fine_doppler_cc::general_work(int noutput_items,
 				d_state=7; //negative acquisition
 			}
 		}
-		d_sample_counter += ninput_items[0]; // sample counter
+		//d_sample_counter += ninput_items[0]; // sample counter
-		consume_each(ninput_items[0]);
+		//consume_each(ninput_items[0]);
+		consume_each(0);
 		break;
 	case 4: // RedefineGrid
 		free_grid_memory();
@@ -475,7 +519,7 @@ int pcps_acquisition_fine_doppler_cc::general_work(int noutput_items,
 		break;
 	case 5: // Fine doppler estimation
 		DLOG(INFO) << "Performing fine Doppler estimation";
-		//estimate_Doppler(input_items, ninput_items[0]); //disabled in repo
+		estimate_Doppler(input_items, ninput_items[0]); //disabled in repo
 		d_sample_counter += ninput_items[0]; // sample counter
 		consume_each(ninput_items[0]);
 		d_state=6;

src/algorithms/input_filter/adapters/freq_xlating_fir_filter.cc

@@ -181,6 +181,6 @@ void FreqXlatingFirFilter::init()
     for (std::vector<double>::iterator it = taps_d.begin(); it != taps_d.end(); it++)
         {
             taps_.push_back(float(*it));
-            std::cout<<"TAP="<<float(*it)<<std::endl;
+            //std::cout<<"TAP="<<float(*it)<<std::endl;
         }
 }

src/utils/front-end-cal/main.cc

@@ -366,13 +366,19 @@ int main(int argc, char** argv)
     std::map<int,double> f_fs_estimation_Hz_map;
     std::map<int,double> f_ppm_estimation_Hz_map;
 
+	std::cout <<std::setiosflags(std::ios::fixed)<<std::setprecision(2)<<
+			"Doppler analysis results:"<<std::endl;
+
+    std::cout << "SV ID  Measured [Hz]   Predicted [Hz]" <<std::endl;
+
     for (std::map<int,double>::iterator it = doppler_measurements_map.begin() ; it != doppler_measurements_map.end(); ++it)
       {
-    	std::cout << "Doppler measured for (SV=" << it->first<<")="<<it->second<<" [Hz]"<<std::endl;
+    	//std::cout << "Doppler measured for (SV=" << it->first<<")="<<it->second<<" [Hz]"<<std::endl;
     	try{
     		double doppler_estimated_hz;
     		doppler_estimated_hz=front_end_cal.estimate_doppler_from_eph(it->first,current_TOW,lat_deg,lon_deg,altitude_m);
-    		std::cout << "Doppler estimated for (SV=" << it->first<<")="<<doppler_estimated_hz<<" [Hz]"<<std::endl;
+    		//std::cout << "Doppler estimated for (SV=" << it->first<<")="<<doppler_estimated_hz<<" [Hz]"<<std::endl;
+    		std::cout << "  "<<it->first<<"   "<<it->second<<"   "<<doppler_estimated_hz<<std::endl;
     		// 7. Compute front-end IF and sampling frequency estimation
     		// Compare with the measurements and compute clock drift using FE model
     		double estimated_fs_Hz, estimated_f_if_Hz,f_osc_err_ppm;
@@ -384,7 +390,8 @@ int main(int argc, char** argv)
 
     	}catch(int ex)
     	{
-    		std::cout<<"Eph not found for SV "<<it->first<<std::endl;
+    		//std::cout<<"Eph not found for SV "<<it->first<<std::endl;
+    		std::cout << "  "<<it->first<<"   "<<it->second<<"  (Eph not found)"<<std::endl;
     	}
       }
 
@@ -421,7 +428,7 @@ int main(int argc, char** argv)
     delete gnss_synchro;
 
     google::ShutDownCommandLineFlags();
-    std::cout << "GNSS-SDR program ended." << std::endl;
+    std::cout << "GNSS-SDR Front-end calibration program ended." << std::endl;
 
 //    if (global_gps_acq_assist_map.size()>0)
 //    {