メインコンテンツ

このページは機械翻訳を使用して翻訳されました。最新版の英語を参照するには、ここをクリックします。

PlutoSDR を使用したGPS信号伝送、取得、トラッキング

この例では、Satellite Communications Toolbox を使用して生成された Global 測位 System (GPS) 波形の無線送受信に ADALM-PLUTO 無線を使用する方法を示します。

はじめに

この例では、ADALM-PLUTO 無線を使用して、複合GPS信号の無線伝送を実行します。複合GPS信号を生成するには、次の手順に従います。

  1. 複数の衛星からレガシーGPS波形を取得します。GPSベースバンド波形を生成するために必要なさまざまなパラメーターを設定する方法の詳細については、GPS 波形生成 の例を参照してください。

  2. 各衛星の波形にドップラーシフトと遅延を追加し、複合信号を形成します。

この ADALM-PLUTO 無線は、 GPSベースバンド波形を繰り返しモードで送信します。同じ ADALM-PLUTO 無線機は、送信されたGPS信号を受信して取得を実行し、取得操作から検出された衛星のコード位相と搬送周波数を追跡できます。この例に示す取得とトラッキングは、coarse取得 (C/A) コードを含むGPS信号用です。ハードウェアをセットアップするには、インストールと設定 ページの指示に従ってください。

ハードウェアを必要としないシミュレーション例として、C/Aコードを使用したGPS受信機の取得とトラッキング の例を参照することもできます。

この図はGPS信号伝送と受信のプロセスの概要を示しています。

GPS_Pluto_SDR.png

この例は、次のセクションで構成されています。

  • シミュレーションパラメーターの設定 — 波形を生成するためのさまざまなパラメーターを設定し、 GPS受信機を設定します。

  • GPS波形の生成 — 複数の衛星からGPS波形を生成し、遅延、ドップラーシフト、ノイズを追加します。

  • PlutoSDR を構成する — Pluto 無線のサンプル レート、中心周波数、および送受信のその他のパラメーターを構成します。

  • 取得とトラッキング — 可視衛星を見つけて、ドップラー周波数オフセットとタイミングオフセットのcoarse値を計算します。着信信号の変化するドップラー オフセットとコード位相を追跡します。

シミュレーションパラメーターの設定

次の手順に従ってシミュレーションパラメーターを構成します。

1. 制御フラグ ShowVisualizationstrue に設定します。これにより、すべての出力プロットを表示できます。波形に含めるGPS衛星の数を指定します。

ShowVisualizations =  true;
% Transmitter Configuration
% Specify the number of GPS satellites in the waveform.
numSat = 4;

2. 受信機に表示される衛星の疑似ランダム ノイズ (PRN) 識別番号 (ID) を指定します。PRNIDs の長さは、指定された衛星の数 numSat 以上である必要があります。PRINDs の長さが numSat の値よりも大きい場合、この例では、シミュレーションに PRINDs の最初の N 要素のみを考慮します。ここで、NnumSat の値です。

PRNIDs = [2; 4; 11; 8]; 

3. データ ビット数 NumNavDataBits を指定します。NumNavDataBits 値を使用して、指定された数のデータ ビットのナビゲーション データのGPS波形を生成します。一般に、 GPSレガシー ナビゲーション データ全体は 25 フレームで構成されます。各フレームは 1500 ビットで構成されます。このナビゲーション データには、エフェメリス、クロック、アルマナックの情報が含まれています。このサイズの波形を生成するには長い時間がかかる可能性があるため、この例では 200 ビットの NumNavDataBits 値のGPS波形を生成します。

% Set this value to control the number of navigation data bits in the
% generated waveform
NumNavDataBits = 200;

4. ナビゲーション データ NavDataBitStartIndex のビット開始インデックスを設定します。NavDataBitStartIndex を使用して、ナビゲーション データのランダム インデックスからGPS波形を生成します。

NavDataBitStartIndex = 1321;

5. GPS波形を送信するサンプル レートを設定します。

SampleRate = 3.41e6;  % In samples/sec

6.信号遅延値を設定します。各衛星はGPS受信機から異なる距離にあります。したがって、各信号に生じる遅延は衛星ごとに異なります。各衛星の C/A コード チップ遅延の数を列ベクトルとして指定します。送信時間は C/A コード チップ期間 0.9775 マイクロ秒の整数倍にはならない可能性が高いため、遅延値は小数になる可能性があります。

sigdelay = [300.34; 587.21; 425.89; 312.88]; % Number of C/A code chips delay

