# End-to-End GPS Legacy Navigation Receiver Using C/A-Code

This example shows how to estimate the global positioning system (GPS) receiver position using a multi-satellite GPS baseband waveform. You use the receiver independent exchange format (RINEX) and an almanac file to model the GPS constellation and generate a multi-satellite baseband waveform. Simulate the satellite scenario to get relative positions of satellites with respect to the modeled receiver. For this satellite scenario, model the Doppler shift, delay, and received signal power. Based on these calculations, impair the generated baseband signal with Doppler shift, delay, and noise. This example shows how to estimate the simulated receiver position from this impaired GPS baseband signal.

### Introduction

In this example, you start with a RINEX file and use `rinexread`

(Navigation Toolbox) to read the file and provide input to `satelliteScenario`

(Satellite Communications Toolbox) to simulate the GPS constellation. Based on this GPS constellation and a given receiver position, you calculate the Doppler shift, delay, and signal path loss from each visible satellite to the receiver. Based on the ephemeris and clock and almanac data from the RINEX and almanac files, generate data bits as per the IS-GPS-200 [1] standard. For the satellites that are visible to the receiver, generate coarse acquisition code (C/A-code) and precision code (P-code). P-code is placed on the in-phase (I) branch and C/A-code is placed on the quadrature-phase (Q) branch of the baseband waveform. The I-branch signal is attenuated by 3 dB as given in IS-GPS-200 [1]. Generate the baseband waveform for all the visible satellites and pass this baseband waveform through the propagation channel as shown in the following figure. Propagation channel characteristics for each satellite signal are unique because the position and velocity of each satellite with respect to the receiver is unique. In this example, you model the characteristics of the propagation channel — namely, the Doppler shift, delay, and scaling signal power based on the propagation path loss — and add thermal noise to the composite signal. Provide this noisy signal as input to the GPS receiver. This example supports storing this noisy signal in a file so that you can test your receiver.

For a detailed description of properties of GPS waveform generation, see GPS Waveform Generation (Satellite Communications Toolbox).

For more details on how to calculate latency and Doppler, see Calculate Latency and Doppler in a Satellite Scenario (Satellite Communications Toolbox).

The complex baseband waveform generated from the propagation channel is processed through a GPS receiver. The following diagram shows the high-level details of such a GPS receiver. First, the initial synchronization module detects visible satellites. It estimates coarse values of the Doppler offset and delay for the visible satellites. For each of the detected satellites, create separate receiver channels (not to be confused with the propagation channel). Each receiver channel performs tracking, bit synchronization, frame synchronization, data decoding, transmitted time estimation, satellite position estimation, and pseudo-ranges computation. For the estimation of the receiver position, estimate the position of the satellites in the sky and distances from the satellites to the receiver. Estimate the position of the satellites from the decoded ephemeris data as provided in Table 20-IV in IS-GPS-200 [1]. Estimate the distances from the satellites to the receiver by estimating the propagation time and multiplying it with the speed of light. For computation of this propagation time at the receiver, compute the difference between receiver time and transmission time. Although accurately computing the receiver time is ideal, GPS receiver cannot compute it. It can estimate only the signal transmission time with high accuracy [2]. Accurate receiver time is not needed for GPS receiver position estimation. At a given receiver time, compute the transmission time of each satellite signal and estimate the pseudo-range for each satellite from this value. The pseudo-ranges are not true ranges and have large errors, but because the errors are the same for all the satellites, you can account for the error when solving for the receiver position and get the accurate receiver position [2]. For more details on receiver position estimation, see GPS Receiver Position Estimation section.

For a detailed explanation of initial synchronization and tracking, see GPS Receiver Acquisition and Tracking Using C/A-Code (Satellite Communications Toolbox).

For more details of decoding GPS data, see GPS Data Decode (Satellite Communications Toolbox).

For more details on estimating receiver position, see Estimate GNSS Receiver Position with Simulated Satellite Constellations (Satellite Communications Toolbox).

### Initialize Parameters

Initialize the parameters that are necessary to configure and run the end-to-end GPS receiver simulation.

