GPS实验

实验文档

一． 实验模块具体内容

实验1 将载波从输入信号中“剥离”出来

一．函数模块解释

函数功能：将载波从输入信号中剥离出来，将输入信号与本地载波信号分别相乘来实现。 函数名称：[I1, Q1,I2, Q2]=peelcarrier(signal2 ,signal2,sinCarr,cosCarr) 函数输入：

signal1 ——第1ms输入信号；

signal2 ——第2ms输入信号；

? sinCarr ——本地载波的正弦分量； cosCarr ——本地载波的余弦分量；

函数输出：I1 ——第1ms输入信号对应的同相分量； Q1——第1ms输入信号对应的正交分量； ? I2 ——第2ms输入信号对应的同相分量； Q2 ——第2ms输入信号对应的正交分量； 二．函数流程图

输入变量Signal1 ,signal2,sinCarr, cosCarr 将输入信号与本地载波分量分别相乘，将载波从输入信号中剥离出来

输出变量I1，Q1 I2，Q2

三．实验操作

1.实现代码编写：

function [I1, Q1,I2, Q2]=peelcarrier(signal2 ,signal2,sinCarr,cosCarr) I1 = sinCarr .* signal1; Q1 = I2 = Q2 = End

2.实验数据：

（1）signal.mat 里面是第1ms输入信号和第2ms输入信号 （2）localsig.mat 里面是本地载波的正弦分量和余弦分量

（3）Removecar_result.mat 是剥离载波的结果，用来与实验结果对照 3.实验要求：

（1） 完善上述的函数代码；

（2） 设计实验，利用所给的实验数据验证代码的正确性，获得基带信号。

实验2 用粗捕获的结果求相关峰比值并判断卫星号PRN是否存在

一．函数模块解释

函数作用：在所对应PRN号卫星的相关计算结果矩阵中找出相关峰对应的多普勒频移和码相位，并检测卫星的存在性。

函数名称：function [acqResults,peakSize,secondPeakSize]=Findprn(PRN,results,settings) 输入：PRN --PRN 数字：01 至 32 表天空使用中的卫星编号 results--相关计算结果矩阵（对应相应的PRN号的卫星） settings--全局设置值 输出：

peakSize --相关峰

codePhase--相关峰的码相位

secondPeakSize--第二高相关峰的峰值

acqResults.peakMetric(PRN)--对应卫星号相关峰比值 acqResults.codePhase(PRN)--对应卫星号的码相位 二．函数流程图

得到的捕获结果results是一矩阵：行数表示搜索的频率范围；列数表示搜索的采样点数通过二维搜索找到结果中峰值对应的频率序号frequencyBinIndex和采样点序号codePhase根据码片采样点数samplesPerCodeChip和峰值对应的采样点序号codePhase求出码相位范围codePhaseRange在frequencyBinIndex的对应行，在codePhaseRange范围内寻找第二个峰值将第一个峰值比上第二个峰值，得比值peakMetric（对应PRN号卫星）打印信息，搜索下一颗卫星NpeakMetric>acqThreshold?Y进行后续处理判断PRN号卫星信号是否存在

三．实验操作

1.实现函数代码的编写补全：

function [acqResults,peakSize,secondPeakSize]=Findprn(PRN,results,settings) %初始化acqResults

acqResults.codePhase = zeros(1, 32); acqResults.peakMetric = zeros(1, 32); %代码补充

[peakSize frequencyBinIndex] = max(max(results, [], 2)); [peakSize codePhase] =

acqResults.codePhase(PRN) =

samplesPerCodeChip = round(settings.samplingFreq / settings.codeFreqBasis); excludeRangeIndex1 = codePhase - samplesPerCodeChip; excludeRangeIndex2 = codePhase + samplesPerCodeChip; %对于每个扩频码求采样数

samplesPerCode = round(settings.samplingFreq / ...

(settings.codeFreqBasis / settings.codeLength)); % 获得第二相关峰的搜索范围 if excludeRangeIndex1 < 2

codePhaseRange = excludeRangeIndex2 : ...

(samplesPerCode + excludeRangeIndex1); elseif excludeRangeIndex2 >= samplesPerCode

codePhaseRange = (excludeRangeIndex2 - samplesPerCode) : ... excludeRangeIndex1; else

codePhaseRange = [1:excludeRangeIndex1, ...

excludeRangeIndex2 : samplesPerCode]; end