7. ドップラーシフトとドップラーレートを設定します。各衛星の速度と位置は異なるため、各衛星に発生するドップラーシフトとドップラー率も異なります。この例では、ドップラーシフトを正弦波的に変化する周波数オフセットとしてモデル化します。衛星シナリオにおけるレイテンシとドップラーシフトの詳細については、衛星シナリオにおけるレイテンシとドップラー効果の計算 の例を参照してください。

% Initialize peak Doppler shift and Doppler rate for each satellite. This example can track Doppler shift from -10KHz to 10KHz.
peakDoppler = [3589; 4256; 8596; 9568]; % In Hz
dopplerRate = [1000; 500; 700; 500];    % In Hz/sec

受信機の構成

GPS受信機を構成するには、次のパラメーターを設定します。

1. 周波数ロック ループ (FLL)、位相ロック ループ (PLL)、および遅延ロック ループ (DLL) のノイズ帯域幅を設定します。

PLLNoiseBandwidth = 80; % In Hz
FLLNoiseBandwidth = 4;  % In Hz
DLLNoiseBandwidth = 1;  % In Hz

GPS受信機では、トラッキングアルゴリズムは、それぞれ周波数ロック ループ (FLL)、位相ロック ループ (PLL)、遅延ロック ループ (DLL) を使用して、周波数、位相、遅延を追跡します。ループ帯域幅が広いほどトラッキングが速くなりますが、SNR が低い場合にロックが失われる可能性があります。ループ帯域幅が狭いと、SNR が低い場合にパフォーマンスは良好になりますが、位相、周波数、遅延の変化の追跡は遅くなります。これらのプロパティの解釈の詳細については、この例の 取得とトラッキング セクションを参照してください。

2. PLL の積分時間を 1 ミリ秒に設定します。

IntegrationTime = 1e-3; % In seconds

GPS波形を生成する

GPS波形には、同相ブランチに高精度コード (パルス符号) が含まれ、直交位相ブランチに C/A コードが含まれます。一般的に、 GPS衛星はL1 周波数 (1575.42 MHz) で波形を送信します。ただし、この例では、リアルタイムGPS信号との干渉を避けるために、送信周波数を 2.41 GHz に設定しています。この例で生成されたGPSデータの詳細については、IS- GPS-200標準[1]を参照してください。

この例では、各衛星のGPSベースバンド波形を生成し、各衛星信号の遅延を個別にモデル化し、すべての信号を組み合わせて複合信号 . を生成します。次の手順に従ってGPS波形を生成します。

1. 遅延モデリング numBitsForDelay の追加ナビゲーション データ ビットの数を設定します。numBitsForDelay 1 に設定すると、信号に最大 20 ミリ秒が導入されます。より高い遅延値を導入するには、numBitsForDelay を整数値で増やします。

numBitsForDelay = 1;

2. 各衛星の従来のGPS波形を生成するために HelperGPSNavigationConfig オブジェクトを作成します。各 HelperGPSNavigationConfig オブジェクトには、衛星から送信された情報が含まれています。この例では、データのデコードではなく、取得とトラッキングの操作の説明に重点を置いているため、この例では、すべての衛星のパラメーターを一定に保ちます。

%Initialize output waveform
resultsig = 0;