Initialize the data duration for which this example must run. Typically, a GPS receiver needs at least 50 seconds of data to estimate the receiver position. To simulate this example for 50 seconds of data would take a lot of wall-clock time based on the available computer resources. For the purposes of this example, set the data duration to 3 seconds. Also, initialize the sampling rate for baseband waveform generation.

simulatedDataDuration = 3; % In seconds samplingRate = 10.23e6; % In Hz

If needed, enable the `WriteWaveformToFile`

property so that you can process this waveform through the receiver of your choice.

`WriteWaveformToFile = false;`

This example processes data in chunks of 1 millisecond. That is, 1 millisecond of data is generated by the waveform generator and that 1 millisecond of data is processed by the receiver. This example works for a step time of one millisecond only.

```
stepTime = 1e-3; % In seconds
numSteps = (simulatedDataDuration/stepTime) + 1;
```

Initialize the parameters required to set up the channel. This example generates a GPS waveform from a simulated GPS constellation. As a starting point to simulate the GPS constellation, read the RINEX file.

rinexFileName = "GODS00USA_R_20211750000_01D_GN.rnx"; almanacFileName = "gpsAlmanac.txt"; rinexdata = rinexread(rinexFileName);

Initialize the receiver position to model the propagation channel based on its location. To estimate the receiver position, model the propagation delay accurately. The delay depends on the receiver position. This example models a stationary receiver and does not support a moving receiver.

rxlat = 39.021; % Latitude in degrees north rxlon = -76.827; % Longitude in degrees east; negative longitude is degrees west rxel = 19; % Elevation in meters

Because the GPS data bit rate is 50 bits per second and each C/A-code block is 1 millisecond, each bit consists of 20 C/A-code blocks. Initialize this parameter.

numCACodeBlocksPerBit = 20;

Initialize the properties required for the receiver. The Receiver must wait for some time before it can start receiving some meaningful data because of large delays modeled in the example. A typical GPS signal has a delay of 60 to 90 milliseconds. So, if receiver waits for 100 milliseconds, then it starts to process signals rather than pure noise.

rxWaitTime = 100; % Milliseconds performInitSync = true; % Initially this must be set to true % Initialize maximum number of tracking channels. Minimum of 4 tracking % channels are needed for proper functioning of the GPS receiver. maxNumTrackingChannels = 8; % Noise bandwidth of each of the tracking loops PLLNoiseBW = 90; % In Hz FLLNoiseBW = 4; % In Hz DLLNoiseBW = 3; % In Hz % Bit synchronization parameters isBitSyncComplete = zeros(maxNumTrackingChannels,1); numBitsForBitSync = 100; numWaitingStepsForBitSync = numCACodeBlocksPerBit*numBitsForBitSync; rxcntr = 1;

Position estimation requires, at minimum, that subframes 2 and 3 are decoded. The length of a subframe is 6 seconds. For this example, decode 48.5 seconds of data to ensure that subframes 2 and 3 are included.

minTimeForPosEst = 48.5; % In seconds minStepsForPosEst = minTimeForPosEst/stepTime; subframeDuration = 6; % In seconds numStepsPerSubframe = subframeDuration/stepTime;

Initialize the physical constants required for the simulation.

c = physconst("LightSpeed"); % Speed of light in m/sec fe = 1575.42e6; % GPS L1 frequency in Hz Dt = 12; % Directivity of the transmit antenna in dBi DtLin = db2pow(Dt); Dr = 4; % Directivity of the receive antenna in dBi DrLin = db2pow(Dr); Pt = 44.8; % Typical transmission power of a GPS satellite in watts k = physconst("boltzmann"); % Boltzmann constant in Joules/Kelvin T = 300; % Room temperature in Kelvin rxBW = 24e6; % Bandwidth in Hz Nr = k*T*rxBW; % Thermal noise power in watts rng default; % Initializing to default random number generation

### Simulation Configuration

In this section, configure the example based on the parameters that are initialized in the Initialize Parameters section.

Set up the satellite scenario based on the RINEX file.