secondPeakSize = max(results(frequencyBinIndex, codePhaseRange));

% 代码补充

acqResults.peakMetric(PRN) = % 判断卫星

if (peakSize/secondPeakSize) > settings.acqThreshold %save('F:\\test code\\实验2.mat ', PRN); fprintf('\\n\\n卫星存在：d \\n\\n', PRN); else

fprintf('\\n\\n--该卫星不存在-- \\n\\n'); end end

2.实验数据：

（1）文件中的results1.mat；results2.mat；、、、、、、；results32.mat共32个表格，是对应的搜索的卫星PRN号（1到32）的信号相关计算矩阵结果，调用结果时，用的是results2.mat文件，那么对应的PRN号就是2；

（2）文件initSettings.m是用来获得初始化数据的，用settings = initSettings();这个语句就可以获得初始化的设置数值；

（3）结果说明，这个函数是用相关结果来获得峰值比的，然后判断该卫星是否存在，最后，保留卫星信号的载波信息和卫星的相关峰比值等，32 组数据中，只有PRN号为29的卫星是存在的。 3.实验要求：

（1）完善上述的函数代码；

（2）利用所给的实验数据验证自己编写代码的正确性。

实验3 精细搜索找出精细载波频率

一．函数模块解释

函数功能：在粗捕获到卫星prn号前提下，根据连续载波cw(xcarrier)求FFT点数，并将连续载波cw变换到频域；去除冗余频率分量，求最大频率分量及其位置；生成频率向量（数组），找出输入信号中的载波频率。

函数名称：function [acqResults]=findcarrFreq(PRN,xCarrier,samplingFreq) 函数输入：

PRN ——捕获到的卫星号； xCarrier——连续载波信号； samplingFreq—— 采样频率； 函数输出：

acqResults——结构体，这里用来记录卫星PRN捕获到的精细载波频率acqResults.carrFreq(PRN)； 二．函数实现流程图

输入连续载波和采样频率，及捕捉的卫星PRN1.用8*(2^(nextpow2(信号长度)))获得对该信号进行fft变换的点数2.将连续载波cw变换到频域并取模值3.在fftNumPts各频率分量中，只取一半就可以，uniqFftPts = ceil((fftNumPts + 1) / 2)4.在fftxc中，对 从第5个频率分量开始长度为uniqFftPts的一段频率分量计算结果进行搜索，找出其中的最大值，返回对应下标号fftMaxIndex。5.生成频率向量（数组）fftFreqBins = (0 : uniqFftPts-1) * samplingFreq/fftNumPts;找出对应标号fftMaxIndex载波频率 三．实验操作 1.实验代码补充：

function [acqResults]=findcarrFreq(PRN,xCarrier,samplingFreq) %

% %求已经解扩的零直流信号的FFT点数 fftNumPts =

acqResults.carrFreq=zeros(1,32);

%对xCarrier进行fftNumPts点的FFT运算，并求模值 fftxc =

%功能：去除冗余频率分量，求最大频率分量及其位置 uniqFftPts =

[fftMax, fftMaxIndex] =

fftFreqBins = (0 : uniqFftPts-1)*samplingFreq/fftNumPts; acqResults.carrFreq(PRN) =

输入相应的跟踪变量和每颗卫星的c/a码起点及初值1.初始化接收机接收卫星信号的时间为无穷大2.根据通道c/a码记录结果，求每个卫星的传输时间3.然后求所有传输时间相对最先到达的卫星的时间差3.用相对时延差加上设置的一般到达时间作为该卫星到达接收机的时间pseudoranges：用到达时间乘以光速获得卫星到接收机的观测伪距，作为输出（注意统一单位，ms转化为s） 三．实验操作 1.实验代码补充：

function [pseudoranges] = calculatePseudoranges(absoluteSample,settings) % absoluteSample是导航C/A码起始索引点找到对应的数字点，除以每码采样点得到传输时间

% absoluteSample是一个四维的数组，包含四颗星的数据 %初始化各卫星传输时间为无穷大

travelTime = inf(1,length(absoluteSample));

% 求每码采样数据点（samplesPerCode =16368） samplesPerCode = round(settings.samplingFreq / ...

(settings.codeFreqBasis / settings.codeLength)); %补充代码

%求每一个卫星的数据点对应的时间 travelTime=

minimum = floor(min(travelTime)); %求相对时延 dt=

