From e4ff3581789cb2c0c4693540474659914f8d0bdc Mon Sep 17 00:00:00 2001 From: Cillian O'Driscoll Date: Thu, 6 Aug 2015 09:31:39 +0100 Subject: [PATCH 1/9] Added seconds_to_skip and header_size Use these configuration parameters to skip ahead in the FileSignalSource Fixed bug in determining file duration --- .../adapters/file_signal_source.cc | 46 +++++++++++++++++-- 1 file changed, 41 insertions(+), 5 deletions(-) diff --git a/src/algorithms/signal_source/adapters/file_signal_source.cc b/src/algorithms/signal_source/adapters/file_signal_source.cc index 780d93ed0..161afd0fe 100644 --- a/src/algorithms/signal_source/adapters/file_signal_source.cc +++ b/src/algorithms/signal_source/adapters/file_signal_source.cc @@ -58,6 +58,8 @@ FileSignalSource::FileSignalSource(ConfigurationInterface* configuration, std::string default_item_type = "short"; std::string default_dump_filename = "./my_capture.dat"; + double default_seconds_to_skip = 0.0; + size_t header_size = 0; samples_ = configuration->property(role + ".samples", 0); sampling_frequency_ = configuration->property(role + ".sampling_frequency", 0); filename_ = configuration->property(role + ".filename", default_filename); @@ -72,6 +74,11 @@ FileSignalSource::FileSignalSource(ConfigurationInterface* configuration, enable_throttle_control_ = configuration->property(role + ".enable_throttle_control", false); std::string s = "InputFilter"; //double IF = configuration->property(s + ".IF", 0.0); + double seconds_to_skip = configuration->property(role + ".seconds_to_skip", default_seconds_to_skip ); + header_size = configuration->property( role + ".header_size", 0 ); + long samples_to_skip = 0; + + bool is_complex = false; if (item_type_.compare("gr_complex") == 0) { @@ -88,6 +95,7 @@ FileSignalSource::FileSignalSource(ConfigurationInterface* configuration, else if (item_type_.compare("ishort") == 0) { item_size_ = sizeof(int16_t); + is_complex = true; } else if (item_type_.compare("byte") == 0) { @@ -96,6 +104,7 @@ FileSignalSource::FileSignalSource(ConfigurationInterface* configuration, else if (item_type_.compare("ibyte") == 0) { item_size_ = sizeof(int8_t); + is_complex = true; } else { @@ -107,6 +116,30 @@ FileSignalSource::FileSignalSource(ConfigurationInterface* configuration, { file_source_ = gr::blocks::file_source::make(item_size_, filename_.c_str(), repeat_); + if( seconds_to_skip > 0 ) + { + samples_to_skip = static_cast< long >( + seconds_to_skip * sampling_frequency_ ); + + if( is_complex ) + { + samples_to_skip *= 2; + } + } + if( header_size > 0 ) + { + samples_to_skip += header_size; + } + + if( samples_to_skip > 0 ) + { + LOG(INFO) << "Skipping " << samples_to_skip << " samples of the input file"; + if( not file_source_->seek( samples_to_skip, SEEK_SET ) ) + { + LOG(INFO) << "Error skipping bytes!"; + } + } + } catch (const std::exception &e) { @@ -174,7 +207,9 @@ FileSignalSource::FileSignalSource(ConfigurationInterface* configuration, if (size > 0) { - samples_ = floor(static_cast(size) / static_cast(item_size()) - ceil(0.002 * static_cast(sampling_frequency_))); //process all the samples available in the file excluding at least the last 1 ms + long bytes_to_skip = samples_to_skip*item_size_; + long bytes_to_process = static_cast(size) - bytes_to_skip; + samples_ = floor(static_cast(bytes_to_process) / static_cast(item_size()) - ceil(0.002 * static_cast(sampling_frequency_))); //process all the samples available in the file excluding at least the last 1 ms } } @@ -182,10 +217,11 @@ FileSignalSource::FileSignalSource(ConfigurationInterface* configuration, double signal_duration_s; signal_duration_s = static_cast(samples_) * ( 1 / static_cast(sampling_frequency_)); - if ((item_type_.compare("gr_complex") != 0) || (item_type_.compare("ishort") != 0) || (item_type_.compare("ibyte") != 0) ) // signal is complex (interleaved) - { - signal_duration_s /= 2; - } + if( is_complex ) + { + signal_duration_s /= 2.0; + } + DLOG(INFO) << "Total number samples to be processed= " << samples_ << " GNSS signal duration= " << signal_duration_s << " [s]"; std::cout << "GNSS signal recorded time to be processed: " << signal_duration_s << " [s]" << std::endl; From 0948e90c613af3f890b7fec88d412acb7752940b Mon Sep 17 00:00:00 2001 From: Cillian O'Driscoll Date: Thu, 6 Aug 2015 09:33:33 +0100 Subject: [PATCH 2/9] Updated gps_l1_ca_dll_pll_read_tracking_dump --- .../matlab/libs/gps_l1_ca_dll_pll_read_tracking_dump.m | 7 ++++++- 1 file changed, 6 insertions(+), 1 deletion(-) diff --git a/src/utils/matlab/libs/gps_l1_ca_dll_pll_read_tracking_dump.m b/src/utils/matlab/libs/gps_l1_ca_dll_pll_read_tracking_dump.m index 913692690..b8e3424e7 100644 --- a/src/utils/matlab/libs/gps_l1_ca_dll_pll_read_tracking_dump.m +++ b/src/utils/matlab/libs/gps_l1_ca_dll_pll_read_tracking_dump.m @@ -35,7 +35,7 @@ function [GNSS_tracking] = gps_l1_ca_dll_pll_read_tracking_dump (filename, count m = nargchk (1,2,nargin); num_float_vars=16; - num_double_vars=1; + num_double_vars=2; double_size_bytes=8; float_size_bytes=4; skip_bytes_each_read=float_size_bytes*num_float_vars+double_size_bytes*num_double_vars; @@ -100,6 +100,9 @@ function [GNSS_tracking] = gps_l1_ca_dll_pll_read_tracking_dump (filename, count bytes_shift=bytes_shift+float_size_bytes; fseek(f,bytes_shift,'bof'); % move to next interleaved float v17 = fread (f, count, 'float64',skip_bytes_each_read-double_size_bytes); + bytes_shift=bytes_shift+double_size_bytes; + fseek(f,bytes_shift,'bof'); % move to next interleaved float + v18 = fread (f, count, 'float64',skip_bytes_each_read-double_size_bytes); fclose (f); %%%%%%%% output vars %%%%%%%% @@ -155,6 +158,7 @@ function [GNSS_tracking] = gps_l1_ca_dll_pll_read_tracking_dump (filename, count carrier_lock_test=v15; var1=v16; var2=v17; + var3=v18; GNSS_tracking.E=E; GNSS_tracking.P=P; @@ -173,5 +177,6 @@ function [GNSS_tracking] = gps_l1_ca_dll_pll_read_tracking_dump (filename, count GNSS_tracking.carrier_lock_test=carrier_lock_test; GNSS_tracking.var1=var1; GNSS_tracking.var2=var2; + GNSS_tracking.var3=var3; end From 5b9683baad55acdbaba775cfc53b105f2ec0b0ad Mon Sep 17 00:00:00 2001 From: Cillian O'Driscoll Date: Thu, 6 Aug 2015 09:34:25 +0100 Subject: [PATCH 3/9] Multiple changes pcps acquisition: including linear corr. Reduced zero_padding_factor in pcps acquisition This was 16, I'm processing data at about 16MHz, this oversampling factor caused memory allocation issues Previously dump file name was not being read from the configuration Added fix for non-zero IF in pcps_acquisition_cc Added linear correlation to pcps acquisition. --- .../galileo_e1_pcps_ambiguous_acquisition.cc | 4 + .../gnuradio_blocks/pcps_acquisition_cc.cc | 83 +++++++++++++++++-- .../pcps_acquisition_fine_doppler_cc.cc | 2 +- 3 files changed, 79 insertions(+), 10 deletions(-) diff --git a/src/algorithms/acquisition/adapters/galileo_e1_pcps_ambiguous_acquisition.cc b/src/algorithms/acquisition/adapters/galileo_e1_pcps_ambiguous_acquisition.cc index 27b602dd8..d063f35ed 100644 --- a/src/algorithms/acquisition/adapters/galileo_e1_pcps_ambiguous_acquisition.cc +++ b/src/algorithms/acquisition/adapters/galileo_e1_pcps_ambiguous_acquisition.cc @@ -97,6 +97,10 @@ GalileoE1PcpsAmbiguousAcquisition::GalileoE1PcpsAmbiguousAcquisition( vector_length_ = sampled_ms_ * samples_per_ms; + if( bit_transition_flag_ ){ + vector_length_ *= 2; + } + code_ = new gr_complex[vector_length_]; if (item_type_.compare("gr_complex") == 0) diff --git a/src/algorithms/acquisition/gnuradio_blocks/pcps_acquisition_cc.cc b/src/algorithms/acquisition/gnuradio_blocks/pcps_acquisition_cc.cc index e6307f74e..fe78e7fcb 100644 --- a/src/algorithms/acquisition/gnuradio_blocks/pcps_acquisition_cc.cc +++ b/src/algorithms/acquisition/gnuradio_blocks/pcps_acquisition_cc.cc @@ -40,6 +40,7 @@ #include #include "gnss_signal_processing.h" #include "control_message_factory.h" +#include using google::LogMessage; @@ -65,8 +66,8 @@ pcps_acquisition_cc::pcps_acquisition_cc( gr::msg_queue::sptr queue, bool dump, std::string dump_filename) : gr::block("pcps_acquisition_cc", - gr::io_signature::make(1, 1, sizeof(gr_complex) * sampled_ms * samples_per_ms), - gr::io_signature::make(0, 0, sizeof(gr_complex) * sampled_ms * samples_per_ms)) + gr::io_signature::make(1, 1, sizeof(gr_complex) * sampled_ms * samples_per_ms * ( bit_transition_flag ? 2 : 1 )), + gr::io_signature::make(0, 0, sizeof(gr_complex) * sampled_ms * samples_per_ms * ( bit_transition_flag ? 2 : 1 )) ) { d_sample_counter = 0; // SAMPLE COUNTER d_active = false; @@ -92,6 +93,24 @@ pcps_acquisition_cc::pcps_acquisition_cc( d_channel = 0; d_doppler_freq = 0.0; + //set_relative_rate( 1.0/d_fft_size ); + + // COD: + // Experimenting with the overlap/save technique for handling bit trannsitions + // The problem: Circular correlation is asynchronous with the received code. + // In effect the first code phase used in the correlation is the current + // estimate of the code phase at the start of the input buffer. If this is 1/2 + // of the code period a bit transition would move all the signal energy into + // adjacent frequency bands at +/- 1/T where T is the integration time. + // + // We can avoid this by doing linear correlation, effectively doubling the + // size of the input buffer and padding the code with zeros. + if( d_bit_transition_flag ) + { + d_fft_size *= 2; + d_max_dwells = 1; + } + d_fft_codes = static_cast(volk_malloc(d_fft_size * sizeof(gr_complex), volk_get_alignment())); d_magnitude = static_cast(volk_malloc(d_fft_size * sizeof(float), volk_get_alignment())); @@ -135,7 +154,17 @@ pcps_acquisition_cc::~pcps_acquisition_cc() void pcps_acquisition_cc::set_local_code(std::complex * code) { - memcpy(d_fft_if->get_inbuf(), code, sizeof(gr_complex) * d_fft_size); + // COD + // Here we want to create a buffer that looks like this: + // [ 0 0 0 ... 0 c_0 c_1 ... c_L] + // where c_i is the local code and there are L zeros and L chips + int offset = 0; + if( d_bit_transition_flag ) + { + std::fill_n( d_fft_if->get_inbuf(), d_samples_per_code, gr_complex( 0.0, 0.0 ) ); + offset = d_samples_per_code; + } + memcpy(d_fft_if->get_inbuf() + offset, code, sizeof(gr_complex) * d_samples_per_code); d_fft_if->execute(); // We need the FFT of local code volk_32fc_conjugate_32fc(d_fft_codes, d_fft_if->get_outbuf(), d_fft_size); } @@ -157,7 +186,7 @@ void pcps_acquisition_cc::init() { d_grid_doppler_wipeoffs[doppler_index] = static_cast(volk_malloc(d_fft_size * sizeof(gr_complex), volk_get_alignment())); int doppler = -static_cast(d_doppler_max) + d_doppler_step * doppler_index; - complex_exp_gen(d_grid_doppler_wipeoffs[doppler_index], d_freq - doppler, d_fs_in, d_fft_size); + complex_exp_gen(d_grid_doppler_wipeoffs[doppler_index], -d_freq - doppler, d_fs_in, d_fft_size); } } @@ -222,6 +251,8 @@ int pcps_acquisition_cc::general_work(int noutput_items, d_sample_counter += d_fft_size * ninput_items[0]; // sample counter consume_each(ninput_items[0]); + //DLOG(INFO) << "Consumed " << ninput_items[0] << " items"; + break; } @@ -232,7 +263,13 @@ int pcps_acquisition_cc::general_work(int noutput_items, unsigned int indext = 0; float magt = 0.0; const gr_complex *in = (const gr_complex *)input_items[0]; //Get the input samples pointer - float fft_normalization_factor = static_cast(d_fft_size) * static_cast(d_fft_size); + + int effective_fft_size = ( d_bit_transition_flag ? d_fft_size/2 : d_fft_size ); + size_t offset = ( d_bit_transition_flag ? effective_fft_size : 0 ); + + float fft_normalization_factor = static_cast(d_fft_size) + * static_cast(d_fft_size); + d_input_power = 0.0; d_mag = 0.0; @@ -273,8 +310,9 @@ int pcps_acquisition_cc::general_work(int noutput_items, d_ifft->execute(); // Search maximum - volk_32fc_magnitude_squared_32f(d_magnitude, d_ifft->get_outbuf(), d_fft_size); - volk_32f_index_max_16u(&indext, d_magnitude, d_fft_size); + size_t offset = ( d_bit_transition_flag ? effective_fft_size : 0 ); + volk_32fc_magnitude_squared_32f(d_magnitude, d_ifft->get_outbuf() + offset, effective_fft_size); + volk_32f_index_max_16u(&indext, d_magnitude, effective_fft_size); // Normalize the maximum value to correct the scale factor introduced by FFTW magt = d_magnitude[indext] / (fft_normalization_factor * fft_normalization_factor); @@ -309,9 +347,19 @@ int pcps_acquisition_cc::general_work(int noutput_items, std::stringstream filename; std::streamsize n = 2 * sizeof(float) * (d_fft_size); // complex file write filename.str(""); - filename << "../data/test_statistics_" << d_gnss_synchro->System + + boost::filesystem::path p = d_dump_filename; + filename << p.parent_path().string() + << boost::filesystem::path::preferred_separator + << p.stem().string() + << "_" << d_gnss_synchro->System <<"_" << d_gnss_synchro->Signal << "_sat_" - << d_gnss_synchro->PRN << "_doppler_" << doppler << ".dat"; + << d_gnss_synchro->PRN << "_doppler_" + << doppler + << p.extension().string(); + + DLOG(INFO) << "Writing ACQ out to " << filename.str(); + d_dump_file.open(filename.str().c_str(), std::ios::out | std::ios::binary); d_dump_file.write((char*)d_ifft->get_outbuf(), n); //write directly |abs(x)|^2 in this Doppler bin? d_dump_file.close(); @@ -346,6 +394,8 @@ int pcps_acquisition_cc::general_work(int noutput_items, consume_each(1); + DLOG(INFO) << "Done. Consumed 1 item."; + break; } @@ -402,3 +452,18 @@ int pcps_acquisition_cc::general_work(int noutput_items, return noutput_items; } + + +//void pcps_acquisition_cc::forecast (int noutput_items, gr_vector_int &ninput_items_required) +//{ + //// COD: + //// For zero-padded case we need one extra code period + //if( d_bit_transition_flag ) + //{ + //ninput_items_required[0] = noutput_items*(d_samples_per_code * d_max_dwells + d_samples_per_code); + //} + //else + //{ + //ninput_items_required[0] = noutput_items*d_fft_size*d_max_dwells; + //} +//} diff --git a/src/algorithms/acquisition/gnuradio_blocks/pcps_acquisition_fine_doppler_cc.cc b/src/algorithms/acquisition/gnuradio_blocks/pcps_acquisition_fine_doppler_cc.cc index b3f18bead..0d4743ef2 100644 --- a/src/algorithms/acquisition/gnuradio_blocks/pcps_acquisition_fine_doppler_cc.cc +++ b/src/algorithms/acquisition/gnuradio_blocks/pcps_acquisition_fine_doppler_cc.cc @@ -316,7 +316,7 @@ int pcps_acquisition_fine_doppler_cc::estimate_Doppler(gr_vector_const_void_star { // Direct FFT - int zero_padding_factor = 16; + int zero_padding_factor = 2; 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); From af3d706fd38103a434cf394ff8526f7e645ad556 Mon Sep 17 00:00:00 2001 From: Cillian O'Driscoll Date: Sun, 9 Aug 2015 22:23:10 +0100 Subject: [PATCH 4/9] Updated info on homebrew install --- CMakeLists.txt | 3 +++ README.md | 39 +++++++++++++++++++++++++++++++++++++-- 2 files changed, 40 insertions(+), 2 deletions(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index e728ba8a5..3dd64a0bb 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -364,6 +364,9 @@ if(NOT GNURADIO_RUNTIME_FOUND) message("You can install it easily via Macports.") message("Open a terminal and type:") message("sudo port install gnuradio ") + message("Alternatively you can use homebrew.") + message("brew tap odrisci/gnuradio") + message("brew install gnuradio" ) message(FATAL_ERROR "GNU Radio 3.7.3 or later is required to build gnss-sdr") endif(OS_IS_MACOSX) endif(NOT GNURADIO_RUNTIME_FOUND) diff --git a/README.md b/README.md index 34c93c509..7311bfa5f 100644 --- a/README.md +++ b/README.md @@ -390,7 +390,11 @@ Agree to Xcode license: $ sudo xcodebuild -license ~~~~~~ -Then, you need a package manager. For example, you can [install Macports](http://www.macports.org/install.php "Macports"). If you are upgrading from a previous installation, please follow the [migration rules](http://trac.macports.org/wiki/Migration). +Software pre-requisites can be installed using either [Macports](#macports) or [Homebrew](#homebrew). + +####Macports + +First, [install Macports](http://www.macports.org/install.php). If you are upgrading from a previous installation, please follow the [migration rules](http://trac.macports.org/wiki/Migration). In a terminal, type: @@ -416,12 +420,43 @@ and you can activate a certain version (2.7 works well) by typing: $ sudo port select --set python python27 ~~~~~~ +#### Homebrew + +Instructions for installing gnuradio using [homebrew](http://www.brew.sh) can be found [here](http://github.com/odrisci/homebrew-gnuradio) - please ensure to install all dependencies as required. + +Install Armadillo and dependencies: + +~~~~~~ +$ brew tap homebrew/science +$ brew install cmake hdf5 arpack superlu +$ brew install armadillo +$ brew install glog gflags +~~~~~~ + +#### Build GNSS-SDR + Finally, you are ready to clone the GNSS-SDR repository and build the software: ~~~~~~ $ git clone https://github.com/gnss-sdr/gnss-sdr $ cd gnss-sdr/build +~~~~~~ + +If using Macports, run: + +~~~~~~ +$ cmake -DCMAKE_CXX_COMPILER=/usr/bin/clang++ -DCMAKE_PREFIX_PATH=/opt/local -DUSE_MACPORTS_PYTHON=/opt/local/bin/python ../ +~~~~~~ + +If using homebrew, run: + +~~~~~~ $ cmake ../ +~~~~~~ + +Finally, run: + +~~~~~~ $ make ~~~~~~ @@ -431,7 +466,7 @@ This will create three executables at gnss-sdr/install, namely ```gnss-sdr```, ` $ sudo make install ~~~~~~ - +Note, it is advisable not to run the install step in a homebrew environment. The documentation can be built by: From 086118c4aaf268e64c581904fadbd98197b25584 Mon Sep 17 00:00:00 2001 From: Cillian O'Driscoll Date: Sat, 24 Oct 2015 09:39:10 +0100 Subject: [PATCH 5/9] Improvements to gnss_flowgraph 1) Fixed a bug in gnss_signal::compare The signal string comparison was not correct. 2) Fix for gnss_flowgraph when satellites are specified 3) Better handling of acquisition in flow graph Now PRNs are alternated amongst channels, so we don't get 'stuck' on a particular PRN, and we don't miss out on PRNs that are assigned to channels that come after a channel that is always in the acquisition state. 4) Added SYS.prns to config file. Now we can specify which PRNs are available in the configuration, rather than re-compiling. Particularly useful for debugging, or experimenting with files where you already know which signals are present --- src/core/receiver/gnss_flowgraph.cc | 73 +++++++++++++++++++++-- src/core/system_parameters/gnss_signal.cc | 2 +- 2 files changed, 70 insertions(+), 5 deletions(-) diff --git a/src/core/receiver/gnss_flowgraph.cc b/src/core/receiver/gnss_flowgraph.cc index 1d1c142db..3cdacde9f 100644 --- a/src/core/receiver/gnss_flowgraph.cc +++ b/src/core/receiver/gnss_flowgraph.cc @@ -43,6 +43,8 @@ #include "gnss_block_interface.h" #include "channel_interface.h" #include "gnss_block_factory.h" +#include +#include #define GNSS_SDR_ARRAY_SIGNAL_CONDITIONER_CHANNELS 8 @@ -313,12 +315,12 @@ void GNSSFlowgraph::connect() } channels_.at(i)->set_signal(available_GNSS_signals_.front()); LOG(INFO) << "Channel " << i << " assigned to " << available_GNSS_signals_.front(); - available_GNSS_signals_.pop_front(); channels_.at(i)->start(); if (channels_state_[i] == 1) { channels_.at(i)->start_acquisition(); + available_GNSS_signals_.pop_front(); LOG(INFO) << "Channel " << i << " connected to observables and ready for acquisition"; } else @@ -413,13 +415,20 @@ void GNSSFlowgraph::apply_action(unsigned int who, unsigned int what) LOG(INFO) << "Channel " << who << " ACQ SUCCESS satellite " << channels_.at(who)->get_signal().get_satellite(); channels_state_[who] = 2; acq_channels_count_--; - if (acq_channels_count_ < max_acq_channels_) + if (!available_GNSS_signals_.empty() && acq_channels_count_ < max_acq_channels_) { for (unsigned int i = 0; i < channels_count_; i++) { if (channels_state_[i] == 0) { channels_state_[i] = 1; + while (channels_.at(i)->get_signal().get_signal_str().compare(available_GNSS_signals_.front().get_signal_str()) != 0 ) + { + available_GNSS_signals_.push_back(available_GNSS_signals_.front()); + available_GNSS_signals_.pop_front(); + } + channels_.at(i)->set_signal(available_GNSS_signals_.front()); + available_GNSS_signals_.pop_front(); acq_channels_count_++; channels_.at(i)->start_acquisition(); break; @@ -442,6 +451,7 @@ void GNSSFlowgraph::apply_action(unsigned int who, unsigned int what) { channels_state_[who] = 0; channels_.at(who)->standby(); + available_GNSS_signals_.push_back( channels_.at(who)->get_signal() ); } // for (unsigned int i = 0; i < channels_count_; i++) @@ -593,8 +603,55 @@ void GNSSFlowgraph::set_signals_list() 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36}; + std::string sv_list = configuration_->property("Galileo.prns", std::string("") ); + if( sv_list.length() > 0 ) + { + // Reset the available prns: + std::set< unsigned int > tmp_set; + boost::tokenizer<> tok( sv_list ); + std::transform( tok.begin(), tok.end(), std::inserter( tmp_set, tmp_set.begin() ), + boost::lexical_cast ); + + if( tmp_set.size() > 0 ) + { + available_galileo_prn = tmp_set; + } + } + + sv_list = configuration_->property("GPS.prns", std::string("") ); + + if( sv_list.length() > 0 ) + { + // Reset the available prns: + std::set< unsigned int > tmp_set; + boost::tokenizer<> tok( sv_list ); + std::transform( tok.begin(), tok.end(), std::inserter( tmp_set, tmp_set.begin() ), + boost::lexical_cast ); + + if( tmp_set.size() > 0 ) + { + available_gps_prn = tmp_set; + } + } + + sv_list = configuration_->property("SBAS.prns", std::string("") ); + + if( sv_list.length() > 0 ) + { + // Reset the available prns: + std::set< unsigned int > tmp_set; + boost::tokenizer<> tok( sv_list ); + std::transform( tok.begin(), tok.end(), std::inserter( tmp_set, tmp_set.begin() ), + boost::lexical_cast ); + + if( tmp_set.size() > 0 ) + { + available_sbas_prn = tmp_set; + } + } + if ((configuration_->property("Channels_1C.count", 0) > 0) or (default_system.find(std::string("GPS")) != std::string::npos) or (default_signal.compare("1C") == 0) or (configuration_->property("Channels_GPS.count", 0) > 0) ) { /* @@ -691,8 +748,16 @@ void GNSSFlowgraph::set_signals_list() } else { - Gnss_Signal signal_value = Gnss_Signal(Gnss_Satellite(gnss_system, gnss_it->get_satellite().get_PRN()), gnss_signal); - available_GNSS_signals_.remove(signal_value); + Gnss_Signal signal_value = Gnss_Signal(Gnss_Satellite(gnss_system, ( sat != 0 ? sat : gnss_it->get_satellite().get_PRN())), gnss_signal); + if( gnss_it == available_GNSS_signals_.begin() ) + { + available_GNSS_signals_.remove(signal_value); + gnss_it = available_GNSS_signals_.begin(); + } + else + { + available_GNSS_signals_.remove(signal_value); + } available_GNSS_signals_.insert(gnss_it, signal_value); } } diff --git a/src/core/system_parameters/gnss_signal.cc b/src/core/system_parameters/gnss_signal.cc index 3230c10dd..feba0e7d1 100644 --- a/src/core/system_parameters/gnss_signal.cc +++ b/src/core/system_parameters/gnss_signal.cc @@ -73,7 +73,7 @@ bool operator==(const Gnss_Signal &sig1, const Gnss_Signal &sig2) if (sig1.get_satellite() == sig2.get_satellite()) { - if (sig1.get_signal_str().compare(sig1.get_signal_str())) + if (sig1.get_signal_str().compare(sig1.get_signal_str())==0) { equal = true; } From e6d5fe35ea87b95ae75fd3c6cadbd0e44719943a Mon Sep 17 00:00:00 2001 From: Cillian O'Driscoll Date: Sun, 22 Nov 2015 08:00:24 +0000 Subject: [PATCH 6/9] Added check for zero denominator in dll discrim. --- src/algorithms/tracking/libs/tracking_discriminators.cc | 9 ++++++++- 1 file changed, 8 insertions(+), 1 deletion(-) diff --git a/src/algorithms/tracking/libs/tracking_discriminators.cc b/src/algorithms/tracking/libs/tracking_discriminators.cc index 119f40a2b..106b3a3ba 100644 --- a/src/algorithms/tracking/libs/tracking_discriminators.cc +++ b/src/algorithms/tracking/libs/tracking_discriminators.cc @@ -101,7 +101,14 @@ float dll_nc_e_minus_l_normalized(gr_complex early_s1, gr_complex late_s1) float P_early, P_late; P_early = std::abs(early_s1); P_late = std::abs(late_s1); - return (P_early - P_late) / ((P_early + P_late)); + if( P_early + P_late == 0.0 ) + { + return 0.0; + } + else + { + return (P_early - P_late) / ((P_early + P_late)); + } } /* From b95219dbc6144dd23ee9cda66616990b6291ee19 Mon Sep 17 00:00:00 2001 From: Cillian O'Driscoll Date: Thu, 26 Nov 2015 13:15:56 +0000 Subject: [PATCH 7/9] Added two Carrier_rotate correlator utilities 1) Carrier_rotate_and_EPL_volk This is a very efficient VOLK based carrier wipeoff and Early/Prompt/Late correlation. The carrier wipe off is done using the VOLK rotator functions, which means there is no need to generate a local replica sinusoid. 2) Carrier_rotate_and_VEPL_Volk This adds Very Early/Early/Prompt/Late/Very Late correlation with the VOLK rotator implementation of carrier wipe-off --- src/algorithms/tracking/libs/correlator.cc | 47 ++++++++++++++++++++++ src/algorithms/tracking/libs/correlator.h | 26 ++++++++++++ 2 files changed, 73 insertions(+) diff --git a/src/algorithms/tracking/libs/correlator.cc b/src/algorithms/tracking/libs/correlator.cc index e44b45051..0f920a961 100644 --- a/src/algorithms/tracking/libs/correlator.cc +++ b/src/algorithms/tracking/libs/correlator.cc @@ -151,3 +151,50 @@ void Correlator::Carrier_wipeoff_and_EPL_volk_custom(int signal_length_samples, volk_cw_epl_corr_u(input, carrier, E_code, P_code, L_code, E_out, P_out, L_out, signal_length_samples); } #endif +void Correlator::Carrier_rotate_and_EPL_volk(int signal_length_samples, + const gr_complex* input, + gr_complex *phase_as_complex, + gr_complex phase_inc_as_complex, + const gr_complex* E_code, + const gr_complex* P_code, + const gr_complex* L_code, + gr_complex* E_out, + gr_complex* P_out, + gr_complex* L_out ) +{ + gr_complex* bb_signal = static_cast(volk_malloc(signal_length_samples * sizeof(gr_complex), volk_get_alignment())); + + volk_32fc_s32fc_x2_rotator_32fc(bb_signal, input, phase_inc_as_complex, phase_as_complex, signal_length_samples); + volk_32fc_x2_dot_prod_32fc(E_out, bb_signal, E_code, signal_length_samples); + volk_32fc_x2_dot_prod_32fc(P_out, bb_signal, P_code, signal_length_samples); + volk_32fc_x2_dot_prod_32fc(L_out, bb_signal, L_code, signal_length_samples); + + volk_free(bb_signal); +} + +void Correlator::Carrier_rotate_and_VEPL_volk(int signal_length_samples, + const gr_complex* input, + gr_complex *phase_as_complex, + gr_complex phase_inc_as_complex, + const gr_complex* VE_code, + const gr_complex* E_code, + const gr_complex* P_code, + const gr_complex* L_code, + const gr_complex* VL_code, + gr_complex* VE_out, + gr_complex* E_out, + gr_complex* P_out, + gr_complex* L_out, + gr_complex* VL_out ) +{ + gr_complex* bb_signal = static_cast(volk_malloc(signal_length_samples * sizeof(gr_complex), volk_get_alignment())); + + volk_32fc_s32fc_x2_rotator_32fc(bb_signal, input, phase_inc_as_complex, phase_as_complex, signal_length_samples); + volk_32fc_x2_dot_prod_32fc(VE_out, bb_signal, VE_code, signal_length_samples); + volk_32fc_x2_dot_prod_32fc(E_out, bb_signal, E_code, signal_length_samples); + volk_32fc_x2_dot_prod_32fc(P_out, bb_signal, P_code, signal_length_samples); + volk_32fc_x2_dot_prod_32fc(L_out, bb_signal, L_code, signal_length_samples); + volk_32fc_x2_dot_prod_32fc(VL_out, bb_signal, VL_code, signal_length_samples); + + volk_free(bb_signal); +} diff --git a/src/algorithms/tracking/libs/correlator.h b/src/algorithms/tracking/libs/correlator.h index 38d2d9b78..a5f2f6623 100644 --- a/src/algorithms/tracking/libs/correlator.h +++ b/src/algorithms/tracking/libs/correlator.h @@ -67,6 +67,32 @@ public: void Carrier_wipeoff_and_EPL_volk_custom(int signal_length_samples, const gr_complex* input, gr_complex* carrier, gr_complex* E_code, gr_complex* P_code, gr_complex* L_code, gr_complex* E_out, gr_complex* P_out, gr_complex* L_out); #endif + void Carrier_rotate_and_EPL_volk(int signal_length_samples, + const gr_complex* input, + gr_complex *phase_as_complex, + gr_complex phase_inc_as_complex, + const gr_complex* E_code, + const gr_complex* P_code, + const gr_complex* L_code, + gr_complex* E_out, + gr_complex* P_out, + gr_complex* L_out ); + + void Carrier_rotate_and_VEPL_volk(int signal_length_samples, + const gr_complex* input, + gr_complex *phase_as_complex, + gr_complex phase_inc_as_complex, + const gr_complex* VE_code, + const gr_complex* E_code, + const gr_complex* P_code, + const gr_complex* L_code, + const gr_complex* VL_code, + gr_complex* VE_out, + gr_complex* E_out, + gr_complex* P_out, + gr_complex* L_out, + gr_complex* VL_out ); + private: unsigned long next_power_2(unsigned long v); }; From 518deb501a1a23e5ccbcdf65583dd45f21e3c6eb Mon Sep 17 00:00:00 2001 From: Cillian O'Driscoll Date: Thu, 26 Nov 2015 13:33:13 +0000 Subject: [PATCH 8/9] Added IF handling to acq and tracking Only GPS L1 CA DLL PLL tracking for now. Still experimental. --- .../gnuradio_blocks/pcps_acquisition_fine_doppler_cc.cc | 2 +- .../tracking/gnuradio_blocks/gps_l1_ca_dll_pll_tracking_cc.cc | 4 ++-- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/src/algorithms/acquisition/gnuradio_blocks/pcps_acquisition_fine_doppler_cc.cc b/src/algorithms/acquisition/gnuradio_blocks/pcps_acquisition_fine_doppler_cc.cc index 0d4743ef2..a0fed41d2 100644 --- a/src/algorithms/acquisition/gnuradio_blocks/pcps_acquisition_fine_doppler_cc.cc +++ b/src/algorithms/acquisition/gnuradio_blocks/pcps_acquisition_fine_doppler_cc.cc @@ -201,7 +201,7 @@ void pcps_acquisition_fine_doppler_cc::update_carrier_wipeoff() doppler_hz = d_config_doppler_min + d_doppler_step*doppler_index; // doppler search steps // compute the carrier doppler wipe-off signal and store it - phase_step_rad = static_cast(GPS_TWO_PI) * doppler_hz / static_cast(d_fs_in); + phase_step_rad = static_cast(GPS_TWO_PI) * ( d_freq + doppler_hz ) / static_cast(d_fs_in); d_grid_doppler_wipeoffs[doppler_index] = new gr_complex[d_fft_size]; fxp_nco(d_grid_doppler_wipeoffs[doppler_index], d_fft_size,0, phase_step_rad); } diff --git a/src/algorithms/tracking/gnuradio_blocks/gps_l1_ca_dll_pll_tracking_cc.cc b/src/algorithms/tracking/gnuradio_blocks/gps_l1_ca_dll_pll_tracking_cc.cc index 3193d335e..5027cbf0a 100644 --- a/src/algorithms/tracking/gnuradio_blocks/gps_l1_ca_dll_pll_tracking_cc.cc +++ b/src/algorithms/tracking/gnuradio_blocks/gps_l1_ca_dll_pll_tracking_cc.cc @@ -297,7 +297,7 @@ void Gps_L1_Ca_Dll_Pll_Tracking_cc::update_local_code() void Gps_L1_Ca_Dll_Pll_Tracking_cc::update_local_carrier() { float sin_f, cos_f; - float phase_step_rad = static_cast(GPS_TWO_PI) * d_carrier_doppler_hz / static_cast(d_fs_in); + float phase_step_rad = static_cast(GPS_TWO_PI) * ( d_if_freq + d_carrier_doppler_hz ) / static_cast(d_fs_in); int phase_step_rad_i = gr::fxpt::float_to_fixed(phase_step_rad); int phase_rad_i = gr::fxpt::float_to_fixed(d_rem_carr_phase_rad); @@ -424,7 +424,7 @@ int Gps_L1_Ca_Dll_Pll_Tracking_cc::general_work (int noutput_items, gr_vector_in //carrier phase accumulator for (K) doppler estimation d_acc_carrier_phase_rad = d_acc_carrier_phase_rad + GPS_TWO_PI * d_carrier_doppler_hz * GPS_L1_CA_CODE_PERIOD; //remanent carrier phase to prevent overflow in the code NCO - d_rem_carr_phase_rad = d_rem_carr_phase_rad + GPS_TWO_PI * d_carrier_doppler_hz * GPS_L1_CA_CODE_PERIOD; + d_rem_carr_phase_rad = d_rem_carr_phase_rad + GPS_TWO_PI * ( d_if_freq + d_carrier_doppler_hz ) * GPS_L1_CA_CODE_PERIOD; d_rem_carr_phase_rad = fmod(d_rem_carr_phase_rad, GPS_TWO_PI); // ################## DLL ########################################################## From 26b18c19ee2f64ef01e1392f7b39138c51297fdc Mon Sep 17 00:00:00 2001 From: Cillian O'Driscoll Date: Sat, 24 Oct 2015 09:42:04 +0100 Subject: [PATCH 9/9] Added a generic tracking_loop_filter class This implements a generic loop filter. Based on the analog PLL filters from Kaplan and Hegarty, with a bilinear (Tustin's) transform from s-plane to z-plane ( 1/s -> T/2 ( 1 + z^-1 )/( 1 - z^-1 ) ) Also added tests. Note the "truth" outputs were derived from an Octave implementation of the loop filter and Octave's builtin filter function --- src/algorithms/tracking/libs/CMakeLists.txt | 1 + .../tracking/libs/tracking_loop_filter.cc | 284 ++++++++++++++++++ .../tracking/libs/tracking_loop_filter.h | 98 ++++++ src/tests/CMakeLists.txt | 3 +- .../arithmetic/tracking_loop_filter_test.cc | 234 +++++++++++++++ 5 files changed, 619 insertions(+), 1 deletion(-) create mode 100644 src/algorithms/tracking/libs/tracking_loop_filter.cc create mode 100644 src/algorithms/tracking/libs/tracking_loop_filter.h create mode 100644 src/tests/arithmetic/tracking_loop_filter_test.cc diff --git a/src/algorithms/tracking/libs/CMakeLists.txt b/src/algorithms/tracking/libs/CMakeLists.txt index a6a51cd01..9374f4cbc 100644 --- a/src/algorithms/tracking/libs/CMakeLists.txt +++ b/src/algorithms/tracking/libs/CMakeLists.txt @@ -39,6 +39,7 @@ set(TRACKING_LIB_SOURCES tracking_2nd_PLL_filter.cc tracking_discriminators.cc tracking_FLL_PLL_filter.cc + tracking_loop_filter.cc ) include_directories( diff --git a/src/algorithms/tracking/libs/tracking_loop_filter.cc b/src/algorithms/tracking/libs/tracking_loop_filter.cc new file mode 100644 index 000000000..3de0521b6 --- /dev/null +++ b/src/algorithms/tracking/libs/tracking_loop_filter.cc @@ -0,0 +1,284 @@ +/*! + * \file tracking_loop_filter.cc + * \brief Generic 1st to 3rd order loop filter implementation + * \author Cillian O'Driscoll, 2015. cillian.odriscoll(at)gmail.com + * + * Class implementing a generic 1st, 2nd or 3rd order loop filter. Based + * on the bilinear transform of the standard Weiner filter. + * + * ------------------------------------------------------------------------- + * + * Copyright (C) 2010-2015 (see AUTHORS file for a list of contributors) + * + * GNSS-SDR is a software defined Global Navigation + * Satellite Systems receiver + * + * This file is part of GNSS-SDR. + * + * GNSS-SDR is free software: you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * GNSS-SDR is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with GNSS-SDR. If not, see . + * + * ------------------------------------------------------------------------- + */ + + +#include "tracking_loop_filter.h" +#include +#include + + +#define MAX_LOOP_ORDER 3 +#define MAX_HISTORY_LENGTH 4 + +Tracking_loop_filter::Tracking_loop_filter( float update_interval, + float noise_bandwidth, + int loop_order, + bool include_last_integrator ) +: d_loop_order( loop_order ), + d_current_index( 0 ), + d_include_last_integrator( include_last_integrator ), + d_noise_bandwidth( noise_bandwidth ), + d_update_interval( update_interval ) +{ + d_inputs.resize( MAX_HISTORY_LENGTH, 0.0 ); + d_outputs.resize( MAX_HISTORY_LENGTH, 0.0 ); + update_coefficients(); +} + +Tracking_loop_filter::Tracking_loop_filter() +: d_loop_order( 2 ), + d_current_index( 0 ), + d_include_last_integrator( false ), + d_noise_bandwidth( 15.0 ), + d_update_interval( 0.001 ) +{ + d_inputs.resize( MAX_HISTORY_LENGTH, 0.0 ); + d_outputs.resize( MAX_HISTORY_LENGTH, 0.0 ); + update_coefficients(); +} + +Tracking_loop_filter::~Tracking_loop_filter() +{ + // Don't need to do anything here +} + +float Tracking_loop_filter::apply( float current_input ) +{ + + // Now apply the filter coefficients: + float result = 0; + + // Hanlde the old outputs first: + for( unsigned int ii=0; ii < d_output_coefficients.size(); ++ii ) + { + result += d_output_coefficients[ii] * d_outputs[ (d_current_index+ii)%MAX_HISTORY_LENGTH ]; + } + + // Now update the index to handle the inputs. + // DO NOT CHANGE THE ORDER OF THE ABOVE AND BELOW CODE + // SNIPPETS!!!!!!! + + // Implementing a sort of circular buffer for the inputs and outputs + // the current input/output is at d_current_index, the nth previous + // input/output is at (d_current_index+n)%d_loop_order + d_current_index--; + if( d_current_index < 0 ) + { + d_current_index += MAX_HISTORY_LENGTH; + } + + d_inputs[d_current_index] = current_input; + + + for( unsigned int ii=0; ii < d_input_coefficients.size(); ++ii ) + { + result += d_input_coefficients[ii] * d_inputs[ (d_current_index+ii)%MAX_HISTORY_LENGTH ]; + } + + + d_outputs[d_current_index] = result; + + + return result; +} + +void Tracking_loop_filter::update_coefficients( void ) +{ + // Analog gains: + float g1; + float g2; + float g3; + + // Natural frequency + float wn; + float T = d_update_interval; + + float zeta = 1/std::sqrt(2); + + // The following is based on the bilinear transform approximation of + // the analog integrator. The loop format is from Kaplan & Hegarty + // Table 5.6. The basic concept is that the loop has a cascade of + // integrators: + // 1 for a 1st order loop + // 2 for a 2nd order loop + // 3 for a 3rd order loop + // The bilinear transform approximates 1/s as + // T/2(1 + z^-1)/(1-z^-1) in the z domain. + + switch( d_loop_order ) + { + case 1: + wn = d_noise_bandwidth*4.0; + g1 = wn; + if( d_include_last_integrator ) + { + d_input_coefficients.resize(2); + d_input_coefficients[0] = g1*T/2.0; + d_input_coefficients[1] = g1*T/2.0; + + d_output_coefficients.resize(1); + d_output_coefficients[0] = 1; + } + else + { + d_input_coefficients.resize(1); + d_input_coefficients[0] = g1; + + d_output_coefficients.resize(0); + } + break; + case 2: + wn = d_noise_bandwidth * (8*zeta)/ (4*zeta*zeta + 1 ); + g1 = wn*wn; + g2 = wn*2*zeta; + if( d_include_last_integrator ) + { + d_input_coefficients.resize(3); + d_input_coefficients[0] = T/2*( g1*T/2 + g2 ); + d_input_coefficients[1] = T*T/2*g1; + d_input_coefficients[2] = T/2*( g1*T/2 - g2 ); + + d_output_coefficients.resize(2); + d_output_coefficients[0] = 2; + d_output_coefficients[1] = -1; + } + else + { + d_input_coefficients.resize(2); + d_input_coefficients[0] = ( g1*T/2.0+g2 ); + d_input_coefficients[1] = g1*T/2-g2; + + d_output_coefficients.resize(1); + d_output_coefficients[0] = 1; + } + break; + + case 3: + wn = d_noise_bandwidth / 0.7845; // From Kaplan + float a3 = 1.1; + float b3 = 2.4; + g1 = wn*wn*wn; + g2 = a3*wn*wn; + g3 = b3*wn; + + if( d_include_last_integrator ) + { + d_input_coefficients.resize(4); + d_input_coefficients[0] = T/2*( g3 + T/2*( g2 + T/2*g1 ) ); + d_input_coefficients[1] = T/2*( -g3 + T/2*( g2 + 3*T/2*g1 ) ); + d_input_coefficients[2] = T/2*( -g3 - T/2*( g2 - 3*T/2*g1 ) ); + d_input_coefficients[3] = T/2*( g3 - T/2*( g2 - T/2*g1 ) ); + + d_output_coefficients.resize(3); + d_output_coefficients[0] = 3; + d_output_coefficients[1] = -3; + d_output_coefficients[2] = 1; + } + else + { + d_input_coefficients.resize(3); + d_input_coefficients[0] = g3 + T/2*( g2 + T/2*g1 ); + d_input_coefficients[1] = g1*T*T/2 -2*g3; + d_input_coefficients[2] = g3 + T/2*( -g2 + T/2*g1 ); + + + d_output_coefficients.resize(2); + d_output_coefficients[0] = 2; + d_output_coefficients[1] = -1; + } + break; + + }; + +} + +void Tracking_loop_filter::set_noise_bandwidth( float noise_bandwidth ) +{ + d_noise_bandwidth = noise_bandwidth; + update_coefficients(); +} + +float Tracking_loop_filter::get_noise_bandwidth( void ) const +{ + return d_noise_bandwidth; +} + +void Tracking_loop_filter::set_update_interval( float update_interval ) +{ + d_update_interval = update_interval; + update_coefficients(); +} + +float Tracking_loop_filter::get_update_interval( void ) const +{ + return d_update_interval; +} + +void Tracking_loop_filter::set_include_last_integrator( bool include_last_integrator ) +{ + d_include_last_integrator = include_last_integrator; + update_coefficients(); +} + +bool Tracking_loop_filter::get_include_last_integrator( void ) const +{ + return d_include_last_integrator; +} + +void Tracking_loop_filter::set_order( int loop_order ) +{ + if( loop_order < 1 || loop_order > MAX_LOOP_ORDER ) + { + LOG(ERROR) << "Ignoring attempt to set loop order to " << loop_order + << ". Maximum allowed order is: " << MAX_LOOP_ORDER + << ". Not changing current value of " << d_loop_order; + + return; + + } + + d_loop_order = loop_order; + update_coefficients(); +} + +int Tracking_loop_filter::get_order( void ) const +{ + return d_loop_order; +} + +void Tracking_loop_filter::initialize( float initial_output ) +{ + d_inputs.assign( MAX_HISTORY_LENGTH, 0.0 ); + d_outputs.assign( MAX_HISTORY_LENGTH, initial_output ); + d_current_index = MAX_HISTORY_LENGTH - 1; +} diff --git a/src/algorithms/tracking/libs/tracking_loop_filter.h b/src/algorithms/tracking/libs/tracking_loop_filter.h new file mode 100644 index 000000000..ac4041f91 --- /dev/null +++ b/src/algorithms/tracking/libs/tracking_loop_filter.h @@ -0,0 +1,98 @@ +/*! + * \file tracking_loop_filter.h + * \brief Generic 1st to 3rd order loop filter implementation + * \author Cillian O'Driscoll, 2015. cillian.odriscoll(at)gmail.com + * + * Class implementing a generic 1st, 2nd or 3rd order loop filter. Based + * on the bilinear transform of the standard Weiner filter. + * + * ------------------------------------------------------------------------- + * + * Copyright (C) 2010-2015 (see AUTHORS file for a list of contributors) + * + * GNSS-SDR is a software defined Global Navigation + * Satellite Systems receiver + * + * This file is part of GNSS-SDR. + * + * GNSS-SDR is free software: you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * GNSS-SDR is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with GNSS-SDR. If not, see . + * + * ------------------------------------------------------------------------- + */ + +#ifndef GNSS_SDR_TRACKING_LOOP_FILTER_H_ +#define GNSS_SDR_TRACKING_LOOP_FILTER_H_ + +#include + + +/*! + * \brief This class implements a generic 1st, 2nd or 3rd order loop filter + * + */ +class Tracking_loop_filter +{ +private: + // Store the last inputs and outputs: + std::vector< float > d_inputs; + std::vector< float > d_outputs; + + // Store the filter coefficients: + std::vector< float > d_input_coefficients; + std::vector< float > d_output_coefficients; + + // The loop order: + int d_loop_order; + + // The current index in the i/o arrays: + int d_current_index; + + // Should the last integrator be included? + bool d_include_last_integrator; + + // The noise bandwidth (in Hz) + // Note this is an approximation only valid when the product of this + // number and the update interval (T) is small. + float d_noise_bandwidth; + + // Loop update interval + float d_update_interval; + + // Compute the filter coefficients: + void update_coefficients(void); + + +public: + float get_noise_bandwidth(void) const; + float get_update_interval(void) const; + bool get_include_last_integrator(void) const; + int get_order(void) const; + + void set_noise_bandwidth( float noise_bandwidth ); + void set_update_interval( float update_interval ); + void set_include_last_integrator( bool include_last_integrator ); + void set_order( int loop_order ); + + void initialize(float initial_output = 0.0); + float apply(float current_input ); + + Tracking_loop_filter(float update_interval, float noise_bandwidth, + int loop_order = 2, + bool include_last_integrator = false ); + + Tracking_loop_filter(); + ~Tracking_loop_filter(); +}; + +#endif diff --git a/src/tests/CMakeLists.txt b/src/tests/CMakeLists.txt index 09093e4ad..2df4ed0ef 100644 --- a/src/tests/CMakeLists.txt +++ b/src/tests/CMakeLists.txt @@ -335,8 +335,9 @@ endif(NOT ${GTEST_DIR_LOCAL}) # add_test(acq_test acq_test) add_executable(trk_test - ${CMAKE_CURRENT_SOURCE_DIR}/single_test_main.cc + ${CMAKE_CURRENT_SOURCE_DIR}/single_test_main.cc ${CMAKE_CURRENT_SOURCE_DIR}/gnss_block/galileo_e1_dll_pll_veml_tracking_test.cc + ${CMAKE_CURRENT_SOURCE_DIR}/arithmetic/tracking_loop_filter_test.cc ) if(NOT ${ENABLE_PACKAGING}) set_property(TARGET trk_test PROPERTY EXCLUDE_FROM_ALL TRUE) diff --git a/src/tests/arithmetic/tracking_loop_filter_test.cc b/src/tests/arithmetic/tracking_loop_filter_test.cc new file mode 100644 index 000000000..10cfe97b0 --- /dev/null +++ b/src/tests/arithmetic/tracking_loop_filter_test.cc @@ -0,0 +1,234 @@ +/*! + * \file tracking_loop_filter_test.cc + * \brief This file implements tests for the general loop filter + * \author Cillian O'Driscoll, 2015. cillian.odriscoll(at)gmail.com + * + * + * ------------------------------------------------------------------------- + * + * Copyright (C) 2010-2015 (see AUTHORS file for a list of contributors) + * + * GNSS-SDR is a software defined Global Navigation + * Satellite Systems receiver + * + * This file is part of GNSS-SDR. + * + * GNSS-SDR is free software: you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * GNSS-SDR is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with GNSS-SDR. If not, see . + * + * ------------------------------------------------------------------------- + */ + +#include "tracking_loop_filter.h" +#include "tracking_2nd_PLL_filter.h" + +#include + +TEST(TrackingLoopFilterTest, FirstOrderLoop) +{ + int loop_order = 1; + float noise_bandwidth = 5.0; + float update_interval = 0.001; + bool include_last_integrator = false; + + Tracking_loop_filter theFilter( update_interval, + noise_bandwidth, + loop_order, + include_last_integrator ); + + EXPECT_EQ( theFilter.get_noise_bandwidth(), noise_bandwidth ); + EXPECT_EQ( theFilter.get_update_interval(), update_interval ); + EXPECT_EQ( theFilter.get_include_last_integrator(), include_last_integrator ); + EXPECT_EQ( theFilter.get_order(), loop_order ); + + std::vector< float > sample_data = { 0, 0, 1.0, 0.0, 0.0, 0.0 }; + + theFilter.initialize( 0.0 ); + + float g1 = noise_bandwidth*4.0; + + float result = 0.0; + for( unsigned int i = 0; i < sample_data.size(); ++i ) + { + result = theFilter.apply( sample_data[i] ); + + ASSERT_FLOAT_EQ( result, sample_data[i]*g1 ); + } +} + +TEST(TrackingLoopFilterTest, FirstOrderLoopWithLastIntegrator) +{ + int loop_order = 1; + float noise_bandwidth = 5.0; + float update_interval = 0.001; + bool include_last_integrator = true; + + Tracking_loop_filter theFilter( update_interval, + noise_bandwidth, + loop_order, + include_last_integrator ); + + EXPECT_EQ( theFilter.get_noise_bandwidth(), noise_bandwidth ); + EXPECT_EQ( theFilter.get_update_interval(), update_interval ); + EXPECT_EQ( theFilter.get_include_last_integrator(), include_last_integrator ); + EXPECT_EQ( theFilter.get_order(), loop_order ); + + std::vector< float > sample_data = { 0, 0, 1.0, 0.0, 0.0, 0.0 }; + std::vector< float > expected_out = { 0.0, 0.0, 0.01, 0.02, 0.02, 0.02 }; + + theFilter.initialize( 0.0 ); + + float g1 = noise_bandwidth*4.0; + + float result = 0.0; + for( unsigned int i = 0; i < sample_data.size(); ++i ) + { + result = theFilter.apply( sample_data[i] ); + ASSERT_NEAR( result, expected_out[i], 1e-4 ); + } + std::cout << std::endl; +} + + + +TEST(TrackingLoopFilterTest, SecondOrderLoop) +{ + int loop_order = 2; + float noise_bandwidth = 5.0; + float update_interval = 0.001; + bool include_last_integrator = false; + + Tracking_loop_filter theFilter( update_interval, + noise_bandwidth, + loop_order, + include_last_integrator ); + + EXPECT_EQ( theFilter.get_noise_bandwidth(), noise_bandwidth ); + EXPECT_EQ( theFilter.get_update_interval(), update_interval ); + EXPECT_EQ( theFilter.get_include_last_integrator(), include_last_integrator ); + EXPECT_EQ( theFilter.get_order(), loop_order ); + + std::vector< float > sample_data = { 0, 0, 1.0, 0.0, 0.0, 0.0 }; + std::vector< float > expected_out = { 0.0, 0.0, 13.37778, 0.0889, 0.0889, 0.0889 }; + + theFilter.initialize( 0.0 ); + + float result = 0.0; + for( unsigned int i = 0; i < sample_data.size(); ++i ) + { + result = theFilter.apply( sample_data[i] ); + + ASSERT_NEAR( result, expected_out[i], 1e-4 ); + } +} + +TEST(TrackingLoopFilterTest, SecondOrderLoopWithLastIntegrator) +{ + int loop_order = 2; + float noise_bandwidth = 5.0; + float update_interval = 0.001; + bool include_last_integrator = true; + + Tracking_loop_filter theFilter( update_interval, + noise_bandwidth, + loop_order, + include_last_integrator ); + + EXPECT_EQ( theFilter.get_noise_bandwidth(), noise_bandwidth ); + EXPECT_EQ( theFilter.get_update_interval(), update_interval ); + EXPECT_EQ( theFilter.get_include_last_integrator(), include_last_integrator ); + EXPECT_EQ( theFilter.get_order(), loop_order ); + + std::vector< float > sample_data = { 0, 0, 1.0, 0.0, 0.0, 0.0 }; + std::vector< float > expected_out = { 0.0, 0.0, 0.006689, 0.013422, 0.013511, 0.013600 }; + + theFilter.initialize( 0.0 ); + + float g1 = noise_bandwidth*4.0; + + float result = 0.0; + for( unsigned int i = 0; i < sample_data.size(); ++i ) + { + result = theFilter.apply( sample_data[i] ); + + ASSERT_NEAR( result, expected_out[i], 1e-4 ); + } + std::cout << std::endl; +} + + +TEST(TrackingLoopFilterTest, ThirdOrderLoop) +{ + int loop_order = 3; + float noise_bandwidth = 5.0; + float update_interval = 0.001; + bool include_last_integrator = false; + + Tracking_loop_filter theFilter( update_interval, + noise_bandwidth, + loop_order, + include_last_integrator ); + + EXPECT_EQ( theFilter.get_noise_bandwidth(), noise_bandwidth ); + EXPECT_EQ( theFilter.get_update_interval(), update_interval ); + EXPECT_EQ( theFilter.get_include_last_integrator(), include_last_integrator ); + EXPECT_EQ( theFilter.get_order(), loop_order ); + + std::vector< float > sample_data = { 0, 0, 1.0, 0.0, 0.0, 0.0 }; + std::vector< float > expected_out = { 0.0, 0.0, 15.31877, 0.04494, 0.04520, 0.04546}; + + theFilter.initialize( 0.0 ); + + float result = 0.0; + for( unsigned int i = 0; i < sample_data.size(); ++i ) + { + result = theFilter.apply( sample_data[i] ); + + ASSERT_NEAR( result, expected_out[i], 1e-4 ); + } +} + +TEST(TrackingLoopFilterTest, ThirdOrderLoopWithLastIntegrator) +{ + int loop_order = 3; + float noise_bandwidth = 5.0; + float update_interval = 0.001; + bool include_last_integrator = true; + + Tracking_loop_filter theFilter( update_interval, + noise_bandwidth, + loop_order, + include_last_integrator ); + + EXPECT_EQ( theFilter.get_noise_bandwidth(), noise_bandwidth ); + EXPECT_EQ( theFilter.get_update_interval(), update_interval ); + EXPECT_EQ( theFilter.get_include_last_integrator(), include_last_integrator ); + EXPECT_EQ( theFilter.get_order(), loop_order ); + + std::vector< float > sample_data = { 0, 0, 1.0, 0.0, 0.0, 0.0 }; + std::vector< float > expected_out = { 0.0, 0.0, 0.007659, 0.015341, 0.015386, 0.015432}; + + theFilter.initialize( 0.0 ); + + float g1 = noise_bandwidth*4.0; + + float result = 0.0; + for( unsigned int i = 0; i < sample_data.size(); ++i ) + { + result = theFilter.apply( sample_data[i] ); + ASSERT_NEAR( result, expected_out[i], 1e-4 ); + } + std::cout << std::endl; +} + +