Main Content

Esta página se ha traducido mediante traducción automática. Haga clic aquí para ver la última versión en inglés.

Receptor de navegación GPS heredado de extremo a extremo mediante código C/A

Este ejemplo utiliza:

Abrir script en vivo

Este ejemplo muestra cómo estimar la posición del receptor del sistema de posicionamiento global (GPS) utilizando una forma de onda de banda base GPS multisatélite. Se utiliza el formato de intercambio independiente del receptor (RINEX) y un archivo de almanaque para modelar la constelación GPS y generar una forma de onda de banda base multisatélite. Simule el escenario del satélite para obtener posiciones relativas de los satélites con respecto al receptor modelado. Para este escenario de satélite, modele el desplazamiento Doppler, el retraso y la potencia de la señal recibida. Según estos cálculos, deteriore la señal de banda base generada con desplazamiento Doppler, retraso y ruido. Este ejemplo muestra cómo estimar la posición del receptor simulada a partir de esta señal de banda base de GPS deteriorada.

Introducción

En este ejemplo, comienza con un archivo RINEX y utiliza rinexread (Navigation Toolbox) para leer el archivo y proporcionar información a satelliteScenario ( Communications Toolbox por satélite) para simular el Constelación GPS. Con base en esta constelación de GPS y una posición determinada del receptor, se calcula el desplazamiento Doppler, el retraso y la pérdida de ruta de la señal desde cada satélite visible hasta el receptor. Con base en las efemérides y los datos de reloj y almanaque de los archivos RINEX y almanaque, genere bits de datos según el estándar IS-GPS-200 [1] . Para los satélites que son visibles para el receptor, genere un código de adquisición aproximado (código C/A) y un código de precisión (código P). El código P se coloca en la rama en fase (I) y el código C/A se coloca en la rama de fase en cuadratura (Q) de la forma de onda de banda base. La señal de la rama I se atenúa en 3 dB como se indica en IS-GPS-200 [1]. Genere la forma de onda de banda base para todos los satélites visibles y pase esta forma de onda de banda base a través del canal de propagación como se muestra en la siguiente figura. Las características del canal de propagación para cada señal de satélite son únicas porque la posición y velocidad de cada satélite con respecto al receptor es única. En este ejemplo, modelará las características del canal de propagación (es decir, el desplazamiento Doppler, el retardo y el escalamiento de la potencia de la señal en función de la pérdida de la ruta de propagación) y agregará ruido térmico a la señal compuesta. Proporcione esta señal ruidosa como entrada al receptor GPS. Este ejemplo admite el almacenamiento de esta señal ruidosa en un archivo para que pueda probar su receptor.

La compleja forma de onda de banda base generada a partir del canal de propagación se procesa a través de un receptor GPS. El siguiente diagrama muestra los detalles de alto nivel de dicho receptor GPS. Primero, el módulo de sincronización inicial detecta satélites visibles. Estima valores aproximados del desplazamiento y retraso Doppler para los satélites visibles. Para cada uno de los satélites detectados, cree canales de recepción separados (que no debe confundirse con el canal de propagación). Cada canal receptor realiza seguimiento, sincronización de bits, sincronización de tramas, decodificación de datos, estimación del tiempo de transmisión, estimación de la posición del satélite y cálculo de pseudodistancias. Para estimar la posición del receptor, estime la posición de los satélites en el cielo y las distancias de los satélites al receptor. Calcule la posición de los satélites a partir de los datos de efemérides decodificados como se proporciona en la Tabla 20-IV en IS-GPS-200 [1]. Calcule las distancias de los satélites al receptor estimando el tiempo de propagación y multiplicándolo por la velocidad de la luz. Para calcular este tiempo de propagación en el receptor, calcule la diferencia entre el tiempo del receptor y el tiempo de transmisión. Aunque lo ideal es calcular con precisión la hora del receptor, el receptor GPS no puede calcularla. Solo puede estimar el tiempo de transmisión de la señal con alta precisión [2]. No se necesita una hora exacta del receptor para estimar la posición del receptor GPS. En un momento dado del receptor, calcule el tiempo de transmisión de cada señal de satélite y estime el pseudorango para cada satélite a partir de este valor. Los pseudorangos no son rangos verdaderos y tienen grandes errores, pero debido a que los errores son los mismos para todos los satélites, puede tener en cuenta el error al resolver la posición del receptor y obtener la posición precisa del receptor [2]. Para obtener más detalles sobre la estimación de la posición del receptor, consulte la sección GPS Receiver Position Estimation .