% 求估计的传输时间=时延+设定的初值（到达时间） travelTime = %求伪距

pseudoranges = end

2.实验数据说明：

（1）用initSettings()函数获得settings的设置初值，用于伪距的计算； （2）输入的绝对采样点取值为如下数组：

absoluteSample=[0,35935,47222,-15232]

（3）实验结果见说明文档。 3.实验要求：

（1）补充完整实验代码，并设计程序，根据所给输入数据和结果验证代码的正确性。

实验11 求卫星位置和卫星钟差

一．函数模块解释

函数功能：利用测量的伪距求卫星在ECEF坐标系中的坐标位置和卫星钟差，其中的求解原理见GPS简化版文件。

函数名称：function [satPositions, satClkCorr] = satpos(transmitTime, prnList, ... eph, settings) 函数输入：

transmitTime ——卫星信号到接收机的粗略传输时间； prnList——要处理的伪随机码列； eph—— 卫星星历表；

settings—— 接收机设置的值（包含多项设置）；

函数输出：

satPositions——卫星位置（在ECEF系统中，坐标是[x; y; z]）； satClkCorr—— 卫星时钟校正

二．函数实现流程图

初始化常数：取定地球的公转速度，自转速度，引力常数GM，以及其他常数的设置；同时初始化结果为0；求每颗卫星在无摄情况下的参量：求发射时刻GPS系统时间修正项：1.观测时间的归一化用最终偏近点1.Dt在计算的过程中与2.计算卫星的平均角速度n角E求时间相偏近点角E有关，所以3.计算在t时刻的平近点角M需要最终求得的E来进对修正4.由平近点角计算片偏近点角E（迭代求解）行计算5.计算真近点角V2.利用相关公式求卫星6.计算升交点角矩phi时钟修正项求卫星受摄动时的轨道参数变化：1.由星历数据中获得6个参考变量值2.然后计算升交角距、轨道半径和轨道倾角的摄动修正项3.计算修正后的升交角距，轨道半径，轨道倾角求卫星ECEF坐标：将卫星由参考时间点的轨道平面旋转到ECEF坐标系内：通过卫星的椭圆轨道参数计算卫星的ECEF坐标值结束 三．实验操作

1.结合公式和流程图补充本块实验代码：

%################################################################## % 程序小块：4.2-2

% 作用： 根据校正值计算卫星受到摄动时各个修正项及修正后的参数 % 对应文档中：B计算卫星受摄时的轨道修正量

%################################################################## u = phi + ...

eph(prn).C_uc * cos(2*phi) + ... eph(prn).C_us * sin(2*phi); % 校正后的卫星地心相径

r = a * (1 - eph(prn).e*cos(E)) + ... eph(prn).C_rc * cos(2*phi) + ... eph(prn).C_rs * sin(2*phi); % 校正后的轨道倾角

i = eph(prn).i_0 + eph(prn).iDot * tk + ... eph(prn).C_ic * cos(2*phi) + ... eph(prn).C_is * sin(2*phi); % 计算升交点的经度

Omega = eph(prn).omega_0 + (eph(prn).omegaDot - Omegae_dot)*tk - ... Omegae_dot * eph(prn).t_oe;

Omega = rem(Omega + 2*gpsPi, 2*gpsPi); %求卫星ECEF坐标

satPositions(1, satNr) = cos(u)*r * cos(Omega) - sin(u)*r * cos(i)*sin(Omega); satPositions(2, satNr) = cos(u)*r * sin(Omega) + sin(u)*r * cos(i)*cos(Omega); satPositions(3, satNr) = sin(u)*r * sin(i); %本块作用：求该卫星的时钟修正项

satClkCorr(satNr) = (eph(prn).a_f2 * dt + eph(prn).a_f1) * dt + ... eph(prn).a_f0 - ... eph(prn).T_GD + dtr; 2.实验要求：

（1）根据原理将上述代码的原理在书上找到相应的公式，并加以注释和理解；

（2）根据satpos.m把实验代码与流程图联系起来理解，并结合书上原理，理解求卫星位置的函数过程。

实验12最小二乘法求接收机位置

一．函数模块解释

函数功能： 根据卫星的估计位置和钟差，以及测量的伪距来求接收机位置和时钟误差，用的算法是最小二乘法，原理见书。

函数名称：function [pos, el, az, dop] = leastSquarePos(satpos, obs, settings) 函数输入：

satpos ——ECEF坐标系中卫星坐标； obs——测量的伪距结果；