% Initialize satellite scenario sc = satelliteScenario; % Set up the satellites based on the RINEX data sat = satellite(sc,rinexdata,"OrbitPropagator","gps"); rx = groundStation(sc,rxlat,rxlon); % Set up the receiver ac = access(sat,rx); % Calculate access between satellites and the receiver % Get the list of satellites that are considered in satellite scenario from % the RINEX data indices = ones(length(sat),1); for isat = 1:length(sat) ele = orbitalElements(sat(isat)); % Check for match of time of applicability and pseudo-random noise % (PRN) IDs so that data from RINEX file is considered for waveform % generation indices(isat) = find(rinexdata.GPS(:,:).Toe == ele.GPSTimeOfApplicability & ... rinexdata.GPS(:,:).SatelliteID == ele.PRN); end

Generate the navigation data to transmit. First generate a navigation configuration object using the helper function `HelperGPSRINEX2Config`

. This helper function maps the parameters read from the RINEX file to the configuration parameters needed for navigation data generation.

% Generate navigation configuration object navcfg = HelperGPSRINEX2Config(almanacFileName,rinexdata.GPS(indices,:)); [mintow,locmintow] = min([navcfg(:).HOWTOW]); % The time of week (TOW) value is such that it must be 1 plus a multiple of % 5 so that data of subframe 1 is always generated first. mintow = ceil((mintow-1)/5)*5 + 1; % HOWTOW is a counter that contains integer values. Incrementing by a value % of 1 represents 6 seconds of data. The counter increases by a value of 1 % for each new subframe. [navcfg(:).HOWTOW] = deal(mintow); % Set the starting of frames based on mintow firstsubframeID = mod(mintow-1,125) + 1; frameID = ceil(firstsubframeID/5); allFrameIDs = [frameID:25,1:(frameID-1)]; [navcfg(:).FrameIndices] = deal(allFrameIDs); % Generate GPS navigation data numNavBits = 37500; % Full GPS data length is 37500 bits navdata = zeros(numNavBits,length(navcfg)); for isat = 1:length(navcfg) navdata(:,isat) = HelperGPSNAVDataEncode(navcfg(isat)); end

Get the satellite positions and velocities over time for accurate modeling of the Doppler, delay, and power. For more information, see Calculate Latency and Doppler in a Satellite Scenario (Satellite Communications Toolbox).

sc.StartTime = HelperGPSConvertTime(navcfg(locmintow).CEIDataSet.WeekNumber, ... mintow*subframeDuration); sc.StopTime = sc.StartTime + seconds(simulatedDataDuration); sc.SampleTime = stepTime; acstats = accessStatus(ac); % This example runs for a maximum of 2 minutes of data and in that % duration, access does not change. Hence, consider only the access status % of the first sample time. satindices = find(acstats(:,1)); numsat = length(satindices); % Get the states of satellites in the sky over time [initsatpos,satvel] = states(sat,"CoordinateFrame","ecef"); % Shuffle dimensions such that each row corresponds to one satellite data initsatpos = permute(initsatpos,[3,1,2]); satvel = permute(satvel,[3,1,2]);

This figure shows the satellites in the GPS constellation. The green dotted lines indicate access between the receiver and the GPS satellites. Use `play(sc)`

to get this figure.

% Calculate Doppler offset over time for all the visible satellites [az,el] = aer(sat(satindices),rx); [~,satV] = states(sat(satindices),CoordinateFrame="geographic"); dir = cat(1,permute(cosd(el).*cosd(az),[3 2 1]), ... permute(cosd(el).*sind(az),[3 2 1]), ... permute(-sind(el),[3 2 1])); dopV = permute(dot(satV,dir),[3 2 1]); fo = ((c./(c-dopV))*fe).'; % Observed frequency fShift = (fo - fe); % Doppler shift over time for all the visible satellites % Compute the distance between satellites and receiver over time satdist = zeros(numSteps,numsat); for istep = 1:numSteps satdist(istep,:) = pseudoranges([rxlat,rxlon,rxel], ... initsatpos(satindices,:,istep),"RangeAccuracy",0); end delays = satdist/c; % Power at receiver from free space pathloss equation SqrtPr = sqrt(Pt*DtLin*DrLin)*(c./(4*pi*fo.*satdist)); timeinweek = mintow*6; PRNIDs = [navcfg(:).PRNID]; disp("Available satellites - " + num2str(PRNIDs(satindices)))

Available satellites - 5 10 13 15 18 20 23 24 27 29 30