Inicializar parámetros

Inicialice los parámetros necesarios para configurar y ejecutar la simulación del receptor GPS de un extremo a otro.

Inicialice la duración de los datos durante la cual se debe ejecutar este ejemplo. Normalmente, un receptor GPS necesita al menos 50 segundos de datos para estimar la posición del receptor. Simular este ejemplo con 50 segundos de datos requeriría mucho tiempo según los recursos informáticos disponibles. Para los fines de este ejemplo, establezca la duración de los datos en 3 segundos. Además, inicialice la frecuencia de muestreo para la generación de formas de onda de banda base.

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

Si es necesario, habilite la propiedad WriteWaveformToFile para que pueda procesar esta forma de onda a través del receptor de su elección.

WriteWaveformToFile = false;

Este ejemplo procesa datos en fragmentos de 1 milisegundo. Es decir, el generador de formas de onda genera 1 milisegundo de datos y el receptor procesa ese milisegundo de datos. Este ejemplo funciona solo para un tiempo de paso de un milisegundo.

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

Inicialice los parámetros necesarios para configurar el canal. Este ejemplo genera una forma de onda GPS a partir de una constelación GPS simulada. Como punto de partida para simular la constelación GPS, lea el archivo RINEX.

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

Inicialice la posición del receptor para modelar el canal de propagación en función de su ubicación. Para estimar la posición del receptor, modele con precisión el retardo de propagación. El retraso depende de la posición del receptor. Este ejemplo modela un receptor estacionario y no admite un receptor en movimiento.

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

Debido a que la velocidad de bits de datos del GPS es de 50 bits por segundo y cada bloque de código C/A es de 1 milisegundo, cada bit consta de 20 bloques de código C/A. Inicialice este parámetro.

numCACodeBlocksPerBit = 20;

Inicialice las propiedades necesarias para el receptor. El receptor debe esperar algún tiempo antes de poder comenzar a recibir algunos datos significativos debido a los grandes retrasos modelados en el ejemplo. Una señal GPS típica tiene un retraso de 60 a 90 milisegundos. Entonces, si el receptor espera 100 milisegundos, comienza a procesar señales en lugar de ruido puro.

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;

La estimación de la posición requiere, como mínimo, que las subtramas 2 y 3 estén decodificadas. La duración de una subtrama es de 6 segundos. Para este ejemplo, decodifique 48,5 segundos de datos para garantizar que se incluyan las subtramas 2 y 3.

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

Inicialice las constantes físicas necesarias para la simulación.

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

Configuración de simulación

En esta sección, configure el ejemplo según los parámetros que se inicializan en la sección Initialize Parameters .

Configure el escenario satelital basado en el archivo RINEX.

% 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

Generar los datos de navegación a transmitir. Primero genere un objeto de configuración de navegación usando la función auxiliar HelperGPSRINEX2Config. Esta función auxiliar asigna los parámetros leídos del archivo RINEX a los parámetros de configuración necesarios para la generación de datos de navegación.

% 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

Obtenga las posiciones y velocidades de los satélites a lo largo del tiempo para modelar con precisión el Doppler, el retraso y la potencia. Para obtener más información, consulte Calculate Latency and Doppler in a Satellite Scenario (Satellite Communications Toolbox). 

sc.StartTime = HelperGPSConvertTime(navcfg(locmintow).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);

Esta figura muestra los satélites de la constelación GPS. Las líneas de puntos verdes indican el acceso entre el receptor y los satélites GPS. Utilice play(sc) para obtener esta cifra.

% Calculate Doppler shift over time for all the visible satellites
fShift = dopplershift(sat(satindices),rx,Frequency=fe);

% Get the states of satellites in the sky over time. Also, shuffle the
% dimensions such that each row corresponds to one satellite data
initsatpos = permute(states(sat,"CoordinateFrame","ecef"),[3 1 2]);