% Generate waveform from each satellite
for isat = 1:numSat
    % Create the legacy navigation (LNAV) configuration object
    lnavConfig = HelperGPSNavigationConfig("SignalType","LNAV","PRNID",PRNIDs(isat));

    % Generate the navigation data bits from the configuration object
    lnavData = HelperGPSNAVDataEncode(lnavConfig);

    % Configure the GPS waveform generation properties
    t = lnavConfig.HOWTOW*6; % First get the initial time
    % HOWTOW is an indication of the next subframe starting point. Because
    % each subframe is 300 bits long, you must subtract 300 bits from the
    % initial value to get the starting value of the first subframe. This
    % value can be negative as well. Because bit is of a 20 millisecond
    % duration, to get the time elapsed for bits, you must multiply the bit
    % index by 20e-3.
    bitDuration = 20e-3; % In sec
    pCodeRate = 10.23e6; % In Hz
    numPChipsPerNavBit = bitDuration*pCodeRate;
    navdatalen = length(lnavData);
    offsetTime = mod(NavDataBitStartIndex-301,navdatalen)*bitDuration;
    inittime = t + offsetTime;

    % To model delay, get one extra navigation bit from the previous bit
    navBitIndices = mod(NavDataBitStartIndex+(-1*numBitsForDelay:(NumNavDataBits-1)),navdatalen);
    navBitIndices(navBitIndices==0) = navdatalen;
    navbits = lnavData(navBitIndices);
    navdata = 1 - 2*navbits;
    upSampledNavData = repelem(navdata,numPChipsPerNavBit,1); 

    % Generate P-code and C/A code
    pgen = gpsPCode("PRNID",PRNIDs(isat),"InitialTime",inittime, ...
        "OutputCodeLength",(NumNavDataBits+numBitsForDelay)*numPChipsPerNavBit);
    pcode = 1 - 2*double(pgen());

    % Reduce the power of the I-branch signal by 3 dB, per IS-GPS-200 [1].
    % See table 3-Va in [1].
    isig = pcode/sqrt(2);

    cacode = 1 - 2*double(gnssCACode(PRNIDs(isat),"GPS"));

    numCACodeBlocks = (NumNavDataBits + numBitsForDelay)*bitDuration*1e3;
    caCodeBlocks = repmat(cacode(:),int64(numCACodeBlocks),1);

    % Because C/A code is 10 times slower than P-code, repeat each sample
    % of C/A code 10 times
    qsig = repelem(caCodeBlocks,10,1);

    % Generate the baseband waveform
    gpsBBWaveform = (isig + 1j*qsig).*upSampledNavData;

    % Initialize the number of samples per bit
    numSamplesPerBit = SampleRate*bitDuration;
    
    % Introduce Doppler
    numSamplesGPSBB = length(gpsBBWaveform);
    sampleIndices = (0:(numSamplesGPSBB-1));
    ph = sin(dopplerRate(isat)*sampleIndices/(peakDoppler(isat)*10.23e6));
    phase = 2*pi*(peakDoppler(isat)^2)/dopplerRate(isat)*ph;
    bbwave = gpsBBWaveform(:).*exp(1j*phase(:));

    % Rate match the generated signal to the radio sample rate
    [upfac,downfac] = rat(SampleRate/10.23e6);
    upgcode = repelem(bbwave,upfac,1);
    gpsWaveform = upgcode(1:downfac:end);

    % Get the number of samples for delay
    caCodeRate = 1.023e6;
    numDelaySamples = floor(sigdelay(isat)*SampleRate/caCodeRate);

    % Add delay to the signal by keeping samples of the previous bit at the
    % beginning of the signal
    delayedSig = gpsWaveform(numSamplesPerBit-numDelaySamples+1:end); 

    % Remove the final samples to make all signals of equal length
    delayedSig = delayedSig(1:end-numDelaySamples);

    % Get the composite signal by adding the current satellite signal
    resultsig = resultsig + delayedSig;
end

PlutoSDR を設定する

サンプル レート、送信周波数、送信機ゲインを設定して、Pluto 無線送信機を構成します。

% Configure Pluto radio transmitter
fs = SampleRate;
fc = 2.41e9; % center frequency for transmission and reception
tx = sdrtx('Pluto');
tx.CenterFrequency = fc;
tx.BasebandSampleRate = fs;
tx.Gain = -33;
transmitRepeat(tx,resultsig);
## Establishing connection to hardware. This process can take several seconds.
## Waveform transmission has started successfully and will repeat indefinitely. 
## Call the release method to stop the transmission.

受信中心周波数と各受信フレームのサンプル数を設定して、Pluto 無線受信機を構成します。同じ無線で送信と受信の両方を実行するため、受信機のサンプル レートを送信機のサンプル レートと同じに設定します。

% Configure pluto radio receiver
rx = sdrrx("Pluto");
rx.CenterFrequency = fc;
rx.BasebandSampleRate = fs;
rx.SamplesPerFrame = 102300;
rx.OutputDataType =  "single";
recordDuration = 0.7; % time duration for receiving data, in seconds
rxwaveform = [];
ovrflw_Cnt = 0; % count number of overflows to check discontinuities in reception
loopCnt = round(recordDuration/(rx.SamplesPerFrame/fs));

for i = 1:loopCnt
    [y1,~,of] = rx();
    ovrflw_Cnt = of+ovrflw_Cnt;
    rxwaveform = [rxwaveform; y1];
end
## Establishing connection to hardware. This process can take several seconds.
release(tx);
release(rx);

取得とトラッキング

取得モジュールは、すべての可視衛星を検出し、各衛星からの受信信号のcoarse ドップラー シフトとcoarse時間オフセットを推定します。この例では、受信機がコールド スタート モードで開始し、32 個のGPS衛星すべてを探して検索を開始すると想定しています。

トラッキングモジュールは、取得後に残るドップラーシフトとコード位相オフセットを補正します。搬送周波数、搬送位相、コード遅延の 3 つのパラメーターを追跡する必要があります。各パラメーターを追跡するには、フィードバック ループを使用します。FLL を使用して搬送周波数を追跡し、PLL を使用して位相を追跡し、DLL を使用して遅延 (コード位相オフセット) を追跡します。

