Galileo GNSS Receiver Positioning
This example simulates the Galileo signal generation and reception chain to estimate the location of the receiver. The example creates a multi-satellite scenario using a Galileo Receiver Independent Exchange (RINEX) file, and calculates the Doppler shift and signal latency for a specified receiver position for all the satellites. These effects enable a realistic simulation by impairing the signal. The receiver modules process the superimposed impaired signals of multiple satellites in this order: acquisition, tracking, frame-synchronization, data decoding, and position estimation of the receiver.
Initialize Simulation Parameters
Initialize the relevant simulation parameters.
Check the ShowVisualization box to show the plots in simulation. Enable WriteWaveformToFile to save the generated complex baseband waveform for later use with another receiver. Use the WaveformType setting to specify the band for which to generate the waveform.
ShowVisualizations =false; % Show plots during simulation WriteWaveformToFile =
false; % Save complex baseband waveform to file WaveformType = "E1B";
A Galileo subframe is 30 seconds long and contains 15 nominal pages of 2 seconds each. A Galileo receiver typically needs 35 seconds of transmitted data to collect ephemeris, clock, and other information for position estimation. Therefore, set simulationDuration to at least 35 seconds. Initialize the sampling rate for the baseband waveform generation and other receiver operations. Because the Galileo E1B waveform uses binary offset carrier (BOC) modulation, it requires a higher sampling rate to achieve better position accuracy.
However, if you want faster execution or have limited computational resources, you can set simulationDuration to any value less than 35 seconds since the GalileoRxPosProperties dataset is available. This dataset already includes the necessary data for position estimation, allowing you to conduct quicker simulations.
% Set the simulation duration to a value greater than or equal to 32 seconds. simulationDuration =1; % In seconds minSimulationDuration = 35; % Minimum time to first fix (seconds) % Set the sample rate of the generated waveform. sampleRate =
16000000; % In Hz % Example processes data in chunk of 4ms. stepTime = 4e-3; % In seconds numSteps = floor(simulationDuration/stepTime); % Position update rate. 1 Hz implies 1 position per second. posUpdateFrequency =
1; % In Hz posUpdatePeriod = 1/posUpdateFrequency;
Initialize the receiver position (latitude, longitude, and altitude) to model waveform propagation based on its location. Accurate modeling of propagation delays is crucial for simulating a receiver. This example models a stationary receiver and does not support movement.
rxlat =-15.5725; % In degrees rxlon =
-56.0637; % In degrees rxalt =
100; % In meters minElevationAngle =
10; % In degrees % Initialize the constants used in the example centerFrequency = 1575.42e6; % Needed for Doppler calculation. Units in Hz c = physconst("LightSpeed"); % Speed of light in m/sec Pt = 160; % Typical transmission power of Galileo satellite in watts Dt = 12; % Directivity of the transmit antenna in dBi DtLin = db2pow(Dt); Dr = 4; % Directivity of the receive antenna in dBi DrLin = db2pow(Dr); k = physconst("boltzmann"); % Boltzmann constant in Joules/Kelvin T = 300; % Room temperature in Kelvin Nr = k*T*sampleRate; % Thermal noise power in watts
Because delays are modeled in the propagation channel, a global navigation satellite system (GNSS) receiver receives the actual signal after some delay. Set an appropriate receiver wait time to ensure that receiver captures a meaningful signal instead of noise.
rxWaitSteps = 25; % In steps (one step is 4ms)Configure Simulation
In this section, you specify the configuration parameters to generate the Galileo waveform and corresponding impairments. The simulation uses RINEX files to collect satellite navigation data for Galileo constellation. Firstly, the appropriate data is extracted from the RINEX file. It is used by satelliteScenario to propagate the satellite orbits for a defined simulation duration. The figure outlines the simulation configuration of this section.
.
% Initialize the RINEX and almanac file to use in the example. rinexFileName ="CUIB00BRA_R_20251762130_15M_01S_MN.rnx"; % Specify the RINEX file name almFileName =
"galileoAlmanac.xml"; % Specify the almanac file name fullrinexdata = rinexread(rinexFileName); % Select RINEX rows based on the specified signal type dataSources = fullrinexdata.Galileo.DataSources; if strcmp(WaveformType,"E1B") dataidx = logical(mod(dataSources,2)); % When bit 1 is set for dataSources, it indicates I/NAV E1B data elseif strcmp(WaveformType,"E5a") %#ok<UNRCH> dataidx = bitget(dataSources,2)&1; % When bit 2 is set for dataSources, it indicates F/NAV E5a data else dataidx = bitget(dataSources,3)&1; % When bit 3 is set for dataSources, it indicates I/NAV E5b data end fullSigTypeData = struct; fullSigTypeData.Galileo = fullrinexdata.Galileo(dataidx,:); if isempty(fullSigTypeData) error("RINEX has no "+WaveformType+" data."); end % Get the start time for simulation from the RINEX file. inittow = fullSigTypeData.Galileo.TransmissionTime(1); startTime = HelperGNSSConvertTime(fullSigTypeData.Galileo.GALWeek(1),inittow); % Setup satellite scenario sc = satelliteScenario; sc.StartTime = startTime; sc.SampleTime = stepTime; sc.StopTime = startTime + seconds(simulationDuration); fullSigTypeData.Galileo.Time.TimeZone = "UTC"; % Select valid RINEX rows based on the specified simulation start time. % Only select the data within 30 minutes difference from the start time. validEphemerisIdx = abs(startTime - fullSigTypeData.Galileo.Time) <= minutes(30); validfullSigTypeData.Galileo = sortrows(fullSigTypeData.Galileo(validEphemerisIdx,:),"IODnav"); % Use same issue of data navigation (IODnav) stamp to ensure coherent navigation data across satellites for the selected interval. IODnavCounts = histcounts(validfullSigTypeData.Galileo.IODnav); [maxIODnavCounts,idxmax] = max(IODnavCounts); startidx = sum(IODnavCounts(1:max(1, idxmax - 1))) + 1; endidx = startidx + maxIODnavCounts - 1; rinexdata.Galileo=sortrows(validfullSigTypeData.Galileo(startidx:endidx,:),"SatelliteID"); % Create instances of satellite and stationary receiver sat = satellite(sc,rinexdata); rx = groundStation(sc,rxlat,rxlon,"Altitude",rxalt,"MaskElevationAngle",minElevationAngle); % Initialize the Galileo navigation data class dataobj = HelperGalileoNavigationData(rinexdata,SignalType=WaveformType); setnavdata(dataobj,almFileName,false); % Initialize the almanac data without updating the ephemeris from the almanac file dataobj.DisableWarnings = true; % Disables all warnings from the data generator object % Calculate the visible satellites, Doppler shift, latency, and SNR for all % visible satellites latencyVal = latency(sat,rx); visiblesatindices = find(~isnan(latencyVal(:,1))).'; numSat = length(visiblesatindices); satelliteIDs = sort(dataobj.Ephemeris.SatelliteID(visiblesatindices)); dopplerVal = dopplershift(sat(visiblesatindices),rx,Frequency=centerFrequency); latencyVal = latency(sat(visiblesatindices),rx); Pr = Pt*DtLin*DrLin./ ... ((4*pi*(centerFrequency + dopplerVal).*latencyVal).^2); snrVal = 10*log10(Pr/Nr); disp("Available Galileo satellite IDs: " + num2str(satelliteIDs.'))
Available Galileo satellite IDs: 18 19 21 23 26 27 33
Generate Galileo Navigation Data and Waveform
Initialize the System objects™ and parameters to generate the Galileo waveform. Galileo generates navigation data based on the number of pages that can be transmitted during the specified simulation duration. The example extracts this data from the RINEX file using the Galileo navigation data class.
% Initialize the Galileo waveform generator object gengalwave = galileoWaveformGenerator(SignalType=WaveformType,SampleRate=sampleRate); gengalwave.SVID = satelliteIDs; gengalwave.InitialTime = inittow; % Generate navigation bits for the specified simulation duration numNavDataBits = simulationDuration/gengalwave.BitDuration; numpages = ceil(5*numNavDataBits/250); % The factor of 5 is because INAV will have 5 times more bits in the same duration inavbits = navdata2bits(dataobj,satelliteIDs,inittow,numpages); navbits = inavbits(1:numNavDataBits,:); % Optional waveform writer if WriteWaveformToFile == 1 bbWriter = comm.BasebandFileWriter("BBWaveform.bb",sampleRate,0); end
Transmit-Receive Loop
Initialize the System objects™ and parameters to implement the receiver as a finite state machine (FSM).
% Initialize the receiver operational parameters rxistep = 1; prevwave = 0; poscntr = 1; chipRate = gengalwave.ChipRate; numSymPerPagePart = 250; % Number of symbols in a Galileo INAV page part numStepPerSym = stepTime/gengalwave.BitDuration; % Number of receiver steps to process one symbol numStepsPerPagePart = numSymPerPagePart*numStepPerSym; % Number of receiver steps to process one page part isframesynced = 0; % Frame synchronization flag bit % Initialize maximum number of tracking channels maxNumTrackingChannels = 8; % Noise bandwidth of each of the tracking loops PLLNoiseBW = 40; % In Hz FLLNoiseBW = 4; % In Hz DLLNoiseBW = 1.5; % In Hz % Galileo navigation data class for handling the received data rxDataObj = HelperGalileoNavigationData(SignalType=WaveformType); % Initialize the signal acquisition and signal tracking objects. sigAcquisition = HelperGNSSSignalAcquirer; sigAcquisition.SampleRate = sampleRate; sigAcquisition.GNSSSignalType = "Galileo " + WaveformType; sigAcquisition.FrequencyResolution = 125; sigTracker = HelperGNSSSignalTracker(GNSSSignalType=sigAcquisition.GNSSSignalType, ... SampleRate=sampleRate,PLLNoiseBandwidth=PLLNoiseBW,FLLNoiseBandwidth=FLLNoiseBW,DLLNoiseBandwidth=DLLNoiseBW); trackInfo = struct("PhaseError",[],"PhaseEstimate",[],"FrequencyError",[],"FrequencyEstimate",[], ... "DelayError",[],"DelayEstimate",[]); % Initialize the propagation channel object with the waveform sample rate. gnssChannel = HelperGNSSChannel(RandomStream="mt19937ar with seed",SampleRate=sampleRate); if ShowVisualizations % Object to show the received samples constellation rxconstellation = comm.ConstellationDiagram(1,ShowReferenceConstellation=false, ... Title="Constellation diagram of signal at the " + ... "output of tracking"); %#ok<UNRCH> end % Set the first state of the receiver chain as acquisition of the satellites. nxtState = "acquisition";
The transmitted signal passes through a GNSS channel that introduces impairments such as Doppler shift, propagation delay, and thermal noise. Each impaired signal chunk is then processed by the Galileo receiver, which uses a FSM design. The FSM progresses through five stages: acquisition, tracking, frame synchronization, data decoding, and position estimation. Each state refines synchronization or extracts navigation data, and transitions occur based on predefined conditions that ensure the receiver has sufficient information to move to the next stage. The diagram shows the I/NAV message structure and provides a basis for defining state transition conditions at the level of a page or a page part. In Galileo, a page part represents the smallest message unit and contains 250 symbols.
Finite State Machine Architecture
The receiver FSM consists of five distinct states that represent the sequential processing stages of GNSS signal reception. The diagram outlines the workflow of the receiver FSM.
The following information provides greater detail about the function of each state and its transition logic.
State 1: Acquisition
Perform initial synchronization by detecting visible satellites and obtaining coarse estimates of code phase and Doppler frequency offsets. Process both the previous waveform chunk and the current received waveform to ensure correct correlation peak detection.
Transition Logic:
To Tracking — If four or more satellites are successfully detected, the receiver has sufficient information to proceed. The receiver passes the initial estimates to the tracking module and initializes counters and storage arrays.
To Exit — If fewer than three satellites are detected, acquisition has failed. The simulation terminates as position estimation requires at least four satellites (three for 2D position plus one for time).
State 2: Tracking
Maintain fine synchronization with detected satellites and extract navigation symbols from the received signal. The tracker module refines the code phase and carrier frequency estimates for all detected satellites in parallel.
Transition Logic:
To Frame Sync — When at least 3 page part (3 × number steps per page) are collected from all satellites and frame synchronization is not yet performed.
To Data Decode — When at least one complete page (500 symbols) has been collected and frame synchronization is completed.
To Position Estimate — When the time-to-first-fix (TTFF) threshold has been reached and it is time for a position update (based on position update frequency).
Remain in Tracking — Otherwise, continue accumulating samples.
State 3: Frame Synchronization
Determine the exact start of navigation data page and detect bit inversions. Frame synchronization is successful when the parts of the received samples match sync-code defined in the standard.
Transition Logic:
To Data Decode — When at least 500 symbols (one complete page) are available after frame synchronization.
To Tracking — If insufficient symbols are available, return to tracking to collect more data.
State 4: Data Decode
Extract navigation data such as ephemeris, clock corrections, almanac, and more from the demodulated symbols. The data is decoded from chunks of 500 symbols, corresponding to one page, that is, an even and an odd page.
Transition Logic:
To Position Estimate — When the TTFF threshold has been reached and it is time for a position update
To Tracking — Otherwise, continue collecting and decoding more data.
State 5: Position Estimate
Compute the geographic position of the receiver using pseudorange measurements and satellite positions. Estimate the pseudorange by combining three timing observables obtained during the acquisition, tracking, and frame synchronization stages. The pseudorange represents the time a satellite signal takes to reach the receiver, converted into distance.
Transition Logic:
To Tracking — Always return to the tracking state to process incoming signals and update the position periodically.
tic for istep = 1:numSteps % Transmit Galileo waveform for this step iwave = gengalwave(navbits(istep,:)); % Introduce propagation channel effects to the transmitted signal gnssChannel.SignalToNoiseRatio = snrVal(:,istep); gnssChannel.SignalDelay = latencyVal(:,istep)'; gnssChannel.FrequencyOffset = dopplerVal(:,istep)'; irxwave = gnssChannel(iwave); % Optionally write the waveform to a file if WriteWaveformToFile == 1 bbWriter(irxwave) end if strcmp(nxtState,"exit") break; end % Receiver if istep > rxWaitSteps while true switch(nxtState) case "acquisition" % Initial synchronization [acqd,corrval] = sigAcquisition([prevwave; irxwave],1:50); acqIdx = acqd.IsDetected == 1; numRxSat = min(maxNumTrackingChannels,nnz(acqIdx)); detectedSatTable = sortrows(acqd(1:numRxSat,:)); svidToSearch = detectedSatTable.PRNID; disp("The detected satellite PRN IDs: " + num2str(svidToSearch')) if sum(acqIdx) >= 4 % If four or more satellites are detected sigTracker.InitialFrequencyOffset = detectedSatTable.FrequencyOffset; sigTracker.InitialCodePhaseOffset = detectedSatTable.CodePhaseOffset; sigTracker.PRNID = svidToSearch; invertedBitsSatIdx = zeros(numRxSat,1); refTOW = -ones(numRxSat,1); refrxistep = -ones(numRxSat,1); allrxposest = zeros((simulationDuration-minSimulationDuration)/posUpdatePeriod,3); allposerrors = zeros((simulationDuration-minSimulationDuration)/posUpdatePeriod,1); % To store the frame synchronization index satSFstartIdx = zeros(numRxSat,1); towBitIdx = zeros(numRxSat,1); nxtState = "tracking"; % Property for storing the output of the tracking module. trackedWave = zeros(numSteps - rxWaitSteps,numRxSat); % Property for counting the received samples from each satellite. sampleCntr = zeros(1,numRxSat); proccessedPageCntr = zeros(numRxSat,1); % Estimate code phase offset time. One of the % measurements for pseudorange estimation. codeOffsetTime = sigTracker.InitialCodePhaseOffset/chipRate; if ShowVisualizations % Correlation plot for first detected satellite figure; %#ok<UNRCH> mesh(-5e3:125:5e3,0:size(corrval,1)-1, corrval(:,:,1)) xlabel("Doppler Offset") ylabel("Code Phase Offset") zlabel("Correlation") msg = ["Correlation Plot for PRN ID: " num2str(svidToSearch(1))]; title(msg) end else % If acquisition fails because fewer than three % satellites are detected, the simulation ends. nxtState = "exit"; break end case "tracking" [trackedWave(rxistep,:),trackInfo(rxistep)] = sigTracker(irxwave); % Store the tracked samples for isat = 1:numRxSat sampleCntr(isat) = sampleCntr(isat) + 1; rxSamples{isat}(sampleCntr(isat)) = trackedWave(rxistep,isat); end rxistep = rxistep + 1; % Next state decision logic if min(sampleCntr) == 3*numStepsPerPagePart && ~isframesynced nxtState = "frame-sync"; elseif ~mod(min(sampleCntr),2*numSymPerPagePart) && isframesynced nxtState = "data-decode"; elseif istep*stepTime >= minSimulationDuration if ~(mod(istep*stepTime,posUpdatePeriod)) nxtState = "pos-estimate"; else nxtState = "tracking"; break; end else nxtState = "tracking"; break; end case "frame-sync" for isat = 1:numRxSat % Obtain symbols from the received IQ samples. irxsyms = imag(rxSamples{isat}.') < 0; % Synchronize the frame, and obtain the subframe start index [SFstartIdx,irxsyms,isInverted] = HelperGalileoINAVFrameSynchronize(irxsyms); isFirstPageEven = checkFirstPageIsEven(rxDataObj,irxsyms(1:250)); invertedBitsSatIdx(isat) = isInverted; rxSamples{isat} = rxSamples{isat}(SFstartIdx:end); if ~isFirstPageEven rxSamples{isat} = rxSamples{isat}(251:end); end sampleCntr(isat) = length(rxSamples{isat}); satSFstartIdx(isat) = SFstartIdx; end isframesynced = all(satSFstartIdx); % Observable for pseudorange estimation framesyncTime = (satSFstartIdx)*stepTime; if min(sampleCntr)>=500 nxtState = "data-decode"; else nxtState = "tracking"; break; end case "data-decode" for isat = 1:numRxSat startIdx = proccessedPageCntr(isat)*500+1; endIdx = proccessedPageCntr(isat)*500+500; irxsyms = imag(rxSamples{isat}(startIdx:endIdx).') < 0; proccessedPageCntr(isat) = proccessedPageCntr(isat) + 1; % Decode the symbols towBitIdx(isat) = bits2navdata(rxDataObj,xor(invertedBitsSatIdx(isat),irxsyms),svidToSearch(isat)); if towBitIdx(isat) & ~isempty(rxDataObj.Ephemeris.TimeOfWeek) refTOW(isat) = rxDataObj.Ephemeris.TimeOfWeek(isat); refrxistep(isat) = (rxistep - 1) - (500 - towBitIdx(isat)) - (sampleCntr(isat) - min(sampleCntr)); end end if istep*stepTime >= minSimulationDuration if ~(mod(istep*stepTime,posUpdatePeriod)) nxtState = "pos-estimate"; else nxtState = "tracking"; break; end else nxtState = "tracking"; break; end case "pos-estimate" % Estimate Pseudorange trackingOffsetTime = trackInfo(rxistep - 1).DelayEstimate/chipRate; % Calculate propagation delay from time measurements delayEst = (codeOffsetTime + framesyncTime - trackingOffsetTime'); rho = delayEst*c; % Structure the decoded navigation data as RINEX-style tables eph = rxDataObj.Ephemeris; galweek = rxDataObj.Ephemeris.WeekNumber; tow = rxDataObj.Ephemeris.TimeOfWeek; clock = rxDataObj.Clock; Time = HelperGNSSConvertTime(galweek,clock.Toc); % Construct time-table of the received clock and ephemeris data. rxtimetable = timetable(Time,eph.SatelliteID,clock.SVClockBias,clock.SVClockDrift,clock.SVClockDriftRate,eph.IODnav, ... eph.Crs,eph.MeanMotionDifference*pi,eph.MeanAnomaly*pi,eph.Cuc,eph.Eccentricity,eph.Cus,sqrt(eph.SemiMajorAxis),eph.Toe,eph.Cic, ... eph.LongitudeOfAscendingNode*pi,eph.Cis,eph.Inclination*pi,eph.Crc,eph.ArgumentOfPerigee*pi,eph.RateOfRightAscension*pi,eph.InclinationRate*pi,galweek); rxtimetable.Properties.VariableNames = [rinexdata.Galileo.Properties.VariableNames(1:21) {'GALWeek'}]; transmissionTime = HelperGNSSConvertTime(galweek,tow)+seconds((rxistep - 1 - refrxistep)*stepTime); % Estimate satellite positions using common reception % time algorithm satpos = zeros(numRxSat,3); for isat = 1:numRxSat [satpos(isat,:),~] = gnssconstellation(transmissionTime(isat),rxtimetable(isat,:)); end % Estimate receiver position [rxposest,~,hdop,vdop,info] = receiverposition(rho,satpos); allrxposest(poscntr,:) = rxposest; estRxPosNED = lla2ned(rxposest,[rxlat rxlon rxalt],"ellipsoid"); distanceError = vecnorm(estRxPosNED); % In meters allposerrors(poscntr) = distanceError; % Display the error of the estimated position disp("Position estimation error is " + distanceError + " meters") poscntr = poscntr + 1; if hdop > 20 warning("Dilution of Precision (DOP) ratings are poor. The position " + ... "estimation error can be high.") end nxtState = "tracking"; break; end end end prevwave = irxwave; if ~mod(istep,250) disp("Processed " + (istep/250) + " sec of data at the receiver.") end end
The detected satellite PRN IDs: 18 19 21 23 26 27 33
Processed 1 sec of data at the receiver.
toc
Elapsed time is 47.013303 seconds.
% When simulation duration is lesser than minimum simulation duration or % decoded data is insufficient if length(rxDataObj.Ephemeris.SatelliteID)<4 load("GalileoRxPosProperties.mat") disp("Position estimation error is " + distanceError + " meters") end
Position estimation error is 34.2454 meters
if WriteWaveformToFile % Release the waveform writer object release(bbWriter) %#ok<UNRCH> end
Supporting Files
This example uses these data and helper files:
HelperGNSSConvertTime.m — Convert GNSS week and time of week into datetime object and the other way around.
HelperGNSSChannel.m — Provide a GNSS propagation channel with Doppler shift, signal delay, and random noise
HelperGalileoINAVFrameSynchronize.m — Synchronize the received symbols to find the start of the page
HelperGalileoNavigationData.m — Generate Galileo F/NAV or I/NAV data. Convert Galileo Navigation message to bits and vice versa
HelperGNSSSignalAcquirer.m — Acquire visible satellites and estimate coarse frequency and phase offsets
HelperGNSSSignalTracker.m — Carrier frequency and code phase tracker
GalileoRxPosProperties.mat — Contains the properties required for calculation of the receiver position
galileoAlmanac.xml — Sample Galileo almanac file downloaded from gsc-europa on November 29 2024
References
[1] European GNSS Service Centre (GSC). Galileo Open Service Signal-In-Space Interface Control Document. OS SIS ICD v2.1. GSC, November 2023. https://www.gsc-europa.eu/sites/default/files/sites/all/files/Galileo_OS_SIS_ICD_v2.1.pdf.
See Also
galileoWaveformGenerator | galileoCodes