% Compute the distance between satellites and receiver over time
satdist = zeros(numsat,numSteps);
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)*(1./(4*pi*(fe+fShift).*delays));
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

Genere códigos C/A para todos los satélites de la constelación GPS. Utilice los códigos de los satélites visibles mientras genera la forma de onda.

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);
idx = 1:downfac:upfac*size(cacodesig1,1);
cacodesig = cacodesig1(ceil(idx/upfac),:);

Debido a que la generación de código P es un proceso que requiere mucho tiempo y este ejemplo no utiliza el código P para ninguna operación del receptor, use valores alternos de 1 y 0 en lugar del código P real. Divida el valor de la señal por2porque la rama I en la que se encuentra el código P está atenuada en 3 dB según IS-GPS-200 [1].

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

Inicialice el objeto necesario para modelar el desplazamiento Doppler.

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

Inicialice el objeto necesario para modelar retrasos dinámicos. El retardo dinámico de la señal del satélite se modela en dos pasos. Primero, se modela un retraso estático que no cambia con el tiempo. Luego, modele el retraso variable usando dsp.VariableFractionalDelay. Inicialmente, mientras modela el retraso estático, no introduzca todo el retraso para que el modelado de retraso dinámico maneje el retraso restante.

% Initialize static delay object
dynamicDelayRange = 20000;
staticdelay = round(delays(:,1)*samplingRate - dynamicDelayRange);
if nnz(staticdelay<0)~=0
    staticdelay = zeros(numsat,1);
end
staticDelayObj = dsp.Delay("Length",staticdelay);

% Initialize Variable Fractional Delay object for modeling dynamic delay
vfd = dsp.VariableFractionalDelay("InterpolationMethod","Farrow", ...
    "MaximumDelay",65535);

Inicialice el objeto del diagrama de constelación.

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

Inicialice el objeto que realiza la sincronización inicial en el receptor GPS.

initialsync = gnssSignalAcquirer;
initialsync.SampleRate = samplingRate;
initialsync.IntermediateFrequency = 0        % Baseband signal
initialsync = 
  gnssSignalAcquirer with properties:

              GNSSSignalType: "GPS C/A"
                  SampleRate: 10230000
       IntermediateFrequency: 0
              FrequencyRange: [-10000 10000]
         FrequencyResolution: 500
    DetectionThresholdFactor: 1.9000

Inicialice el escritor de archivos de banda base.

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

Inicialice las propiedades necesarias para la cadena de simulación de un extremo a otro.

% 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;

Cadena de simulación de extremo a extremo

Esta sección contiene tres pasos principales:

  1. Generación de forma de onda

  2. Canal de propagación

  3. Algoritmos de procesamiento de señales del receptor.

Como parte del canal de propagación, modele el desplazamiento, el retraso y el ruido Doppler. Como parte del receptor, los procesos involucrados incluyen sincronización inicial, seguimiento, sincronización de bits, sincronización de gráficas y decodificación de datos.

tic % Start of simulation
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,1:32); % Initial synchronization

            PRNIDsToSearch = y(y(:,4).IsDetected==1,1).PRNID.';
            doppleroffsets = y(y(:,4).IsDetected==1,2).FrequencyOffset;
            codephoffsets = y(y(:,4).IsDetected==1,3).CodePhaseOffset;

            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  24  18  27   5  20  23  15  29  10  30

Figure contains an axes object. The axes object with title Correlation Plot for PRN ID: 13, xlabel Doppler Offset, ylabel Code Phase Offset contains an object of type surface.

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 85.4106 seconds.

Estimación de la posición del receptor GPS