GPS受信機の起動モード、FLL、PLL、DLL を使用した取得アルゴリズムの詳細については、C/Aコードを使用したGPS受信機の取得とトラッキング の例を参照してください。

取得を実行するには、gnssSignalAcquirer オブジェクトを初期化し、そのプロパティを構成します。

initialsync = gnssSignalAcquirer;
initialsync.SampleRate = SampleRate
initialsync = 
  gnssSignalAcquirer with properties:

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

トラッキングに必要なパラメーターを初期化します。可視衛星を見つけて追跡するために取得を実行します。

% Consider data that is 1 millisecond long.
numSamples = ceil(SampleRate*IntegrationTime);
[allRxInput,prevSamples] = buffer(rxwaveform,numSamples);
nFrames = size(allRxInput,2);
numdetectsat = 0;
PRNIDsToSearch = 1:32;

for iBuffer = 1:nFrames
    bufferWave = allRxInput(:,iBuffer);

    if iBuffer == 1
        % This example assumes a hot start for all the satellites. Hence,
        % acquisition performed only once in this example. When decoding
        % the almanac data, based on the available satellites, you can
        % perform acquisition for the visible satellites only.
        numSamplesForInitSync = SampleRate*1e-3; % 1 milliseccond
        [y,corrval] = initialsync(bufferWave(1:numSamplesForInitSync),1:32);
        PRNIDsToSearch = (y(y(:,4).IsDetected,1).PRNID).';
        doppleroffsets = (y(y(:,4).IsDetected,2).FrequencyOffset).';
        codephoffsets = (y(y(:,4).IsDetected,3).CodePhaseOffset).';

        % In general, almanac files offer information about available
        % satellites. Because this example does not perform data decoding,
        % it depends on the acquisition results for satellite detection.
        numdetectsat = length(PRNIDsToSearch);

        disp(['The detected satellite PRN IDs: ' num2str(PRNIDsToSearch)])
        if(numdetectsat~=0)
            if ShowVisualizations == 1
                figure;
                mesh(initialsync.FrequencyRange(1):initialsync.FrequencyResolution:initialsync.FrequencyRange(2), ...
                    0:size(corrval,1)-1,corrval(:,:,1));
                xlabel("Doppler Offset")
                ylabel("Code Phase Offset")
                zlabel("Correlation")
                msg = "Correlation Plot for PRN ID: " + PRNIDsToSearch(1);
                title(msg)
            end
        end

        % Initialize all the properties which must be accumulated.
        accuph = zeros(nFrames,numdetectsat); % Each column represents data from a satellite
        accufqy = zeros(nFrames,numdetectsat);
        accufqyerr = zeros(nFrames,numdetectsat);
        accupherr = zeros(nFrames,numdetectsat);
        accuintegwave = zeros(nFrames,numdetectsat);
        accudelay = zeros(nFrames,numdetectsat);
        accudelayerr = zeros(nFrames,numdetectsat);

        % Create the signal tracker object that tracks phase, frequency,
        % and delay in the signal
        carrierCodeTrack = gnssSignalTracker;
        carrierCodeTrack.SampleRate = SampleRate;
        carrierCodeTrack.IntermediateFrequency = 0;
        carrierCodeTrack.PLLNoiseBandwidth = PLLNoiseBandwidth;
        carrierCodeTrack.FLLNoiseBandwidth = FLLNoiseBandwidth;
        carrierCodeTrack.DLLNoiseBandwidth = DLLNoiseBandwidth;
        carrierCodeTrack.IntegrationTime = IntegrationTime;
        carrierCodeTrack.PRNID = PRNIDsToSearch;
        carrierCodeTrack.InitialFrequencyOffset = doppleroffsets;
        carrierCodeTrack.InitialCodePhaseOffset = codephoffsets;
    end


    [integwave,trackInfo] = carrierCodeTrack(bufferWave);

    % Accumulate the values to see the results at the end
    accuintegwave(iBuffer,:) = integwave;
    accufqyerr(iBuffer,:) = trackInfo.FrequencyError;
    accufqy(iBuffer,:) = trackInfo.FrequencyEstimate;
    accupherr(iBuffer,:) = trackInfo.PhaseError;
    accuph(iBuffer,:) = trackInfo.PhaseEstimate;
    accudelayerr(iBuffer,:) = trackInfo.DelayError;
    accudelay(iBuffer,:) = trackInfo.DelayEstimate;
end
The detected satellite PRN IDs: 2   8  11   4

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

