使用NI USRP 生成 GNSS 信号

东枫科技

于 2025-08-09 02:04:45 发布

阅读量467

点赞数 9

CC 4.0 BY-SA版权

文章标签： USRP GNSS GPS

本文链接：https://blog.csdn.net/qq_36666115/article/details/150087880

本示例演示如何生成全球定位系统 (GPS) 波形并使用软件定义无线电 (SDR) 进行传输。您可以使用此示例生成的 GPS 波形来测试您的接收器。虽然本示例描述的是 GPS 波形的生成，但您可以对其进行扩展，以适应其他全球导航卫星系统 (GNSS) 技术。

有关如何为单个卫星生成 GPS 波形的更多详细信息，请参阅GPS 波形生成示例。

此示例遵循以下步骤。

使用 GPS 年历文件来模拟卫星星座。
指定接收器位置以创建真实的 GPS 波形场景。
模拟卫星环境来计算卫星相对于建模接收器的相对位置。
计算此场景中接收信号的多普勒频移、延迟和功率。
gpsWaveformGenerator使用System object™生成基带 GPS 波形。
使用步骤 4 的结果将多普勒频移和延迟等损伤引入到生成的基带信号中。
通过SDR传输受损的基带波形。

该图说明了波形的产生过程。

在这里插入图片描述
使用此示例，您可以将 GPS 波形保存到文件中，或者使用支持的 SDR 通过无线方式传输信号。

在这里插入图片描述

初始化和配置参数

GPS波形生成示例展示了仅从一颗卫星生成波形，且没有任何损伤。在此示例中，您可以通过模拟卫星场景来生成多颗卫星波形，并使用 SDR 进行无线传输。

指定配置参数来生成波形。

useSDR            = false;                  % In the default case, do not use SDR
WriteWaveToFile   = false;
waveformFileName  = "gpsBBWaveform.bb"; % File extension is always .bb
signalType        = "GPS C/A"; % Possible values: "GPS C/A" | "GPS L1C", "GPS L2C" | "GPS L5"
centerFrequency   = 1575.42e6; % Possible values: L1 (1575.42 MHz), L2 (1254.35 MHz), and L5 (1176.45 MHz)
sampleRate        = 5000000; % In Hz
rxlla             = [17.4349,78.3827,19]; % Receiver position in latitude (degrees), longitude (degrees), and altitude (m)
waveDuration      = 10; % Waveform duration in seconds
enableImpairments =  true;                 % Option to enable or disable impairments
minElevationAngle = 10; % In degrees. Below this elevation angle satellites are not considered in simulation
seed              = 73; % Seed for getting reproducible results
useSEMAlmanac     = "stored"; % File to derive ephemeris data and satellite simulation
if useSEMAlmanac == "stored"
    % Provide the system effectiveness model (SEM) almanac file name
    almFileName   = "gpsAlmanac.txt"; % Almanac file name
    % Provide the startTime
    startTime = datetime(2021,6,24,0,0,48,TimeZone="UTC");
else
    % Download the latest SEM almanac file from the Navigation Center
    % website and store the file
    url = "https://www.navcen.uscg.gov/sites/default/files/gps/almanac/current_sem.al3";
    currentDay = string(datetime("today",TimeZone="UTC"));
    almFileName = "gps_SEM_" + currentDay + "_UTC.txt";
    websave(almFileName,url)

    % Uses the latest SEM almanac file. Set the startTime to current date
    % and time.
    startTime = datetime("now",TimeZone="UTC");
end

根据指定的signalType，初始化一些配置参数。

wavegenobj = gpsWaveformGenerator(SampleRate=sampleRate);
switch(signalType)
    case "GPS C/A"
        wavegenobj.SignalType = "legacy";
        navDataType = "LNAV";
    case "GPS L1C"
        wavegenobj.SignalType = "l1C";
        navDataType = "CNAV2";
    case "GPS L2C"
        wavegenobj.SignalType = "l2C";
        navDataType = "CNAV";
    case "GPS L5"
        wavegenobj.SignalType = "l5";
        navDataType = "L5";
end
stepTime = wavegenobj.BitDuration;          % Generate waveform in the step size corresponding to bit duration

计算给定卫星场景的多普勒频移和延迟。

% Initialize satellite scenario
sc = satelliteScenario;
% Set up the satellites based on the RINEX data
sat = satellite(sc,almFileName,OrbitPropagator="gps");
rx = groundStation(sc,rxlla(1),rxlla(2),Altitude=rxlla(3)); % Set up the receiver which doesn't move
rx.MinElevationAngle = minElevationAngle;

sc.StartTime = startTime;
sc.StopTime = sc.StartTime + seconds(waveDuration-stepTime);
sc.SampleTime = stepTime;

% Calculate Doppler shift and latency over time for all the visible satellites
dopShifts = dopplershift(sat,rx,Frequency=centerFrequency).';
ltncy = latency(sat,rx).';

初始化有助于信噪比 (SNR) 计算的值。

c = physconst("LightSpeed"); % Speed of light in m/sec
Pt = 44.8;                   % Typical transmission power of GPS 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

根据自由空间路径损耗方程计算接收器的功率。

Pr   = Pt*DtLin*DrLin./ ...
    ((4*pi*(centerFrequency+dopShifts).*ltncy).^2);