Para estimar la posición del receptor, necesita conocer las distancias de al menos 4 satélites al receptor y la ubicación de cada satélite en el cielo en el momento de la transmisión. El cálculo de la posición del receptor a partir de la salida del decodificador de datos implica estos pasos.

  1. Estimación del tiempo de transmisión de los satélites visibles a una hora determinada del receptor: Las estimaciones del tiempo de transmisión son medidas naturales de un receptor GPS, no pseudodistancias. El tiempo de transmisión de una señal es la culminación de los cálculos de temporización realizados en la sincronización inicial, el seguimiento, la sincronización de bits, la sincronización de gráficas y la decodificación de datos.

  2. Calcular pseudorangos de todos los satélites: Debido a que el reloj del receptor no es preciso, se pueden esperar grandes desfases en los pseudorangos estimados. Por la misma razón que estas estimaciones no son rangos verdaderos, se denominan pseudorangos.

  3. Estimar la ubicación de los satélites visibles en el cielo en un momento de transmisión determinado: De los datos decodificados, se obtienen los datos de efemérides. Estos datos de efemérides son válidos por 2-6 horas [2] y para calcular la posición exacta de un satélite en el cielo a partir de estos datos de efemérides en un momento dado, use gnssconstellation (Navigation Toolbox).

  4. Estimar posición del receptor: Calcule la posición del receptor usando receiverposition (Navigation Toolbox).

Estimar el tiempo de transmisión de los satélites visibles en un momento dado del receptor.

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);

Estimar los pseudorangos.

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;

Estime la posición del satélite en el tiempo de transmisión calculado.

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},reqCEIFields)));

        if isFullDataDecoded(isat)
            SatelliteID(isat) = deccfg1{isat}.PRNID;
            GPSWeek(isat) = deccfg1{isat}.WeekNumber;
            Toe(isat) = deccfg1{isat}.ReferenceTimeOfEphemeris;
            Toc(isat) = deccfg1{isat}.ReferenceTimeOfClock;
            sqrtA(isat) = sqrt(deccfg1{isat}.SemiMajorAxisLength);
            Delta_n(isat) = deccfg1{isat}.MeanMotionDifference*pi;
            M0(isat) = deccfg1{isat}.MeanAnomaly*pi;
            Eccentricity(isat) = deccfg1{isat}.Eccentricity;
            omega(isat) = deccfg1{isat}.ArgumentOfPerigee*pi;
            i0(isat) = deccfg1{isat}.Inclination*pi;
            IDOT(isat) = deccfg1{isat}.InclinationRate*pi;
            hterms = num2cell(deccfg1{isat}.HarmonicCorrectionTerms);
            [Cis(isat),Cic(isat),Crs(isat),Crc(isat),Cus(isat),Cuc(isat)] = deal(hterms{:});
            OMEGA_DOT(isat) = deccfg1{isat}.RateOfRightAscension*pi;
            OMEGA0(isat) = deccfg1{isat}.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

Calcule la posición del receptor, lo que requiere dos entradas: pseudodistancias y la posición exacta del satélite en el momento de la transmisión, que se estima en un momento determinado del receptor. Puede calcular estos valores solo si el ejemplo se ejecuta durante al menos 48 segundos de datos. Debido a que la simulación de dichos datos llevaría mucho tiempo, los pseudorangos y las propiedades de posición del satélite (obtenidas al ejecutar el ejemplo durante 50 segundos) se almacenan en los archivos MAT. Si ejecuta la simulación por menos de minStepsForPosEst, entonces estos archivos MAT se cargan y se calcula la posición del receptor. De lo contrario, las propiedades requeridas se toman directamente de la cadena de simulación.

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.0211  -76.8270   21.2084


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

Exploración adicional

Este ejemplo muestra cómo realizar el procesamiento del receptor GPS para la constelación de satélites simulada durante solo 3 segundos. Amplíe el ejemplo a 50 segundos y calcule la posición del receptor para ver cómo funciona el receptor GPS.

Utilice su propio archivo RINEX para configurar el ejemplo y estimar la posición del receptor a partir de la forma de onda generada en este ejemplo.

Apéndice

Este ejemplo utiliza estos datos y archivos auxiliares:

Referencias

[1] IS-GPS-200, Rev: M. Interfaces de Usuario de Navegación/Segmento Espacial GPS NAVSTAR. 21 de mayo de 2021; Código de identificación: 66RP1.

[2] Elliott D. Kaplan y C. Hegarty, eds., Understanding GPS/GNSS: Principios y Aplicaciones, Tercera edición, Serie de Aplicaciones y Tecnología GNSS (Boston; Londres: Casa Artech, 2017).