settings—— 接收机设置的值（包含多项设置）； 函数输出：

pos——接收机位置和时钟误差（在ECEF系统中，坐标是[x, y, z, dt]）； el—— 卫星仰角角度； az——方位角；

dop—— 误差放大因子； 二．函数实现流程图

已知数据是卫星位置（通过函数调用得出），并设置迭代次数为7次初始化：将初始的用户位置初值和时钟偏差初值以及伪距初值都设置为0，（为由下面的迭代法求用户位置做准备）判断是否是首次迭代？是得到所有已知的卫星位置，并将结果放置于矩阵Rot_X中，然后并令trop=2（trop代表对流层误差修正值）否（1）求该次迭代得到的卫星的预测伪距（2）求第I颗卫星的传输时间（3）根据传输时间修正卫星坐标系（4）求卫星的仰角和方位角判断是否需要进行对流层误差修正？是计算对流层延迟误差torp值否设定对流层误差trop=0根据上述修正后的卫星位置和伪距误差修正求伪距变化量（程序中omc(i)值）否根据定位方程式，确定系数矩阵A求矩阵A的秩，判定方程组是否有解否用左除法求得用户位置变化量，然后叠加在用户位置的初值上，产生新的用户位置是判断是否已经达到已定的迭代次数是定位结果：得到用户的坐标位置和用户时钟偏差,由最终的矩阵A计算误差放大因子DOP

三．实验操作 1.实验代码补充：

function [pos, el, az, dop] = leastSquarePos(satpos, obs, settings) %用最小二乘的迭代算法求接收机的位置，原理参照课本

%迭代次数设置 nmbOfIterations = 7; dtr = pi/180; pos = zeros(4, 1); X = satpos;

nmbOfSatellites = size(satpos, 2); A = zeros(nmbOfSatellites, 4); omc = zeros(nmbOfSatellites, 1); az = zeros(1, nmbOfSatellites); el = az; % 开始迭代运算

for iter = 1:nmbOfIterations

for i = 1:nmbOfSatellites if iter == 1

%--- Initialize variables at the first iteration -------------- Rot_X = X(:, i); trop = 2; else

%--- Update equations ----------------------------------------- rho2 = (X(1, i) - pos(1))^2 + (X(2, i) - pos(2))^2 + ... (X(3, i) - pos(3))^2;

traveltime = sqrt(rho2) / settings.c ;%补充code

%--- Correct satellite position (do to earth rotation) --------

Rot_X = e_r_corr(traveltime, X(:, i));%修正第i颗卫星的位置

%求仰角，方位角（有函数代码）

[az(i), el(i), dist] = topocent(pos(1:3, :), Rot_X - pos(1:3, :));

if (settings.useTropCorr == 1)

%--- Calculate tropospheric correction -------------------- trop = tropo(sin(el(i) * dtr), ...

0.0, 1013.0, 293.0, 50.0, 0.0, 0.0, 0.0); else

% Do not calculate or apply the tropospheric corrections trop = 0; end

end % if iter == 1 ... ... else

% 伪距的修正 补充

omc(i) = (obs(i) - norm(Rot_X - pos(1:3), 'fro') - pos(4) - trop);

% 建立A设计矩阵 补充

A(i, :) = [ (-(Rot_X(1) - pos(1))) / obs(i) ... ];

end % for i = 1:nmbOfSatellites

%判断A的稚是否为4，否则方程无解 补充 if rank(A) ~= 4

pos = zeros(1, 4); return end

x =

% 位置迭代变换更新--补充 pos = end

pos = pos';

%计算精度因子--见原理 if nargout == 4

dop = zeros(1, 5); Q = inv(A'*A);

dop(1) = sqrt(trace(Q)); % GDOP dop(2) = % PDOP dop(3) = % HDOP dop(4) = sqrt(Q(3,3)); % VDOP dop(5) = sqrt(Q(4,4)); % TDOP end

2.实验数据说明：

（1）伪距.mat 和 卫星坐标.mat 是用于测试的卫星文件和相应的卫星对应的观测伪距 （2）结果在实验12的说明文档中。 3.实验要求：

（1）结合流程图和求解原理补充完整函数代码，然后设计实验，将实验文件中所给的数据载入作为上输入；

（2）验证实验结果是否和文件说明中的结果是一致的。

本文来源：https://www.bwwdw.com/article/hwy8.html