diff --git a/conf/mercurio.conf b/conf/mercurio.conf index 1144a9d11..e4f0724ea 100644 --- a/conf/mercurio.conf +++ b/conf/mercurio.conf @@ -7,11 +7,11 @@ ControlThread.wait_for_flowgraph=false ;######### SIGNAL_SOURCE CONFIG ############ SignalSource.implementation=File_Signal_Source -;SignalSource.filename=/media/My Passport/spirent scenario 2/data/sc2_d8.dat -;SignalSource.filename=/home/luis/Project/signals/GPS_cap1.dat -SignalSource.filename=/media/My Passport/KINGSTON (G)/Project Luis/GPSL1_Fs_8MHz_ID_1_CN0_60.dat +;SignalSource.filename=/media/DATALOGGER/signals/spirent scenario 2/data/sc2_d8.dat +;SignalSource.filename=/media/My Passport/KINGSTON (G)/Project Luis/GPSL1_Fs_8MHz_ID_1_CN0_60.dat +SignalSource.filename=/media/DATALOGGER/signals/Agilent GPS Generator/cap2/agilent_cap2.dat SignalSource.item_type=gr_complex -SignalSource.sampling_frequency=8000000 +SignalSource.sampling_frequency=4000000 SignalSource.samples=0 SignalSource.repeat=false SignalSource.dump=false @@ -20,58 +20,72 @@ SignalSource.enable_throttle_control=true ;######### SIGNAL_CONDITIONER CONFIG ############ SignalConditioner.implementation=Pass_Through SignalConditioner.item_type=gr_complex -SignalConditioner.sample_freq_in=8000000 -SignalConditioner.sample_freq_out=8000000 +SignalConditioner.sample_freq_in=4000000 +SignalConditioner.sample_freq_out=4000000 SignalConditioner.dump=false ;######### CHANNELS CONFIGURATION CONFIG ############ -Channels.count=1 +Channels.count=6 ;######### ACQUISITION CONFIG ############ Acquisition.dump=false Acquisition.dump_filename=./acq_dump.dat Acquisition.item_type=gr_complex -Acquisition.fs_in=8000000 +Acquisition.fs_in=4000000 Acquisition.if=0 Acquisition.sampled_ms=1 ;######### ACQUISITION 0 CONFIG ############ Acquisition0.implementation=GPS_L1_CA_PCPS_Acquisition -Acquisition0.threshold=0.006 +Acquisition0.threshold=0.06 Acquisition0.doppler_max=10000 -Acquisition0.doppler_step=250 -Acquisition0.satellite=1 -Acquisition0.repeat_satellite=true +Acquisition0.doppler_step=500 +Acquisition0.satellite=11 +Acquisition0.repeat_satellite=false ;######### ACQUISITION 1 CONFIG ############ -Acquisition1.implementation=GPS_L1_CA_GPS_SDR_Acquisition -Acquisition1.threshold=30 +Acquisition1.implementation=GPS_L1_CA_PCPS_Acquisition +Acquisition1.threshold=0.06 Acquisition1.doppler_max=10000 -Acquisition1.doppler_step=250 -Acquisition1.satellite=16 -Acquisition1.repeat_satellite=true +Acquisition1.doppler_step=500 +Acquisition1.repeat_satellite=false ;######### ACQUISITION 2 CONFIG ############ -Acquisition2.implementation=GPS_L1_CA_TONG_PCPS_Acquisition -Acquisition2.threshold= +Acquisition2.implementation=GPS_L1_CA_PCPS_Acquisition +Acquisition2.threshold=0.06 Acquisition2.doppler_max=10000 -Acquisition2.doppler_step=250 -Acquisition2.satellite=21 -Acquisition2.repeat_satellite=true - +Acquisition2.doppler_step=500 +Acquisition2.repeat_satellite=false ;######### ACQUISITION 3 CONFIG ############ Acquisition3.implementation=GPS_L1_CA_PCPS_Acquisition -Acquisition3.threshold=105000 +Acquisition3.threshold=0.06 +Acquisition3.doppler_max=10000 +Acquisition3.doppler_step=500 +Acquisition3.repeat_satellite=false + +;######### ACQUISITION 3 CONFIG ############ +Acquisition4.implementation=GPS_L1_CA_PCPS_Acquisition +Acquisition4.threshold=0.06 +Acquisition4.doppler_max=10000 +Acquisition4.doppler_step=500 +Acquisition4.repeat_satellite=false + +;######### ACQUISITION 3 CONFIG ############ +Acquisition5.implementation=GPS_L1_CA_PCPS_Acquisition +Acquisition5.threshold=0.06 +Acquisition5.doppler_max=10000 +Acquisition5.doppler_step=500 +Acquisition5.repeat_satellite=false ;######### TRACKING CONFIG ############ Tracking.implementation=GPS_L1_CA_DLL_PLL_Tracking Tracking.item_type=gr_complex -Tracking.vector_length=8000 -Tracking.fs_in=8000000 +Tracking.vector_length=4000 +Tracking.fs_in=4000000 Tracking.if=0 Tracking.dump=true Tracking.dump_filename=./trk_dump.dat @@ -82,7 +96,7 @@ TelemetryDecoder.item_type=gr_complex ;######### OBSERVABLES CONFIG ############ Observables.implementation=GPS_L1_CA_Observables -Observables.fs_in=8000000; +Observables.fs_in=4000000; ;######### PVT CONFIG ############ PVT.implementation=Pass_Through diff --git a/jamroot.jam b/jamroot.jam index dcad194b7..a603700ba 100644 --- a/jamroot.jam +++ b/jamroot.jam @@ -8,10 +8,12 @@ lib gnuradio-usrp ; lib profiler ; lib gsl ; lib gslcblas ; +lib itpp ; project : requirements OMNITHREAD_POSIX -"-std=c++0x" +"-std=c++0x `pkg-config --cflags itpp`" +"`pkg-config --libs itpp`" # src/utils # src/utils/INIReader src/algorithms/acquisition/adapters diff --git a/src/algorithms/acquisition/gnuradio_blocks/gps_l1_ca_pcps_acquisition_cc.cc b/src/algorithms/acquisition/gnuradio_blocks/gps_l1_ca_pcps_acquisition_cc.cc index 513b6b33f..3fa04d36b 100644 --- a/src/algorithms/acquisition/gnuradio_blocks/gps_l1_ca_pcps_acquisition_cc.cc +++ b/src/algorithms/acquisition/gnuradio_blocks/gps_l1_ca_pcps_acquisition_cc.cc @@ -195,6 +195,7 @@ int gps_l1_ca_pcps_acquisition_cc::general_work(int noutput_items, int doppler; unsigned int indext = 0; float magt = 0.0; + float tmp_magt = 0.0; const gr_complex *in = (const gr_complex *)input_items[0]; //Get the input samples pointer bool positive_acquisition = false; int acquisition_message = -1; //0=STOP_CHANNEL 1=ACQ_SUCCEES 2=ACQ_FAIL @@ -215,8 +216,8 @@ int gps_l1_ca_pcps_acquisition_cc::general_work(int noutput_items, { d_input_power += std::abs(in[i]) * std::abs(in[i]); } - d_input_power = d_input_power / ((float)d_fft_size - * (float)d_fft_size); + //d_input_power = d_input_power / ((float)d_fft_size * (float)d_fft_size); + d_input_power = d_input_power / (float)d_fft_size; // 2º- Doppler frequency search loop for (doppler = (int)(-d_doppler_max); doppler < (int)d_doppler_max; doppler @@ -238,11 +239,23 @@ int gps_l1_ca_pcps_acquisition_cc::general_work(int noutput_items, } d_ifft->execute(); - x86_gr_complex_mag(d_ifft->get_outbuf(), d_fft_size); // d_ifft->get_outbuf()=|abs(·)|^2 and the array is converted from CPX->Float - x86_float_max((float*)d_ifft->get_outbuf(), &indext, &magt, - d_fft_size); // find max of |abs(·)|^2 -> index and magt + // ASM accelerated version + //x86_gr_complex_mag(d_ifft->get_outbuf(), d_fft_size); // d_ifft->get_outbuf()=|abs(·)|^2 and the array is converted from CPX->Float + //x86_float_max((float*)d_ifft->get_outbuf(), &indext, &magt, + // d_fft_size); // find max of |abs(·)|^2 -> index and magt - magt = magt / (float)d_fft_size; + // C++ version + indext=0; + magt=0; + for (i = 0; i < d_fft_size; i++) + { + tmp_magt=std::norm(d_ifft->get_outbuf()[i]); + if (tmp_magt>magt){ + magt=tmp_magt; + indext=i; + } + } + //magt = magt / (float)d_fft_size; // Record results to files if (d_dump) { diff --git a/src/algorithms/channel/adapters/channel.cc b/src/algorithms/channel/adapters/channel.cc index 9fcf3aac9..67f1e562b 100644 --- a/src/algorithms/channel/adapters/channel.cc +++ b/src/algorithms/channel/adapters/channel.cc @@ -128,6 +128,8 @@ void Channel::connect(gr_top_block_sptr top_block) top_block->connect(trk_->get_right_block(), 0, nav_->get_left_block(), 0); // channel 1 top_block->connect(trk_->get_right_block(), 1, nav_->get_left_block(), 1); // channel 2 top_block->connect(trk_->get_right_block(), 2, nav_->get_left_block(), 2); // channel 3 + top_block->connect(trk_->get_right_block(), 3, nav_->get_left_block(), 3); // channel 4 + top_block->connect(trk_->get_right_block(), 4, nav_->get_left_block(), 4); // channel 5 DLOG(INFO) << "tracking -> telemetry_decoder"; connected_ = true; diff --git a/src/algorithms/observables/adapters/gps_l1_ca_observables.cc b/src/algorithms/observables/adapters/gps_l1_ca_observables.cc index f212287eb..b09c35f0b 100644 --- a/src/algorithms/observables/adapters/gps_l1_ca_observables.cc +++ b/src/algorithms/observables/adapters/gps_l1_ca_observables.cc @@ -45,13 +45,6 @@ #include #include -/* The extern keyword declares a variable or function and specifies - * that it has external linkage (its name is visible from files other - * than the one in which it's defined). When modifying a variable, - * extern specifies that the variable has static duration (it is allocated - * when the program begins and deallocated when the program ends). - * The variable is defined in main.cc - */ extern concurrent_queue global_gps_nav_msg_queue; using google::LogMessage; diff --git a/src/algorithms/observables/gnuradio_blocks/gps_l1_ca_observables_cc.cc b/src/algorithms/observables/gnuradio_blocks/gps_l1_ca_observables_cc.cc index 943f693d9..4aeda549e 100644 --- a/src/algorithms/observables/gnuradio_blocks/gps_l1_ca_observables_cc.cc +++ b/src/algorithms/observables/gnuradio_blocks/gps_l1_ca_observables_cc.cc @@ -52,18 +52,34 @@ gps_l1_ca_observables_cc::gps_l1_ca_observables_cc(unsigned int nchannels, gr_ms d_dump = dump; d_nchannels = nchannels; d_rinex_printer.set_headers("GNSS-SDR"); //TODO: read filename from config + d_ls_pvt=new gps_l1_ca_ls_pvt(nchannels); + d_ephemeris_clock_s=0.0; + if (d_dump==true) + { + std::stringstream d_dump_filename_str;//create a stringstream to form the dump filename + d_dump_filename_str<<"./data/obs.dat"; //TODO: get filename and path from config in the adapter + d_dump_filename=d_dump_filename_str.str(); + d_dump_file.open(d_dump_filename.c_str(), std::ios::out | std::ios::binary); + } } gps_l1_ca_observables_cc::~gps_l1_ca_observables_cc() { + delete d_ls_pvt; } -bool pairCompare( pair a, pair b) +bool pairCompare_min( pair a, pair b) { - return a.second.last_preamble_index < b.second.last_preamble_index; + return (a.second.preamble_delay_ms+a.second.prn_delay_ms) < (b.second.preamble_delay_ms+b.second.prn_delay_ms); } +bool pairCompare_max( pair a, pair b) +{ + return (a.second.preamble_delay_ms+a.second.prn_delay_ms) > (b.second.preamble_delay_ms+b.second.prn_delay_ms); +} + + int gps_l1_ca_observables_cc::general_work (int noutput_items, gr_vector_int &ninput_items, gr_vector_const_void_star &input_items, gr_vector_void_star &output_items) { @@ -71,7 +87,7 @@ int gps_l1_ca_observables_cc::general_work (int noutput_items, gr_vector_int &ni // CONSTANTS TODO: place in a file const float code_length_s=0.001; //1 ms - const float C_m_ms= GPS_C_m_s/1000; // The speed of light, [m/ms] + const float C_m_ms= GPS_C_m_s/1000.0; // The speed of light, [m/ms] const float nav_sol_period_ms=1000; //--- Find number of samples per spreading code ---------------------------- @@ -84,46 +100,60 @@ int gps_l1_ca_observables_cc::general_work (int noutput_items, gr_vector_int &ni map pseudoranges; map::iterator pseudoranges_iter; + map pseudoranges_dump; + float min_preamble_delay_ms; + float max_preamble_delay_ms; + bool flag_valid_pseudoranges=false; - unsigned long int min_preamble_index; - - float prn_delay_ms; + float pseudoranges_timestamp_ms; float traveltime_ms; float pseudorange_m; + int pseudoranges_reference_sat_ID=0; + for (unsigned int i=0; i(in[i][0].satellite_PRN,in[i][0])); //record the word structure in a map for pseudoranges + gps_words.insert(pair(in[i][0].channel_ID,in[i][0])); //record the word structure in a map for pseudoranges } } if(gps_words.size()>0) { - // find the minimum index (nearest satellite, will be the reference) - gps_words_iter=min_element(gps_words.begin(),gps_words.end(),pairCompare); - min_preamble_index=gps_words_iter->second.last_preamble_index; - //compute the pseudoranges - for(gps_words_iter = gps_words.begin(); gps_words_iter != gps_words.end(); gps_words_iter++) - { - prn_delay_ms=gps_words_iter->second.prn_delay/(float)samples_per_code; - traveltime_ms=(float)(1000*(gps_words_iter->second.last_preamble_index-min_preamble_index))/(float)samples_per_code+GPS_STARTOFFSET_ms+prn_delay_ms; - //cout<<"traveltime ms"<first<<" ="<(gps_words_iter->first,pseudorange_m)); //record the preamble index in a map - } + // find the minimum index (nearest satellite, will be the reference) + gps_words_iter=min_element(gps_words.begin(),gps_words.end(),pairCompare_min); + min_preamble_delay_ms=gps_words_iter->second.preamble_delay_ms+gps_words_iter->second.prn_delay_ms; //[ms] + pseudoranges_timestamp_ms=min_preamble_delay_ms; //save the shortest pseudorange timestamp to compute the RINEX timestamp + pseudoranges_reference_sat_ID=gps_words_iter->second.satellite_PRN; + gps_words_iter=max_element(gps_words.begin(),gps_words.end(),pairCompare_max); + max_preamble_delay_ms=gps_words_iter->second.preamble_delay_ms; //[ms] + if ((max_preamble_delay_ms-min_preamble_delay_ms)<70) flag_valid_pseudoranges=true; + + + //compute the pseudoranges + for(gps_words_iter = gps_words.begin(); gps_words_iter != gps_words.end(); gps_words_iter++) + { + std::cout<<"prn_delay="<second.prn_delay_ms<second.preamble_delay_ms+gps_words_iter->second.prn_delay_ms-min_preamble_delay_ms+GPS_STARTOFFSET_ms; //[ms] + pseudorange_m=traveltime_ms*C_m_ms; // [m] + pseudoranges.insert(pair(gps_words_iter->second.satellite_PRN,pseudorange_m)); //record the preamble index in a map + if (d_dump==true) + { + pseudoranges_dump.insert(pair(gps_words_iter->second.channel_ID,pseudorange_m)); + } + } // write the pseudoranges to RINEX OBS file // 1- need a valid clock - if (d_last_nav_msg.d_satellite_PRN>0) + if (flag_valid_pseudoranges==true and d_last_nav_msg.d_satellite_PRN>0) { - std::cout<<"d_inter_frame_sec_counter="<get_PVT(pseudoranges,d_ephemeris_clock_s+((double)pseudoranges_timestamp_ms-d_ephemeris_timestamp_ms)/1000.0); } - d_inter_frame_sec_counter+=((float)NAVIGATION_OUTPUT_RATE_MS)/1000.0; //TODO: synchronize the gps time of the ephemeris with the obs } //debug @@ -134,29 +164,36 @@ int gps_l1_ca_observables_cc::general_work (int noutput_items, gr_vector_int &ni { cout<<"Pseudoranges =("<first<<","<second<<")"<try_pop(nav_msg)==true) { - - cout<<"New ephemeris record has arrived!"<d_ephemeris[nav_msg.d_channel_ID]=nav_msg; + // **** update pseudoranges clock **** + if (nav_msg.d_satellite_PRN==pseudoranges_reference_sat_ID) + { + d_ephemeris_clock_s=d_last_nav_msg.d_TOW; + d_ephemeris_timestamp_ms=d_last_nav_msg.d_subframe1_timestamp_ms; + } + // **** write ephemeris to RINES NAV file + //d_rinex_printer.LogRinex2Nav(nav_msg); } + if (d_dump==true) + { + float tmp_float=0.0; + for (int i=0;isecond, sizeof(float)); + }else{ + d_dump_file.write((char*)&tmp_float, sizeof(float)); + } + } + } consume_each(1); //one by one return 0; } diff --git a/src/algorithms/observables/gnuradio_blocks/gps_l1_ca_observables_cc.h b/src/algorithms/observables/gnuradio_blocks/gps_l1_ca_observables_cc.h index 8de66f77b..bcf6b2627 100644 --- a/src/algorithms/observables/gnuradio_blocks/gps_l1_ca_observables_cc.h +++ b/src/algorithms/observables/gnuradio_blocks/gps_l1_ca_observables_cc.h @@ -24,6 +24,7 @@ #include "gps_navigation_message.h" #include "rinex_2_1_printer.h" +#include "gps_l1_ca_ls_pvt.h" #include "GPS_L1_CA.h" @@ -56,10 +57,14 @@ private: rinex_printer d_rinex_printer; // RINEX printer class - double d_inter_frame_sec_counter; // counter for seconds between GPS frames - gps_navigation_message d_last_nav_msg; //last navigation message + double d_ephemeris_clock_s; + double d_ephemeris_timestamp_ms; + + gps_l1_ca_ls_pvt *d_ls_pvt; + + public: ~gps_l1_ca_observables_cc (); diff --git a/src/algorithms/observables/libs/gps_l1_ca_ls_pvt.cc b/src/algorithms/observables/libs/gps_l1_ca_ls_pvt.cc new file mode 100644 index 000000000..b3d90c625 --- /dev/null +++ b/src/algorithms/observables/libs/gps_l1_ca_ls_pvt.cc @@ -0,0 +1,316 @@ + +/** + * Copyright notice + */ + +/** + * Author: Javier Arribas, 2011. jarribas(at)cttc.es + */ + +#include +#include +#include + +using namespace itpp; + +#include "gps_l1_ca_ls_pvt.h" + +gps_l1_ca_ls_pvt::gps_l1_ca_ls_pvt(int nchannels) +{ + // init empty ephemerids for all the available GNSS channels + d_nchannels=nchannels; + d_ephemeris=new gps_navigation_message[nchannels]; +} + +gps_l1_ca_ls_pvt::~gps_l1_ca_ls_pvt() +{ + delete d_ephemeris; +} + +vec gps_l1_ca_ls_pvt::e_r_corr(double traveltime, vec X_sat) { + /* + %E_R_CORR Returns rotated satellite ECEF coordinates due to Earth + %rotation during signal travel time + % + %X_sat_rot = e_r_corr(traveltime, X_sat); + % + % Inputs: + % travelTime - signal travel time + % X_sat - satellite's ECEF coordinates + % + % Outputs: + % X_sat_rot - rotated satellite's coordinates (ECEF) + */ + + double Omegae_dot = 7.292115147e-5; // rad/sec + + //--- Find rotation angle -------------------------------------------------- + double omegatau; + omegatau = Omegae_dot * traveltime; + + //--- Make a rotation matrix ----------------------------------------------- + mat R3=zeros(3,3); + R3(0, 0)= cos(omegatau); + R3(0, 1)= sin(omegatau); + R3(0, 2)= 0.0; + R3(1, 0)=-sin(omegatau); + R3(1, 1)= cos(omegatau); + R3(1, 2)=0.0; + R3(2, 0)= 0.0; + R3(2, 1)= 0.0; + R3(2, 2)= 1; + //--- Do the rotation ------------------------------------------------------ + vec X_sat_rot; + X_sat_rot = R3 * X_sat; + return X_sat_rot; + +} +vec gps_l1_ca_ls_pvt::leastSquarePos(mat satpos, vec obs, mat w) { + //Function calculates the Least Square Solution. + + //pos= leastSquarePos(satpos, obs, w); + + // Inputs: + // satpos - Satellites positions (in ECEF system: [X; Y; Z;] + // obs - Observations - the pseudorange measurements to each + // satellite + // w - weigths vector + + // Outputs: + // pos - receiver position and receiver clock error + // (in ECEF system: [X, Y, Z, dt]) + + ////////////////////////////////////////////////////////////////////////// + // el - Satellites elevation angles (degrees) + // az - Satellites azimuth angles (degrees) + // dop - Dilutions Of Precision ([GDOP PDOP HDOP VDOP TDOP]) + + //=== Initialization ======================================================= + int nmbOfIterations = 10; // TODO: include in config + + //double dtr = GPS_PI / 180.0; + + int nmbOfSatellites; + nmbOfSatellites = satpos.cols(); + vec pos = "0.0 0.0 0.0 0.0"; + mat A; + mat omc; + mat az; + mat el; + A=zeros(nmbOfSatellites,4); + omc=zeros(nmbOfSatellites,1); + az=zeros(1,nmbOfSatellites); + el=zeros(1,nmbOfSatellites); + for (int i = 0; i < nmbOfSatellites; i++) { + for (int j = 0; j < 4; j++) { + A.set(i, j, 0.0); + } + omc(i, 0)=0.0; + az(0, i)=0.0; + } + el = az; + mat X = satpos; + + vec Rot_X; + double rho2; + double traveltime; + double trop; + mat mat_tmp; + vec x; + //=== Iteratively find receiver position =================================== + for (int iter = 0; iter < nmbOfIterations; iter++) { + for (int i = 0; i < nmbOfSatellites; i++) { + if (iter == 0) { + //--- Initialize variables at the first iteration -------------- + Rot_X=X.get_col(i); + trop = 0.0; + } else { + //--- Update equations ----------------------------------------- + rho2 = (X(0, i) - pos(0))*(X(0, i) - pos(0)) + (X(1, i) - pos(1))*(X(1, i) - pos(1))+ (X(2,i) - pos(2))*(X(2,i) - pos(2)); + traveltime = sqrt(rho2) / GPS_C_m_s; + //--- Correct satellite position (do to earth rotation) -------- + Rot_X = e_r_corr(traveltime, X.get_col(i)); + + //--- Find the elevation angel of the satellite ---------------- + //[az(i), el(i), dist] = topocent(pos(1:3, :), Rot_X - pos(1:3, :)); + + } + + //--- Apply the corrections ---------------------------------------- + omc(i) = (obs(i) - norm(Rot_X - pos.left(3)) - pos(4) - trop); + //--- Construct the A matrix --------------------------------------- + A.set(i, 0, (-(Rot_X(0) - pos(0))) / obs(i)); + A.set(i, 1, (-(Rot_X(1) - pos(1))) / obs(i)); + A.set(i, 2, (-(Rot_X(2) - pos(2))) / obs(i)); + A.set(i, 3, 1.0); + } + + + // These lines allow the code to exit gracefully in case of any errors + //if (rank(A) != 4) { + // pos.clear(); + // return pos; + //} + + //--- Find position update --------------------------------------------- + x = ls_solve_od(w*A,w*omc); + + //--- Apply position update -------------------------------------------- + pos = pos + x; + + } + return pos; +} + +void gps_l1_ca_ls_pvt::get_PVT(std::map pseudoranges,double GPS_current_time) +{ + std::map::iterator pseudoranges_iter; + mat satpos; + mat W=eye(d_nchannels); //channels weights matrix + vec obs=zeros(d_nchannels); // pseudoranges observation vector + satpos=zeros(3,d_nchannels); //satellite positions matrix + int valid_obs=0; //valid observations counter + for (int i=0; isecond+d_ephemeris[i].d_satClkCorr*GPS_C_m_s; + valid_obs++; + }else{ + // no valid pseudorange for the current channel + W(i,i)=0; // channel de-activated + obs(i)=1; // to avoid algorithm problems (divide by zero) + } + }else{ + // no valid ephemeris for the current channel + W(i,i)=0; // channel de-activated + obs(i)=1; // to avoid algorithm problems (divide by zero) + } + } + std::cout<<"PVT: valid observations="<=4) + { + vec mypos; + mypos=leastSquarePos(satpos,obs,W); + std::cout << "Position ("< +#include +#include + +using namespace itpp; + +class gps_l1_ca_ls_pvt +{ +private: + vec leastSquarePos(mat satpos, vec obs, mat w); + vec e_r_corr(double traveltime, vec X_sat); + //void cart2geo(); + //void topocent(); +public: + int d_nchannels; + gps_navigation_message* d_ephemeris; + double d_pseudoranges_time_ms; + double d_latitude_d; + double d_longitude_d; + double d_height_m; + double d_x_m; + double d_y_m; + double d_z_m; + + gps_l1_ca_ls_pvt(int nchannels); + ~gps_l1_ca_ls_pvt(); + + void get_PVT(std::map pseudoranges,double GPS_current_time); + void cart2geo(double X, double Y, double Z, int elipsoid_selection); +}; + +#endif diff --git a/src/algorithms/observables/libs/jamfile.jam b/src/algorithms/observables/libs/jamfile.jam index 17c357eeb..bd8fde309 100644 --- a/src/algorithms/observables/libs/jamfile.jam +++ b/src/algorithms/observables/libs/jamfile.jam @@ -1,3 +1,4 @@ project : build-dir ../../../../build ; -obj rinex_2_1_printer : rinex_2_1_printer.cc ; \ No newline at end of file +obj rinex_2_1_printer : rinex_2_1_printer.cc ; +obj gps_l1_ca_ls_pvt : gps_l1_ca_ls_pvt.cc ; \ No newline at end of file diff --git a/src/algorithms/observables/libs/rinex_2_1_printer.cc b/src/algorithms/observables/libs/rinex_2_1_printer.cc index f9eadecd3..5f46a9843 100644 --- a/src/algorithms/observables/libs/rinex_2_1_printer.cc +++ b/src/algorithms/observables/libs/rinex_2_1_printer.cc @@ -192,7 +192,7 @@ void rinex_printer::Rinex2ObsHeader() } -void rinex_printer::LogRinex2Obs(gps_navigation_message nav_msg,double interframe_seconds, std::map pseudoranges) +void rinex_printer::LogRinex2Obs(gps_navigation_message nav_msg,double pseudoranges_clock, std::map pseudoranges) { int ss; char sat_vis[36]; @@ -224,7 +224,7 @@ void rinex_printer::LogRinex2Obs(gps_navigation_message nav_msg,double interfram double decimalday,daydecimalhour,decimalhour,decimalmin,decimalsec; double day,hour,minutes,seconds,enterseconds,a; double gpstime; - gpstime=nav_msg.d_master_clock; + gpstime=pseudoranges_clock; //[s] //calculate date of gps time: //Days & weeks between 00h 1 Jan 1970 and 00h 6 Jan 1980 //520 weeks and 12 days. @@ -233,7 +233,7 @@ void rinex_printer::LogRinex2Obs(gps_navigation_message nav_msg,double interfram tmPtr = gmtime(&temps); strftime( cad1, 20, " %y %m %d", tmPtr ); strftime( cad2, 20, " %Y %m %d", tmPtr ); - decimalday=((gpstime+interframe_seconds)/(24*3600));//Dies dins de la semana + decimalday=(gpstime/(24*3600));//Dies dins de la semana daydecimalhour=modf(decimalday,&day);//day=#dies sencers, daydecimalhour=porcio de dia daydecimalhour=daydecimalhour*24;//porcio de dia en hores decimalhour=modf(daydecimalhour,&hour);//hour=hora del dia; decimalhour=porcio d'hora diff --git a/src/algorithms/telemetry_decoder/adapters/gps_l1_ca_telemetry_decoder.cc b/src/algorithms/telemetry_decoder/adapters/gps_l1_ca_telemetry_decoder.cc index e6b342eb1..7e8d96caa 100644 --- a/src/algorithms/telemetry_decoder/adapters/gps_l1_ca_telemetry_decoder.cc +++ b/src/algorithms/telemetry_decoder/adapters/gps_l1_ca_telemetry_decoder.cc @@ -45,13 +45,6 @@ #include #include -/* The extern keyword declares a variable or function and specifies - * that it has external linkage (its name is visible from files other - * than the one in which it's defined). When modifying a variable, - * extern specifies that the variable has static duration (it is allocated - * when the program begins and deallocated when the program ends). - * The variable is defined in main.cc - */ extern concurrent_queue global_gps_nav_msg_queue; using google::LogMessage; diff --git a/src/algorithms/telemetry_decoder/gnuradio_blocks/gps_l1_ca_telemetry_decoder_cc.cc b/src/algorithms/telemetry_decoder/gnuradio_blocks/gps_l1_ca_telemetry_decoder_cc.cc index 6ed358729..7619cb70f 100644 --- a/src/algorithms/telemetry_decoder/gnuradio_blocks/gps_l1_ca_telemetry_decoder_cc.cc +++ b/src/algorithms/telemetry_decoder/gnuradio_blocks/gps_l1_ca_telemetry_decoder_cc.cc @@ -47,7 +47,7 @@ void gps_l1_ca_telemetry_decoder_cc::forecast (int noutput_items, gps_l1_ca_telemetry_decoder_cc::gps_l1_ca_telemetry_decoder_cc(unsigned int satellite, long if_freq, long fs_in, unsigned int vector_length, gr_msg_queue_sptr queue, bool dump) : - gr_block ("gps_navigation_cc", gr_make_io_signature (3, 3, sizeof(float)), + gr_block ("gps_navigation_cc", gr_make_io_signature (5, 5, sizeof(float)), gr_make_io_signature(1, 1, sizeof(gnss_synchro))) { // initialize internal vars d_queue = queue; @@ -106,10 +106,6 @@ int gps_l1_ca_telemetry_decoder_cc::general_work (int noutput_items, gr_vector_i d_sample_counter++; //count for the processed samples const float **in = (const float **) &input_items[0]; //Get the input samples pointer - //std::cout<<"CH1="<d_satellite+1<d_satellite+1<<" with preamble start at "<7000){ - std::cout<<"lost of frame sync"<d_satellite+1< reactions; gps_subframe_fsm_S0(my_context ctx): my_base( ctx ){ - std::cout<<"Enter S0 "< { @@ -21,7 +21,7 @@ public: sc::transition< Ev_gps_word_valid, gps_subframe_fsm_S2 > > reactions; gps_subframe_fsm_S1(my_context ctx): my_base( ctx ){ - std::cout<<"Enter S1 "< > reactions; gps_subframe_fsm_S2(my_context ctx): my_base( ctx ){ - std::cout<<"Enter S2 "<().gps_word_to_subframe(0); } }; @@ -41,7 +41,7 @@ public: sc::transition< Ev_gps_word_valid, gps_subframe_fsm_S4 > > reactions; gps_subframe_fsm_S3(my_context ctx): my_base( ctx ){ - std::cout<<"Enter S3 "<().gps_word_to_subframe(1); } }; @@ -51,7 +51,7 @@ public: sc::transition< Ev_gps_word_valid, gps_subframe_fsm_S5 > > reactions; gps_subframe_fsm_S4(my_context ctx): my_base( ctx ){ - std::cout<<"Enter S4 "<().gps_word_to_subframe(2); } }; @@ -61,7 +61,7 @@ public: sc::transition< Ev_gps_word_valid, gps_subframe_fsm_S6 > > reactions; gps_subframe_fsm_S5(my_context ctx): my_base( ctx ){ - std::cout<<"Enter S5 "<().gps_word_to_subframe(3); } }; @@ -71,7 +71,7 @@ public: sc::transition< Ev_gps_word_valid, gps_subframe_fsm_S7 > > reactions; gps_subframe_fsm_S6(my_context ctx): my_base( ctx ){ - std::cout<<"Enter S6 "<().gps_word_to_subframe(4); } }; @@ -81,7 +81,7 @@ public: sc::transition< Ev_gps_word_valid, gps_subframe_fsm_S8 > > reactions; gps_subframe_fsm_S7(my_context ctx): my_base( ctx ){ - std::cout<<"Enter S7 "<().gps_word_to_subframe(5); } }; @@ -91,7 +91,7 @@ public: sc::transition< Ev_gps_word_valid, gps_subframe_fsm_S9 > > reactions; gps_subframe_fsm_S8(my_context ctx): my_base( ctx ){ - std::cout<<"Enter S8 "<().gps_word_to_subframe(6); } }; @@ -101,7 +101,7 @@ public: sc::transition< Ev_gps_word_valid, gps_subframe_fsm_S10 > > reactions; gps_subframe_fsm_S9(my_context ctx): my_base( ctx ){ - std::cout<<"Enter S9 "<().gps_word_to_subframe(7); } }; @@ -111,7 +111,7 @@ public: sc::transition< Ev_gps_word_valid, gps_subframe_fsm_S11 > > reactions; gps_subframe_fsm_S10(my_context ctx): my_base( ctx ){ - std::cout<<"Enter S10 "<().gps_word_to_subframe(8); } }; @@ -120,11 +120,11 @@ public: typedef sc::transition< Ev_gps_word_preamble, gps_subframe_fsm_S1 > reactions; gps_subframe_fsm_S11(my_context ctx): my_base( ctx ){ - std::cout<<"Completed GPS Subframe!"<().gps_word_to_subframe(9); context< GpsL1CaSubframeFsm >().gps_subframe_to_nav_msg(); //decode the subframe // DECODE SUBFRAME - std::cout<<"Enter S11"<d_preamble_time_ms-6002; + std::cout<<"FSM: set subframe1 preamble timestamp for sat "<push(d_nav); } diff --git a/src/algorithms/telemetry_decoder/libs/gps_l1_ca_subframe_fsm.h b/src/algorithms/telemetry_decoder/libs/gps_l1_ca_subframe_fsm.h index 32847f827..a8f859dac 100644 --- a/src/algorithms/telemetry_decoder/libs/gps_l1_ca_subframe_fsm.h +++ b/src/algorithms/telemetry_decoder/libs/gps_l1_ca_subframe_fsm.h @@ -61,6 +61,8 @@ public: char d_subframe[GPS_SUBFRAME_LENGTH]; char d_GPS_frame_4bytes[GPS_WORD_LENGTH]; + double d_preamble_time_ms; + void gps_word_to_subframe(int position); void gps_subframe_to_nav_msg(); 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 01433c874..14754c599 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 @@ -19,12 +19,11 @@ #include "gps_l1_ca_dll_pll_tracking_cc.h" #include "gps_sdr_signal_processing.h" -//#include "gps_sdr_signaldef.h" #include "gps_sdr_simd.h" #include "gps_sdr_x86.h" #include "control_message_factory.h" - +#include #include #include #include @@ -37,6 +36,11 @@ #include #include +/*! + * \todo Include in definition header file + */ +#define CN0_ESTIMATION_SAMPLES 10 + using google::LogMessage; gps_l1_ca_dll_pll_tracking_cc_sptr @@ -47,11 +51,18 @@ gps_l1_ca_dll_pll_make_tracking_cc(unsigned int satellite, long if_freq, long fs fs_in, vector_length, queue, dump)); } +void gps_l1_ca_dll_pll_tracking_cc::forecast (int noutput_items, + gr_vector_int &ninput_items_required){ + ninput_items_required[0] =d_vector_length*2; //set the required available samples in each call +} gps_l1_ca_dll_pll_tracking_cc::gps_l1_ca_dll_pll_tracking_cc(unsigned int satellite, long if_freq, long fs_in, unsigned int vector_length, gr_msg_queue_sptr queue, bool dump) : - gr_sync_decimator ("gps_l1_ca_dll_pll_tracking_cc", gr_make_io_signature (1, 1, sizeof(gr_complex)), - gr_make_io_signature(3, 3, sizeof(float)),vector_length) { + gr_block ("gps_l1_ca_dll_pll_tracking_cc", gr_make_io_signature (1, 1, sizeof(gr_complex)), + gr_make_io_signature(5, 5, sizeof(float))) { + + //gr_sync_decimator ("gps_l1_ca_dll_pll_tracking_cc", gr_make_io_signature (1, 1, sizeof(gr_complex)), + // gr_make_io_signature(3, 3, sizeof(float)),vector_length) { // initialize internal vars d_queue = queue; d_dump = dump; @@ -61,38 +72,39 @@ gps_l1_ca_dll_pll_tracking_cc::gps_l1_ca_dll_pll_tracking_cc(unsigned int satell d_vector_length = vector_length; // Initialize tracking variables ========================================== - //TODO: Include this setting in configuration file + /*! + * \todo Include PLL and DLL filter setting in configuration file + */ //--- DLL variables -------------------------------------------------------- d_early_late_spc = 0.5; // Define early-late offset (in chips) d_pdi_code = 0.001;// Summation interval for code - d_dllnoisebandwidth=2.0; //Hz + d_dllnoisebandwidth=1; //Hz d_dlldampingratio=0.7; calculate_lopp_coef(&d_tau1_code, &d_tau2_code, d_dllnoisebandwidth, d_dlldampingratio,1.0);// Calculate filter coefficient values //--- PLL variables -------------------------------------------------------- d_pdi_carr = 0.001;// Summation interval for carrier d_plldampingratio=0.7; - d_pllnoisebandwidth=25; + d_pllnoisebandwidth=50; //Calculate filter coefficient values calculate_lopp_coef(&d_tau1_carr, &d_tau2_carr, d_pllnoisebandwidth, d_plldampingratio,0.25);// Calculate filter coefficient values // Initialization of local code replica - // Get a vector with the C/A code sampled 1x/chip d_code_length=1023; - d_ca_code=(gr_complex*)malloc(sizeof(gr_complex)*(d_code_length+2)); + // Get space for a vector with the C/A code replica sampled 1x/chip + d_ca_code=new gr_complex[d_code_length+2]; // Get space for the resampled early / prompt / late local replicas - d_early_code=(gr_complex*)malloc(sizeof(gr_complex)*d_vector_length); - // Get space for the resampled early / prompt / late local replicas - d_prompt_code=(gr_complex*)malloc(sizeof(gr_complex)*d_vector_length); - // Get space for the resampled early / prompt / late local replicas - d_late_code=(gr_complex*)malloc(sizeof(gr_complex)*d_vector_length); + d_early_code= new gr_complex[d_vector_length*2]; + d_prompt_code=new gr_complex[d_vector_length*2]; + d_late_code=new gr_complex[d_vector_length*2]; - d_carr_sign=(gr_complex*)malloc(sizeof(gr_complex)*d_vector_length); - d_bb_sign=(gr_complex*)malloc(sizeof(gr_complex)*d_vector_length); + // space for carrier wipeoff and signal baseband vectors + d_carr_sign=new gr_complex[d_vector_length*2]; + d_bb_sign=new gr_complex[d_vector_length*2]; //--- Perform initializations ------------------------------ @@ -114,167 +126,194 @@ gps_l1_ca_dll_pll_tracking_cc::gps_l1_ca_dll_pll_tracking_cc(unsigned int satell d_old_carr_nco = 0.0; d_old_carr_error = 0.0; - d_absolute_code_phase_chips = 0; + d_absolute_code_phase_samples = 0; // sample synchronization d_sample_counter=0; d_acq_sample_stamp=0; - d_dump_filename="./data/trk_epr.dat"; d_enable_tracking=false; - d_last_seg=0; + d_pull_in=false; + d_last_seg=0; + + d_blksize=d_vector_length; + d_loops_count=0; + + // CN0 estimation and lock detector buffers + d_cn0_estimation_counter=0; + d_P_I_buffer=new float[CN0_ESTIMATION_SAMPLES]; + d_P_Q_buffer=new float[CN0_ESTIMATION_SAMPLES]; + d_carrier_lock_test=1; + d_SNR_SNV=0; + d_SNR_SNV_dB_Hz=0; + d_SNR_MM=0; + d_carrier_lock_fail_counter=0; + d_carrier_lock_threshold=5; + } void gps_l1_ca_dll_pll_tracking_cc::calculate_lopp_coef(float* tau1,float* tau2, float lbw, float zeta, float k){ // Solve natural frequency float Wn; Wn = lbw*8*zeta / (4*zeta*zeta + 1); - // solve for t1 & t2 *tau1 = k / (Wn * Wn); *tau2 = (2.0 * zeta) / Wn; - //std::cout<<"Tau1= "<<*tau1<d_satellite<<" received "; + + if (d_dump==true) + { + //std::stringstream d_dump_filename_str;//create a stringstream to form the dump filename + //d_dump_filename_str<<"./data/trk_epl_CH_"<d_channel<<"_SAT_"<d_satellite<<".dat"; + //d_dump_filename=d_dump_filename_str.str(); + if (d_dump_file.is_open()==false) + { + try { + d_dump_filename="track_ch"; //base path and name for the tracking log file + d_dump_filename.append(boost::lexical_cast(d_channel)); + d_dump_filename.append(".dat"); + d_dump_file.exceptions ( std::ifstream::failbit | std::ifstream::badbit ); + d_dump_file.open(d_dump_filename.c_str(), std::ios::out | std::ios::binary); + std::cout<<"Tracking dump enabled on channel "<d_satellite+1 << std::endl; } void gps_l1_ca_dll_pll_tracking_cc::update_local_code_refs() { float tcode; - float max_tcode; - unsigned int i=0; - int tmp_index; - float block_correction; + int associated_chip_index; + // unified loop for E, P, L code vectors + for (unsigned int i=0;id_satellite<<" received "; - d_dump_file.open(d_dump_filename.c_str(), std::ios::out | std::ios::binary); - d_enable_tracking=true; - std::cout<<"Tracking start on channel "<d_satellite << std::endl; } /*! Tracking signal processing * Notice that this is a class derived from gr_sync_decimator, so each of the ninput_items has vector_length samples */ -int gps_l1_ca_dll_pll_tracking_cc::work (int noutput_items,gr_vector_const_void_star &input_items, gr_vector_void_star &output_items) { - d_sample_counter+=d_vector_length; //count for the processed samples +int gps_l1_ca_dll_pll_tracking_cc::general_work (int noutput_items, gr_vector_int &ninput_items, + gr_vector_const_void_star &input_items, gr_vector_void_star &output_items) { + d_loops_count++; if (d_enable_tracking==true){ + + if (d_pull_in==true) + { + int samples_offset=ceil(d_absolute_code_phase_samples); + d_code_phase_ms=(d_absolute_code_phase_samples*1000.0)/(float)d_fs_in; + consume_each(samples_offset); //shift input to perform alignement with local replica + d_sample_counter+=samples_offset; //count for the processed samples + d_pull_in=false; + return 1; + } float carr_error; float carr_nco; float code_error; float code_nco; float tmp_E,tmp_P,tmp_L; - const gr_complex* in = (gr_complex*) input_items[0]; + const gr_complex* in = (gr_complex*) input_items[0]; //PRN start block alignement float **out = (float **) &output_items[0]; - // Find the size of a "block" or code period in whole samples - // Update the phasestep based on code freq (variable) and // sampling frequency (fixed) - d_code_phase_step = d_code_freq / d_fs_in; + // code_phase_step_per_sample = T_sample/T_chip + d_code_phase_step = d_code_freq / (float)d_fs_in; //[chips] + // variable code PRN sample block size + d_blksize=ceil((1023.0-d_rem_code_phase) / d_code_phase_step); //[samples] - // original variable block size (MATLAB) - //d_blk_size = ceil(((float)d_code_length-d_rem_code_phase) / d_code_phase_step); - // fixed block size (Javi) - d_blk_size = ceil((float)d_code_length*(1/1023000.0)*(float)d_fs_in); + //float rem_code_phase_samples=d_rem_code_phase/d_code_phase_step; + //d_absolute_code_phase_samples-=floor(rem_code_phase_samples); + + //d_rem_code_phase=d_code_phase_step*(rem_code_phase_samples-floor(rem_code_phase_samples)); this->update_local_code_refs(); this->update_local_carrier(); @@ -287,7 +326,8 @@ int gps_l1_ca_dll_pll_tracking_cc::work (int noutput_items,gr_vector_const_void_ d_P_Q=0.0; d_L_Q=0.0; - for(unsigned int i=0;i30) + { + d_carrier_lock_fail_counter++; + }else{ + d_carrier_lock_fail_counter--; + } + if (d_carrier_lock_fail_counter>200) + { + std::cout<<"Channel "<push(tracking_message); + d_carrier_lock_fail_counter=0; + d_enable_tracking=false; // TODO: check if disabling tranking is consistent with the channel state machine + + } + //std::cout<<"d_carrier_lock_fail_counter"<handle(cmf->GetQueueMessage(200,0)); //send stop to the control_thread - delete cmf; - std::cout<<"stop sent from tracking" << std::endl; + + if (d_sample_counter>round((float)this->d_fs_in*70)){ //stop after some seconds debug only! + d_enable_tracking=false; + std::cout<<"Stop tracking at sample "<handle(cmf->GetQueueMessage(200,0)); //send stop to the control_thread + delete cmf; + std::cout<<"stop sent from tracking"; + } } + }else + { + if (floor(d_sample_counter/d_fs_in)!=d_last_seg) + { + d_last_seg=floor(d_sample_counter/d_fs_in); + std::cout<<"TRK CH "< bits, const bits_slice *slices, int num_of_slices) { unsigned long int value; @@ -117,7 +118,6 @@ unsigned long int gps_navigation_message::read_navigation_unsigned(std::bitset bits, const bits_slice *slices, int num_of_slices) { signed long int value=0; @@ -145,13 +145,7 @@ signed long int gps_navigation_message::read_navigation_signed(std::bitsetrun(); delete control_thread; - google::ShutDownCommandLineFlags(); + //google::ShutDownCommandLineFlags(); }