Generate C/A-codes for all the satellites in the GPS constellation. Use the codes for visible satellites while generating the waveform.

cacodes = gnssCACode(PRNIDs,"GPS"); caCodeBlockDuration = 1e-3; % Constant value numcacodeblocks = stepTime/caCodeBlockDuration; % Each C/A-code block is of 1 millisecond duration % To rate match the C/A-code with P-code, repeat each element 10 times cacodesig1 = double(repelem(1-2*repmat(cacodes,numcacodeblocks,1),10,1)); % Rate match the C/A-code with the sampling rate [upfac,downfac] = rat(samplingRate/10.23e6); upcacodesig = repelem(cacodesig1,upfac,1); cacodesig = upcacodesig(1:downfac:end,:);

Because P-code generation is a time-consuming process, and this example does not use the P-code for any receiver operations, use alternating values of `1`

and `0`

instead of the actual P-code. Divide the signal value by $\sqrt{2}$ because the I-branch on which P-code sits is attenuated by 3 dB as per IS-GPS-200 [1].

numSamples = stepTime*samplingRate; pcode = (1-2*repmat([1;0],numSamples/2,1))/sqrt(2);

Initialize the object needed to model the Doppler shift.

pfo = comm.PhaseFrequencyOffset(FrequencyOffsetSource = "Input port", ... SampleRate = samplingRate);

Initialize the object needed to model dynamic delays. The dynamic delay on the satellite signal is modeled in two steps. First, a static delay is modeled which does not change with time. Then, model the variable delay using `dsp.VariableFractionalDelay`

. Initially while modeling the static delay, do not introduce the entire delay so that dynamic delay modeling handles the remaining delay. This dynamic delay is set to 2000 in this example.

% Initialize static delay object dynamicDelayRange = 20000; staticdelay = round(delays(1,:)*samplingRate - dynamicDelayRange); if nnz(staticdelay<0)~=0 staticdelay = zeros(1,numsat); end staticDelayObj = dsp.Delay("Length",staticdelay); % Initialize Variable Fractional Delay object for modeling dynamic delay vfd = dsp.VariableFractionalDelay("InterpolationMethod","Farrow", ... "MaximumDelay",65535);

Initialize the constellation diagram object.

rxconstellation = comm.ConstellationDiagram(1,ShowReferenceConstellation=false, ... Title="Constellation diagram of signal at the output of tracking");

Initialize the object that performs initial synchronization in the GPS receiver.

initialsync = HelperGPSInitialSynchronization; % Initially search all 32 satellites initialsync.PRNIDsToSearch = 1:32; initialsync.SampleRate = samplingRate; initialsync.CenterFrequency = 0 % Baseband signal

initialsync = HelperGPSInitialSynchronization with properties: PRNIDsToSearch: [1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 ... ] CenterFrequency: 0 SampleRate: 10230000 DetectionThresholdFactor: 1.2000

Initialize the baseband file writer.

if WriteWaveformToFile == 1 bbWriter = comm.BasebandFileWriter("gpsBBWaveform.bb",samplingRate,0); end

Initialize the properties that are necessary for the end-to-end simulation chain.

% Properties that store outputs of initial synchronization [doppleroffsets,codephoffsets] = deal(zeros(1,maxNumTrackingChannels)); % Properties required for storing outputs from tracking module [accuPh,accuFqy,accuFqyErr,accuPhErr,accuIntegWave,accuDelay,accuDelayErr] = ... deal(zeros(numSteps,maxNumTrackingChannels)); % Properties to store outputs of bit synchronization and frame % synchronization [maxTransitionLocation,sampleCounter] = deal(zeros(maxNumTrackingChannels,1)); syncidx = zeros(maxNumTrackingChannels,1); % Property to store output of data decoder deccfg = cell(maxNumTrackingChannels,1); % Initialize maximum number of steps for which the simulation chain runs maxSimSteps = 50/stepTime;

### End-to-End Simulation Chain

This section contains three main steps:

Waveform generation

Propagation channel

Receiver signal processing algorithms

As part of the propagation channel, model the Doppler offset, delay, and noise. As part of the receiver, the processes involved include initial synchronization, tracking, bit synchronization, frame synchronization, and data decoding.