trackedSignal = accuintegwave; % Useful for further GPS receiver processing

トラッキングループの結果をプロットします。推定位相誤差は、周波数ロック ループが収束するまで収束しません。周波数の誤差によって位相の誤差が発生するため、このような動作が予想されます。また、ドップラーシフトの変化により、PLL と FLL の出力は時間とともに変化するように見えます。

if ShowVisualizations == 1
    for isat = 1 % To see tracking results for all the detected satellites by using above line
        groupTitle = ['Tracking Loop Results for Satellite PRN ID:', ...
            num2str(PRNIDsToSearch(isat))];

        figure

        % Plot the frequency discriminator output
        subplot(2,1,1)
        plot(accufqyerr(:,isat))
        xlabel('Milliseconds')
        ylabel('Frequency Error')
        title('Frequency Discriminator Output')

        % Plot the FLL output
        subplot(2,1,2)
        plot(accufqy(:,isat))
        xlabel('Milliseconds')
        ylabel('Estimated Frequency Offset')
        title('FLL Output')
        sgtitle(['FLL ' groupTitle])

        figure

        % Plot the phase discriminator output
        subplot(2,1,1)
        plot(accupherr(:,isat))
        xlabel('Milliseconds')
        ylabel('Phase Error')
        title('Phase Discriminator Output')

        % Plot the PLL output
        subplot(2,1,2)
        plot(accuph(:,isat))
        xlabel('Milliseconds')
        ylabel('Estimated Phase')
        title('PLL Output')
        sgtitle(['PLL ' groupTitle])

        figure

        % Plot the delay discriminator output
        subplot(2,1,1)
        plot(accudelayerr(:,isat))
        xlabel('Milliseconds')
        ylabel('Delay Error')
        title('Delay Discriminator Output')
        
        % Plot the DLL output
        subplot(2,1,2)
        plot(accudelay(:,isat))
        xlabel('Milliseconds')
        ylabel('Estimated Delay')
        title('DLL output in the units of number of C/A-code chips')
        sgtitle(['DLL ' groupTitle]) 
    end
end

Figure contains 2 axes objects and another object of type subplottext. Axes object 1 with title Frequency Discriminator Output, xlabel Milliseconds, ylabel Frequency Error contains an object of type line. Axes object 2 with title FLL Output, xlabel Milliseconds, ylabel Estimated Frequency Offset contains an object of type line.

Figure contains 2 axes objects and another object of type subplottext. Axes object 1 with title Phase Discriminator Output, xlabel Milliseconds, ylabel Phase Error contains an object of type line. Axes object 2 with title PLL Output, xlabel Milliseconds, ylabel Estimated Phase contains an object of type line.

Figure contains 2 axes objects and another object of type subplottext. Axes object 1 with title Delay Discriminator Output, xlabel Milliseconds, ylabel Delay Error contains an object of type line. Axes object 2 with title DLL output in the units of number of C/A-code chips, xlabel Milliseconds, ylabel Estimated Delay contains an object of type line.

トラッキング後、信号のコンスタレーションをプロットします。

if ShowVisualizations == 1
    rxconstellation = comm.ConstellationDiagram(1,"ShowReferenceConstellation",false);
    demodsamples = accuintegwave(301:end,:)/rms(accuintegwave(:));
    if ~isempty(demodsamples)
        rxconstellation(demodsamples(:));
    end
end 

コンスタレーション図の中央にコンスタレーション点が表示される場合、それは捕捉アルゴリズムによる衛星の誤検出に相当します。

さらなる探究

この例では、4 つのGPS衛星のGPS信号の取得とトラッキングを実行し、ADALM-Pluto 無線を使用して無線伝送を行う方法を説明しました。さらに詳しく調べるには、次のバリエーションを試してください。

  • 可視GPS衛星の数を増やして、結果を観察します。ドップラー、SNR、遅延を適切に初期化します。

  • GPS信号の振幅をノイズフロア以下にスケーリングしてリアルタイムGPS信号をシミュレートし、ADALM-PLUTO 無線を使用して取得とトラッキングを確認します。

  • SDR デバイスを使用して位置推定を実行します。2 つの PlutoSDR デバイス (送信用と受信用に 1 つずつ) を使用する必要がある場合があります。位置推定の詳細については、End-to-End GPS Legacy Navigation Receiver Using C/A-Code の例を参照してください。

付録

この例では、次のデータとヘルパー ファイルを使用します。

参考文献

[1] IS- GPS-200、改訂:L. NAVSTAR GPS宇宙セグメント/ナビゲーション ユーザー セグメント インターフェース。2020年5月14日; コード識別:66RP1。