计算每颗可见卫星的信噪比。添加一些信号功率（此处为 3 dB），以弥补使用 SDR 传输时的损耗。

snrs = 10*log10(Pr/(k*T*sampleRate)) + 3;

生成波形并添加损伤

从年历初始化配置文件。

satIndices = find(~isnan(ltncy(1,:)));
navcfg = HelperGPSAlmanac2Config(almFileName,navDataType,satIndices,startTime);
visiblesatPRN = [navcfg(:).PRNID]

visiblesatPRN = 1×8

    10    13    15    20    21    24    29    32

生成导航数据。

% Generate GPS navigation data
tempnavdata = HelperGPSNAVDataEncode(navcfg(1));
navdata = zeros(length(tempnavdata),length(navcfg));
navdata(:,1) = tempnavdata;
for isat = 2:length(navcfg)
    navdata(:,isat) = HelperGPSNAVDataEncode(navcfg(isat));
end

配置GPS波形生成对象。

wavegenobj.PRNID = visiblesatPRN

wavegenobj = 
  gpsWaveformGenerator with properties:

           SignalType: "legacy"
                PRNID: [10 13 15 20 21 24 29 32]
          EnablePCode: false
    HasDataWithCACode: true
           SampleRate: 5000000

  Show all properties

gnsschannelobj = HelperGNSSChannel(FrequencyOffset=dopShifts(1,satIndices), ...
    SignalDelay=ltncy(1,satIndices), ...
    SignalToNoiseRatio=snrs(1,satIndices), ...
    SampleRate=sampleRate,RandomStream="mt19937ar with seed",Seed=seed)

gnsschannelobj = 
  HelperGNSSChannel with properties:

               SampleRate: 5000000
    IntermediateFrequency: 0
       DisableImpairments: false
          FrequencyOffset: [956.0901 -3.2174e+03 -2.5058e+03 316.0502 1.7671e+03 1.7725e+03 -1.8776e+03 3.4268e+03]
              SignalDelay: [0.0767 0.0808 0.0715 0.0728 0.0744 0.0771 0.0691 0.0803]
       SignalToNoiseRatio: [-11.2812 -11.7307 -10.6680 -10.8213 -11.0074 -11.3176 -10.3708 -11.6799]
             RandomStream: "mt19937ar with seed"
                     Seed: 73

numsteps = round(waveDuration/stepTime);
samplesPerStep = sampleRate*stepTime;

初始化波形。如果启用损伤，波形中的通道数将设置为 1，因为所有卫星信号都会组合形成单个数据流。如果禁用损伤，则通道数等于可见卫星的数量。

numchannels = 1*enableImpairments + length(visiblesatPRN)*(~enableImpairments);
gpswaveform = zeros(numsteps*samplesPerStep,numchannels);

可选地，初始化将波形写入文件的对象。

if WriteWaveToFile == true
    bbwriter = comm.BasebandFileWriter(waveformFileName,sampleRate,0);
end

生成波形。

for istep = 1:numsteps
    idx = (istep-1)*samplesPerStep + (1:samplesPerStep);
    navbit = navdata(istep,:);
    tempWaveform = wavegenobj(navbit);
    if enableImpairments == true
        gpswaveform(idx,:) = gnsschannelobj(tempWaveform);
    else
        gpswaveform(idx,:) = tempWaveform;
    end
    if WriteWaveToFile == true
        bbwriter(gpswaveform(idx,:))
    end
    
    % Update the properties of the channel object based on the values
    % calculated from the satellite scenario
    gnsschannelobj.SignalToNoiseRatio = snrs(istep,satIndices);
    gnsschannelobj.FrequencyOffset = dopShifts(istep,satIndices);
    gnsschannelobj.SignalDelay = ltncy(istep,satIndices);
end

可视化波形的频谱。

scope = spectrumAnalyzer(SampleRate=sampleRate,AveragingMethod="exponential",ForgettingFactor=1);
scope.SpectrumUnits = "dBW";
% Plot only 1 sec of the waveform spectrum.
wlen = sampleRate; % Number of samples corresponding to one second
if length(gpswaveform) > wlen
    scope(gpswaveform(end-wlen+1:end))
else
    % If the length of the generated waveform is less than one second,
    % visualize the spectrum of the entire waveform
    scope(gpswaveform)
end

在这里插入图片描述

设置电台

调用(无线测试台)函数。该函数返回您使用无线电设置(无线测试台)向导保存的所有可用无线电设置配置。radioConfigurations

if useSDR == 1
    savedRadioConfigurations = radioConfigurations;

要使用已保存的无线电设置配置名称更新下拉菜单，请点击“更新”。然后选择本示例要使用的无线电。

    savedRadioConfigurationNames = [string({savedRadioConfigurations.Name})];
    if ~isempty(savedRadioConfigurationNames)
        radio = savedRadioConfigurationNames(1) ;
    end
end

通过 SDR 传输波形

一旦连接 SDR，就通过无线电广播 GPS 信号。

if useSDR == 1 & ~isempty(savedRadioConfigurationNames) % Using SDR
    if ~exist("txradio","var")
        txradio = basebandTransmitter(radio);
        txradio.RadioGain = 10;
        txradio.CenterFrequency = centerFrequency;
        txradio.SampleRate = sampleRate
    end
    transmit(txradio,gpswaveform,"once")
end