tic % Start of simulation numAx = size(initsatpos,2); for istep = 1:numSteps %% Generate waveform bitidx = floor((istep-1)/numCACodeBlocksPerBit)+1; allpcode = repmat(pcode,1,numsat); % Get navigation bit of each satellite at the corresponding step time navbitsig = 1-2*navdata(bitidx,satindices); iqsig = (allpcode + 1j*cacodesig(:,satindices)).*navbitsig; % Implicit expansion %% Propagation channel % Model the channel which models Doppler, delay, path loss, and noise on % the signal % Introduce Doppler to the signal as a frequency offset dopsig = pfo(iqsig,fShift(istep,:)); % Introduce variable fractional delay staticDelayedSignal = staticDelayObj(dopsig); leftoutDelay = delays(istep,:)*samplingRate - staticdelay; % Value must always be positive delayedSig = vfd(staticDelayedSignal,leftoutDelay); % Scale the delayed signal as per received power calculated rmsPow = rms(delayedSig); rmsPow(rmsPow==0) = 1; % To avoid division by zero scaledsig = SqrtPr(istep,:).*delayedSig./rmsPow; % Add up all the signals at the receiver resultsig = sum(scaledsig,2); % Generate noise noisesig = (wgn(numSamples,1,10*log10(Nr)) + 1j*wgn(numSamples,1,10*log10(Nr)))./sqrt(2); % Add constant thermal noise to the composite signal rxwaveform = resultsig + noisesig; % Scale the received signal for having unit power waveform = rxwaveform/rms(rxwaveform); if WriteWaveformToFile == 1 bbWriter(waveform); end timeinweek = timeinweek + stepTime; %% Receiver % Because there are large delays (70 to 80 milliseconds) modeled on the % signal, start the receiver after some time only so as to process % valid signal instead of pure noise. if istep > rxWaitTime if performInitSync == 1 performInitSync = 0; [y,corrval] = initialsync(waveform); % Initial synchronization PRNIDsToSearch = y(y(:,4)==1,1).'; doppleroffsets = y(y(:,4)==1,2).'; codephoffsets = y(y(:,4)==1,3).'; numdetectsat = length(PRNIDsToSearch); if numdetectsat > maxNumTrackingChannels % Set this value to limit the number of tracking channels numdetectsat = maxNumTrackingChannels; end % Perform the required initialization of tracking modules for each % channel and buffers to store the data. disp("The detected satellite PRN IDs: " + num2str(PRNIDsToSearch)) % Plot the correlation plot for the first satellite figure; mesh(-10e3:500:10e3, 0:size(corrval,1)-1, corrval(:,:,PRNIDsToSearch(1))); xlabel("Doppler Offset") ylabel("Code Phase Offset") zlabel("Correlation") msg = ["Correlation Plot for PRN ID: " num2str(PRNIDsToSearch(1))]; title(msg) framesyncbuffer = cell(1,numdetectsat); % Create a cell array, where each element corresponds to a carrier % tracking object. carrierCodeTrack = cell(numdetectsat,1); framesync = cell(numdetectsat,1); % Update properties for each tracking loop for isat = 1:numdetectsat carrierCodeTrack{isat} = HelperGPSCACodeCarrierTracker; carrierCodeTrack{isat}.SampleRate = samplingRate; carrierCodeTrack{isat}.CenterFrequency = 0; carrierCodeTrack{isat}.PLLNoiseBandwidth = PLLNoiseBW; carrierCodeTrack{isat}.FLLNoiseBandwidth = FLLNoiseBW; carrierCodeTrack{isat}.DLLNoiseBandwidth = DLLNoiseBW; carrierCodeTrack{isat}.PLLIntegrationTime = 1; % In milliseconds carrierCodeTrack{isat}.PRNID = PRNIDsToSearch(isat); carrierCodeTrack{isat}.InitialDopplerShift = doppleroffsets(isat); carrierCodeTrack{isat}.InitialCodePhaseOffset = codephoffsets(isat); % Initialize frame synchronization object framesync{isat} = HelperGPSLNAVFrameSynchronizer; end end % Because it would be sufficient to get receiver position after % running the simulation for 50 seconds of data, stop the loop % after executing for 50 seconds of data. In the default run of the % example, only 3 seconds of data is processed and this line is not % needed. This check is when running the example for large data % duration only (for at least 50 seconds of data). if istep > maxSimSteps break; end for isat = 1:numdetectsat % Perform tracking for each satellite [integwave,fqyerr,fqyoffset,pherr,phoffset,derr,dnco] = ... carrierCodeTrack{isat}(waveform); % Accumulate the values to see the results at the end accuFqyErr(rxcntr,isat) = fqyerr; accuFqy(rxcntr,isat) = fqyoffset; accuPhErr(rxcntr,isat) = pherr; accuPh(rxcntr,isat) = phoffset; accuIntegWave(rxcntr,isat) = sum(integwave); accuDelayErr(rxcntr,isat) = derr; accuDelay(rxcntr,isat) = dnco; end % Perform bit synchronization, frame synchronization, and data % decoding if numWaitingStepsForBitSync of receiver steps are % complete. if rxcntr > numWaitingStepsForBitSync % For each detected satellite, perform bit synchronization, % frame synchronization, and data decoding for isat = 1:numdetectsat if ~isBitSyncComplete(isat) maxTransitionLocation(isat) = ... gnssBitSynchronize( ... imag(accuIntegWave(1:numWaitingStepsForBitSync,isat)), ... numCACodeBlocksPerBit); isBitSyncComplete(isat) = 1; sampleCounter(isat) = rxcntr - maxTransitionLocation(isat) + 1; framesyncbuffer{isat} = accuIntegWave( ... maxTransitionLocation(isat):end,isat); else % Perform frame synchronization and data decoding sampleCounter(isat) = sampleCounter(isat) + 1; framesyncbuffer{isat}(sampleCounter(isat)) = accuIntegWave(rxcntr,isat); if mod(sampleCounter(isat),numStepsPerSubframe) == 0 samples = framesyncbuffer{isat}(sampleCounter(isat) - ... numStepsPerSubframe+1:sampleCounter(isat)); sym = mean(reshape(samples,numCACodeBlocksPerBit,[])); bits = imag(sym)<0; [syncidx(isat),rxsubframes,subframeIDs] = framesync{isat}(bits(:)); if ~isempty(rxsubframes) % Then perform data decoding deccfg{isat}.PRNID = PRNIDsToSearch(isat); deccfg{isat} = HelperGPSLNAVDataDecode(rxsubframes,deccfg{isat}); end end end end end if mod(rxcntr,1000) == 0 disp("Processed " + (rxcntr/1000) + " sec of data at the receiver.") rxconstellation(accuIntegWave(rxcntr-999:rxcntr,1)/ ... rms(accuIntegWave(rxcntr-999:rxcntr,1))) end % Update rxcntr rxcntr = rxcntr + 1; end end

The detected satellite PRN IDs: 13 15 24 27 5 29 18 20

Processed 1 sec of data at the receiver. Processed 2 sec of data at the receiver.

if WriteWaveformToFile == 1 release(bbWriter); end e2eTime = toc; disp("End-to-End chain ran for " + e2eTime + " seconds.")

End-to-End chain ran for 74.1389 seconds.

### GPS Receiver Position Estimation

To estimate the receiver position, you need to know the distances from at least 4 satellites to the receiver and the location of each satellite in the sky at the time of transmission. Receiver position calculation from the output of data decoder involves these steps.

Estimate transmission time of visible satellites at given receiver time: Transmission time estimates are the natural measurements of a GPS receiver, not the pseudo-ranges. The transmission time of a signal is a culmination of timing computations performed at initial synchronization, tracking, bit synchronization, frame synchronization, and data decoding.

Compute pseudo-ranges of all satellites: Because the receiver clock is not accurate, you can expect large offsets in the estimated pseudo-ranges. For the same reason that these estimates are not true ranges, they are called pseudo-ranges.

Estimate visible satellites location in sky at given transmission time: From the decoded data, you get the ephemeris data. This ephemeris data is valid for 2-6 hours [2] and to calculate the accurate position of a satellite in the sky from this ephemeris data at a given time, use

`gnssconstellation`

(Navigation Toolbox).Estimate receiver position: Estimate the receiver position using

`receiverposition`

(Navigation Toolbox).

Estimate the transmission time of visible satellites at a given receiver time.

caChipRate = 1.023e6; % In Hz codeOffsetTime = codephoffsets(1:numdetectsat)/caChipRate; rxcntr(rxcntr<=1) = 2; % So that rxcntr-1 is at least 1 trackingOffsetTime = accuDelay(rxcntr-1,1:numdetectsat)/caChipRate; bitsyncTime = (maxTransitionLocation(1:numdetectsat) - 1)*caCodeBlockDuration; framesyncTime = (syncidx(1:numdetectsat)-1)*numCACodeBlocksPerBit*caCodeBlockDuration; % Calculate transmission time from these parameters tt1 = codeOffsetTime(:) - trackingOffsetTime(:) + bitsyncTime + framesyncTime; tt1 = tt1(syncidx~=0);

Estimate the pseudo-ranges.

trtemp = max(tt1); row = (tt1-trtemp)*c; % Include time of week decoded from the received navigation message tow = zeros(maxNumTrackingChannels,1); for isat = 1:maxNumTrackingChannels if isfield(deccfg{isat},"HOWTOW") tow(isat) = (deccfg{isat}.HOWTOW + 1)*6; % Each subframe has a duration of 6 seconds end end tow = tow(syncidx~=0); % Offset the estimated transmission time. This offset is found empirically. timeOffset = 0.74; tt = tt1 + tow + timeOffset;

Estimate the satellite position at the computed transmission time.

deccfg1 = deccfg(syncidx~=0); % Estimate the satellite position only if receiver has processed at least % minStepsForPosEst of data isFullDataDecoded = true(length(tt),1); if rxcntr > minStepsForPosEst [timeofweek,SatelliteID,Delta_n,M0,Eccentricity,sqrtA,Toe,Toc,Cis,Cic,Crs,Crc,Cus,Cuc, ... OMEGA0,i0,omega,OMEGA_DOT,IDOT,GPSWeek] = deal(zeros(4,1)); reqCEIFields = ["WeekNumber","ReferenceTimeOfEphemeris","ReferenceTimeOfClock", ... "SemiMajorAxisLength","MeanMotionDifference","MeanAnomaly","Eccentricity", ... "ArgumentOfPerigee","Inclination","InclinationRate","HarmonicCorrectionTerms", ... "RateOfRightAscension","LongitudeOfAscendingNode"]; for isat = 1:length(tt) timeofweek(isat) = tt(isat); isFullDataDecoded(isat) = (isfield(deccfg1{isat},"PRNID")) && ... (all(isfield(deccfg1{isat}.CEIDataSet,reqCEIFields))); if isFullDataDecoded(isat) SatelliteID(isat) = deccfg1{isat}.PRNID; GPSWeek(isat) = deccfg1{isat}.CEIDataSet.WeekNumber; Toe(isat) = deccfg1{isat}.CEIDataSet.ReferenceTimeOfEphemeris; Toc(isat) = deccfg1{isat}.CEIDataSet.ReferenceTimeOfClock; sqrtA(isat) = sqrt(deccfg1{isat}.CEIDataSet.SemiMajorAxisLength); Delta_n(isat) = deccfg1{isat}.CEIDataSet.MeanMotionDifference*pi; M0(isat) = deccfg1{isat}.CEIDataSet.MeanAnomaly*pi; Eccentricity(isat) = deccfg1{isat}.CEIDataSet.Eccentricity; omega(isat) = deccfg1{isat}.CEIDataSet.ArgumentOfPerigee*pi; i0(isat) = deccfg1{isat}.CEIDataSet.Inclination*pi; IDOT(isat) = deccfg1{isat}.CEIDataSet.InclinationRate*pi; hterms = num2cell(deccfg1{isat}.CEIDataSet.HarmonicCorrectionTerms); [Cis(isat),Cic(isat),Crs(isat),Crc(isat),Cus(isat),Cuc(isat)] = deal(hterms{:}); OMEGA_DOT(isat) = deccfg1{isat}.CEIDataSet.RateOfRightAscension*pi; OMEGA0(isat) = deccfg1{isat}.CEIDataSet.LongitudeOfAscendingNode*pi; end end transmissionTime = HelperGPSConvertTime(GPSWeek(isFullDataDecoded),timeofweek(isFullDataDecoded)); refTime = HelperGPSConvertTime(GPSWeek(:),Toc(:)); rxtimetable = timetable(refTime,SatelliteID,Delta_n,M0,Eccentricity,sqrtA, ... Toe,Cis,Cic,Crs,Crc,Cus,Cuc,OMEGA0,i0,omega, ... OMEGA_DOT,IDOT,GPSWeek); [satpos, satvel] = gnssconstellation(transmissionTime(1),rxtimetable); end

Estimate the receiver position, which requires two inputs — pseudo-ranges and the exact satellite position at the transmission time which is estimated at a given receiver time. You can compute these values only if the example runs for at least 48 seconds of data. Because simulating such data would take a long time, the pseudo-ranges and satellite position properties (obtained by running the example for 50 seconds) are stored in the MAT files. If you run the simulation for less than `minStepsForPosEst`

, then these MAT files are loaded and the receiver position is computed. If you do not, the required properties are taken directly from the simulation chain.

if rxcntr <= minStepsForPosEst % The parameters that are loaded here are valid for the default % configuration of this example load receiverPositionProperties; % When loading the parameters for default configuration, update the % isFullDataDecoded variable with all ones to compute the receiver % position. isFullDataDecoded = true(length(row),1); defaultRINEXFileName = "GODS00USA_R_20211750000_01D_GN.rnx"; defaultRxPos = [39.021, -76.827, 19]; exampleRxPos = [rxlat, rxlon, rxel]; if ~(strcmp(rinexFileName,defaultRINEXFileName) && isequal(defaultRxPos,exampleRxPos)) warning("satcom:EndToEndGPSLNAVReceiverExample:InsufficientData", ... "Estimated receiver position may be different from what you provided" + ... " as the simulation didn't run for entire data." + ... " To get accurate receiver position, run the example" + ... " for at least 50 seconds of navigation data."); end end % Compute the receiver position rxposest = receiverposition(row(isFullDataDecoded),satpos(isFullDataDecoded,:))

`rxposest = `*1×3*
39.0210 -76.8270 18.5258

estRxPosNED = lla2ned(rxposest,[rxlat,rxlon,rxel],'ellipsoid'); distanceError = vecnorm(estRxPosNED) % In meters

distanceError = 4.0582

### Further Exploration

This example shows how to perform GPS receiver processing for the simulated satellite constellation for only 3 seconds. Extend the example to 50 seconds and estimate the receiver position to see how the GPS receiver works.

Use your own RINEX file to configure the example and estimate the receiver position from the waveform generated in this example.

### Appendix

This example uses these data and helper files:

gpsAlmanac.txt — Almanac data file downloaded from Navcen website

HelperGPSAlmanac2Struct.m — Convert text file of almanac to structure

HelperGPSCACodCarrierTracker.m — Carrier frequency and C/A-code phase tracker

HelperGPSCEIConfig.m — Create configuration object for GPS navigation data

HelperGPSConvertTime.m — Convert GPS week and time of week into datetime object and vice-versa

HelperGPSInitialSynchronization.m — GPS Initial synchronization

HelperGPSLNAVDataDecode.m — Decode the LNAV GPS data.

HelperGPSLNAVFrameSynchronizer.m — Perform frame synchronization on the demodulated data.

HelperGPSLNAVWordDecode.m — Decode each word of a subframe.

HelperGPSNAVDataEncode.m — Encode navigation data into bits from data that is in configuration object

HelperGPSNavigationConfig.m — Create configuration object for GPS navigation data

HelperGPSRINEX2Config.m — Convert the properties in RINEX file into navigation configuration object properties

receiverPositionProperties.mat — Contains the properties required for calculation of receiver position

### References

[1] IS-GPS-200, Rev: M. *NAVSTAR GPS Space Segment/Navigation User Interfaces*. May 21, 2021; Code Ident: 66RP1.

[2] Elliott D. Kaplan and C. Hegarty, eds., *Understanding GPS/GNSS: Principles and Applications*, Third edition, GNSS Technology and Applications Series (Boston ; London: Artech House, 2017).