最近为了熟悉GNSS软件接收机的设计流程，打好软件接收机的基础，依据《软件定义的GPS和伽利略接收机》书后所附程序（提取码：1111），独自编写基于MATLAB的GNSS_SDR程序（这是可以跑的完整代码，其中部分重要程序是自己写的）（提取码：1111这是可以跑的完整代码，其中部分重要程序是自己写的）。
如果失效请及时联系我，jldxwsj@163.com
----------------------------------
更新：由于最近我的代码总失效，故将我的代码上传至了阿里OSS保存，单击下方链接可直接下载，给大家带来的不便，敬请谅解！
https://jldxwsj.oss-cn-hangzhou.aliyuncs.com/GNSS_SDR/wsj.rar
----------------------------------
前段时间我也使用CSDN，记录了我阅读GNSS_SDR程序的全部历程，对GNSS_SDR程序有了初步的理解和认知，接下来，我将从零开始，尽量亲手编写一个可以跑通的GNSS_SDR程序，并希望借助CSDN平台，记录下完整的学习过程。
本程序是基于MATLAB2020b版本编写的，相较于书后所附MATLAB版本有所不同，因此部分代码有出入，但是影响不大。
由于本人是初学者，因此对matlab软件接收机的理解非常有限，希望通过此次程序编写对接收机的认识有所深入，文中如有错误，欢迎大家在评论区批评指正，我会定时更新本博客。
那么进入正题，实现matlab软件接收机，整体上可以分为三大部分：射频(RF)前端处理、基带数字信号处理（GPS卫星信号的捕获与跟踪）、定位导航运算（GPS的定位和滤波）三大部分。为了能对GPS进行定位，首先需要分析射频前端处理。
目录
射频前端处理与数据集
射频前端模块
射频前端处理模块通过天线接收所有可见GPS卫星的信号，经前置滤波器和前置放大器的滤波放大后，再与本机振荡器产生的正弦波本振信号进行混频而下变频成中频(IF)信号，最后经模数(A/D)转换器将中频信号转变成离散时间的数字中频信号。
由于卫星与接收机之间的相对运动会引起信号载波频率的多普勒效应，使接收到的卫星信号的载波频率发生偏移。因为这种相对运动状况和相应的多普勒频移量通常是不可预测的，所以射频前端只得将接收信号从射频下变频到中频(或者近基带)，而不是直接下变频到真正的基带。同时，(考虑到数字电路和数字信号处理的优越性，接收机尽可能将对中频信号的处理安排到A/D转换之后进行。
图1 射频前端处理过程
数据集GNSS_signal_records
本文所采用的数据集为经过射频前端处理后的数据，其中：
“gpsdata - discrete teccomponents -fs38_192-if9_55.bin”（提取码：1111），是由美国科罗拉多大学博尔德分校收集的。采样频率:38.192 MHz/中频:9.55 MHz(标称)
“GPS_and_GIOVE_A-NN-fs16_3676-if4_1304.bin”（提取码：1111），它包含了通常的GPS数据，但在当时伽利略GIOVE-A卫星也是可见的，并以1575.42 MHz的频率发射。该数据集由Marco Pini和都灵理工大学提供。
本程序主要依据上述两个数据集展开定位，接下来将逐一编写接收机代码，编写的顺序主要按照：基带数字信号处理（GPS卫星信号的捕获与跟踪）、定位导航运算（GPS的定位和滤波）的顺序，用到哪个函数编写哪个程序，首先从init.m开始。
init.m程序的编写
包含内容：
1、将命令行窗口中的输出显示格式更改为指定的格式（long、compact隐藏过多的空白行以便在一个屏幕上显示更多输出）。2、添加MATLAB路径
3、在命令行显示启动（欢迎）语句
4、初始化参数，此处使用initSettings()函数
5、在命令行输出“数据处理中，请稍后……”字样，并将参量传入probeData(settings)函数中处理，如果文件有误，返回“文件格式错误”，并返回错误类型。
6、在命令行显示“原始中频数据如图所示”，“如需修改参数，请在"initSettings.m"中配置”
7、在命令行确认是否需要运行，输入“1”开始数据处理（调用postProcessing），输入“0”退出。
编写完成后运行结果如图2所示：
图2 init.m的运行结果图
附程序代码（提取码：1111）【请注意，这是可以跑的本博客的完整代码，读者跑之前需要把数据集下载到GNSS_signal_records文件夹下，该代码与文章开头的代码是一个链接，不需要重复下载】：
% GNSS软件接收机程序编写
% 作者：王少靖
% 2021/10/23
%——————————————————————
%clear all
clc
%% 将命令行窗口中的输出显示格式更改为指定的格式
format compact
format long
%% 添加MATLAB路径
addpath D:\GNSS_SDR\wsj\GNSS_SDR\include
addpath D:\GNSS_SDR\wsj\GNSS_SDR\geoFunctions
%% 在命令行显示启动（欢迎）语句
disp("欢迎使用GNSS软件接收机")
disp("程序由wsj编写")
%% 初始化参数
settings=initSettings();
%% 处理数据
try
disp("数据处理中，请稍后……")
probeData(settings);
catch
disp("文件处理异常")
rethrow(ME.identifier)
end
%% 显示处理结果
disp("“原始中频数据如图所示”")
disp("“如需修改参数，请在'initSettings.m'中配置”)")
%% 确认是否运行
startchar=input('输入“1”开始数据处理，输入其他键退出程序：');
if startchar==1
postProcessing
endinitsettings.m程序编写
Settings中包含以下参数，总结如下：
数据处理过程的参数设置
settings.msToProcess 要处理的毫秒数使用36000 +任何瞬态(见下面-在导航参数)，以确保导航子帧被提供
settings.numberOfChannels 用于信号处理的通道数
settings.skipNumberOfBytes 移动处理的起点，可用于在数据记录的任何点开始信号处理(例如长记录)。
原始信号文件名和其他参数设置
settings.fileName 数据文件(信号记录)的名称
settings.dataType 用于存储一个样本的数据类
settings.IF %中频[Hz]
settings.samplingFreq %采样频率[Hz]
settings.codeFreqBasis %C/A码的码率[Hz]
settings.codeLength %码片长度
捕获回路的参数设置
settings.skipAcquisition %是否跳过捕获程序，如果置1则在postProcessing.m中跳过捕获程序
settings.acqSatelliteList %需要捕获的卫星名单，为了加快捕获速度，可以排除一些卫星
settings.acqSearchBand %[kHz]最大多普勒频移的估算过程
settings.acqThreshold %阈值信号的确定
跟踪回路的参数设置
settings.dllDampingRatio %衰减率 %C/A码跟踪回路参数
settings.dllNoiseBandwidth %噪声带宽[Hz]
settings.dllCorrelatorSpacing %相关器间距[chips]
settings.pllDampingRatio %衰减率 %载波跟踪环路参数
settings.pllNoiseBandwidth %噪声带宽[Hz]
导航定位方式的选择
settings.navSolPeriod %计算伪距和位置的周期[ms]
settings.elevationMask %去除低仰角的卫星[degrees 0 - 90] 初始化取10
settings.useTropCorr %启用/禁用对流层校正0 - Off；1 - On
settings.truePosition.E = nan; %UTM系统（UNIVERSAL TRANSVERSE MERCARTOR GRID）
settings.truePosition.N = nan; %中天线的真实位置(如果已知)，否则输入的所有NaN和平均位置将被用作参考。
settings.truePosition.U = nan;
绘图参数设置
settings.plotTracking %启用/禁用每个通道的跟踪结果绘图% 0 - Off
settings.c %光速299792458
settings.startOffset %[ms] Initial sign. travel time
附上本段程序代码（提取码：1111）
本段程序因为是参数设置，因此未作修改，参考原程序。
probeData.m程序编写
通过阅读程序明确，该程序的主要目的是通过该函数，绘制原始数据信息:时域图，频域图和直方图。
包含内容：
1、检测传入变量个数，如果是1个，就将输入参量varargin赋给settings，如果是2个，就分别赋给fileNameStr, settings。varargin是一个cell型，需要使用varargin{1}赋给settings，settings才是结构体，因此settings支持点索引和赋值。如果传入参量不是1或2，那么报错“函数参数数量不正确”。
2、以只允许读写的形式打开二进制文件“gpsdata - discrete teccomponents -fs38_192-if9_55.bin”文件（fileNameStr），判断是否打开成功，如果失败，报错“数据文件打开失败”。
3、设置二进制文件的读取方式（fseek），从头开始，是否跳过数据头（skipNumberOfBytes），计算每一个码周期的样本数量samplesPerCode。
4、读取10个码周期，并以settings中指定的形式存入data中，同时保存样本数量count，同时关闭文件
5、做出判断，判断文件长度是否小于10个码周期，如果是，返回error，“数据文件过小，没有足够的数据”
6、画出时域波形图，按照采样频率进行采样，参数设置为每50个点显示一个
7、使用pwelch()函数画出功率谱密度图
8、用hist()函数绘制柱状图
根据以上内容，画出以下图形：（放上对比图，原程序和自己编写的程序对比）
图3 绘制原始数据信息 (对比图)
左图为原始程序绘制的图形，右图为自己编写的绘图程序，其中左下角的图形是存在问题的，问题的原因是：hist函数没有添加横坐标间距设置，间距设置需要合理，否则会出现图三的情况。
编程过程中，遇到了一个问题，在画功率谱密度函数（图二） 时，参数的设置不很清晰：
pwelch(data-mean(data),16384, 1024, 2048, settings.samplingFreq/1e6)
此语句如果只写成pwelch(data-mean(data))，则计算出的频谱的效果比较差，以下是完整的程序（提取码：1111）：
% 该函数主要任务：绘制原始数据信息:时域图，频域图和直方图。
function probeData(varargin)
% 判断传入参数个数
if nargin==1
settings=varargin{1};
fileNameStr=settings.fileName;
elseif nargin==2
fileNameStr=varargin{1};
settings=varargin{2};
else
error("函数参数不正确")
end
fid=fopen(fileNameStr,'rb');
% 设置读数据格式
if fid>0
fseek(fid,settings.skipNumberOfBytes,-1);
%计算一个码周期样本数
samplesPerCode=settings.samplingFreq/...
(settings.codeFreqBasis/settings.codeLength);
[data,count]=fread(fid,[1,10*samplesPerCode],settings.dataType);
fclose(fid);
if count<10*samplesPerCode
error("数据文件过小，没有足够的数据")
end
T=0:1/settings.samplingFreq:1e-3;
figure(101)
subplot(2,2,1)
plot(T(1:round(samplesPerCode/50)),data(1:round(samplesPerCode/50)))
%每50个点采样一次
grid on;
title ('时域分布情况');
xlabel('时间/微秒');
ylabel('幅值');
subplot(2,2,2)
pwelch(data-mean(data),16384, 1024, 2048, settings.samplingFreq/1e6)
%此处未弄清楚参数为何这样设置？
axis tight;
grid on;
title ('频谱分布情况');
xlabel('频率(MHz)');
ylabel('幅值');
subplot(2,2,3)
hist(data) %此处需要设置间距，间距设置需要合理，否则会出现图三的情况
axis tight;
xmax = max(abs(data)) + 1;
%xlim([-xmax xmax]);
adata = axis;
axis([-xmax xmax adata(3) adata(4)]);
title ('Histogram');
xlabel('Bin'); ylabel('Number in bin');
subplot(2, 2, 4);
hist(data,-10:10)
dmax = max(abs(data)) + 1;
axis tight;
adata = axis;
axis([-dmax dmax adata(3) adata(4)]);
grid on;
title ('Histogram');
xlabel('Bin'); ylabel('Number in bin');
else
error("函数参数数量不正确");
end以上模块是在捕获前的处理程序，下一节将主要记录捕获函数的编写思路。
acquisition.m 捕获函数编写
该函数的功能是：对收集到的“数据”执行冷启动采集，并搜索所有卫星的GPS信号，在参数结构体中的“acqsatellite elist”字段中列出。同时返回一个“acqResults”结构体，其中保存了检测信号的码相位和频率，其相关原理详见基于MATLAB编写的GNSS_SDR（GNSS软件接收机）——自学笔记（2）
程序主要包含以下步骤：
1、输入两个参量，longSignal（长度是11ms的中频信号）和settings，返回一个acqResults结构体，将检测信号的码相位和频率保存下来
2、计算每1码片周期内的采样点数samplesPerCode
3、创建两个1ms的相关数据矩阵（signal1和signal2）和一个零直流数据signal0DC
4、计算采样时间ts，产生相位点phasePoints（不含频率信息），设置需要搜索的频带数目numberOfFrqBins
5、调用makeCaTable（settings）函数，产生caCodesTable码
6、构建一个二维矩阵results，并赋初值，包含所有搜索的频率带宽和每次相关的结果（numberOfFrqBins, samplesPerCode）；构建一个一维矩阵frqBins，并赋初值，包含numberOfFrqBins列
7、初始化捕获结果acqResults.carrFreq、acqResults.codePhase、acqResults.peakMetric分别初始化：检测信号的载波频率、检测信号的C/A码相位、被检测信号的相关峰值比
8、执行搜索所有列出的PRN卫星
9、对C/A码发生器求傅里叶变换，再求共轭得到caCodeFreqDom，然后对每个频带进行相关分析：计算出需要检测的频率frqBins(frqBinIndex)，[中频-最大多普勒频移，中频+最大多普勒频移]，步长为0.5e3
10、产生正弦sinCarr、余弦表cosCarr（频率*相位点phasePoints），并分别与signal1和signal2相乘得到I1、I2、Q1、Q2，然后分别将signal1和signal2变换到频域得到IQfreqDom1和IQfreqDom2，接着分别与经过傅里叶变换、复数共轭的C/A码相乘，得到convCodeIQ1和convCodeIQ2，最后进行傅里叶反变换，取模得到acqRes1和acqRes2
11、比较acqRes1和acqRes2中的最大值，并将结果赋给results
12、求解第一相关峰和第二相关峰，最大相关峰的[peakSize frequencyBinIndex]和[peakSize codePhase]，第一相关峰好求，第二相关峰需要判断一下最大相关是否在边界上，如果是，就将改变相位取值范围codePhaseRange，然后去掉第一相关峰的那组值，求最大值，即为第二相关峰值secondPeakSize。
13、储存第PRN颗星的判决结果acqResults.peakMetric，并检验peakSize/secondPeakSize是否达到门限条件，如果大于门限判别条件，就增加相干长度，准确判断频率值。
14、产生一组10ms的C/A码，具体过程：1.使用generateCAcode（PRN）函数产生caCode，2.确定10ms内，每个码片的序号codeValueIndex，大概是每个码片有38个采样样本，因此前38个序号是0，第39-76个序号就是1，以此类推；3.将其拼接成一个长C/A码longCaCode（10个周期）
15、计算第1个codePhase处的零直流数据与长C/A码相乘，得到xCarrier
16、以下几个语句不清楚其思路附在下面，随后解答（已解答，见文末）
%Find the next highest power of two and increase by 8x fftNumPts = 8*(2^(nextpow2(length(xCarrier)))); %--- Compute the magnitude of the FFT, find maximum and the %associated carrier frequency fftxc = abs(fft(xCarrier, fftNumPts)); uniqFftPts = ceil((fftNumPts + 1) / 2); [fftMax, fftMaxIndex] = max(fftxc(5 : uniqFftPts-5)); fftFreqBins = (0 : uniqFftPts-1) * settings.samplingFreq/fftNumPts; %--- Save properties of the detected satellite signal ------------- acqResults.carrFreq(PRN) = fftFreqBins(fftMaxIndex);该程序的主要思想是：采用10ms的积分时间，有利于提高信噪比，其相邻两个10ms的相干积分内必定有一个没有发生数据跳变，因此，对其所求结果求最大值，确定频率准确值。
17、输出捕获的频率值和相位值acqResults.carrFreq(PRN)和acqResults.codePhase(PRN)
原程序与本程序的对比结果判断，该程序已经实现了捕获的功能。
图4 捕获函数对比图
可以发现，自行编写的程序在 输出上与原程序无差别，但在编写上存在一定差异和想法，详见程序和注释（提取码：1111）。
%输入两个参量，longSignal（长度是11ms的中频信号）和settings，
%返回一个acqResults结构体，将检测信号的码相位和频率保存下来
function acqResults=acquisition(longSignal,settings)
%计算每1码片周期内的采样点数samplesPerCode
samplesPerCode=settings.samplingFreq/...
(settings.codeFreqBasis/settings.codeLength);
%创建两个1ms的相关数据矩阵（signal1和signal2）
%和一个零直流数据signal0DC
signal1=longSignal(1:samplesPerCode);
signal2=longSignal(samplesPerCode+1:2*samplesPerCode);
signal0DC=longSignal-mean(longSignal);
%计算采样时间ts，产生相位点phasePoints（不含频率信息），
%设置需要搜索的频带数目numberOfFrqBins
ts=1/settings.samplingFreq;
phasePoints=(0:samplesPerCode-1)*2*pi*ts;
numberOfFrqBins=floor(settings.acqSearchBand)*2+1;
%调用makeCaTable（settings）函数，产生caCodesTable码
caCodesTable=makeCaTable(settings);
%构建一个二维矩阵results，并赋初值，包含所有搜索的频率带宽
%和每次相关的结果（numberOfFrqBins, samplesPerCode）；
%构建一个一维矩阵frqBins，并赋初值，包含numberOfFrqBins列
results=zeros(numberOfFrqBins, samplesPerCode);
frqBins=zeros(1,numberOfFrqBins);
%初始化捕获结果acqResults.carrFreq、acqResults.codePhase、
%acqResults.peakMetric分别初始化：检测信号的载波频率、
%检测信号的C/A码相位、被检测信号的相关峰值比
acqResults.carrFreq=zeros(1,length(settings.acqSatelliteList));
acqResults.codePhase=zeros(1,length(settings.acqSatelliteList));
acqResults.peakMetric=zeros(1,length(settings.acqSatelliteList));
fprintf('(');
%执行搜索所有列出的PRN卫星
for PRN = settings.acqSatelliteList
% 对C/A码发生器求傅里叶变换，再求共轭，然后对每个频带进行相关分析：
% 计算出需要检测的频率，
% [中频-最大多普勒频移，中频+最大多普勒频移]，步长为0.5e3
caCodeFreqDom=conj(fft(caCodesTable(PRN, :)));
%调试过程修改后，请注意caCodesTable(PRN, :)只写caCodesTable不可以！！
for frqBinIndex=1:numberOfFrqBins
frqBins(frqBinIndex)=settings.IF...
-(settings.acqSearchBand/2) * 1000+...
0.5e3 * (frqBinIndex - 1);
%产生正弦sinCarr、余弦表cosCarr（频率*相位点phasePoints），
%并分别与signal1和signal2相乘得到I1、I2、Q1、Q2，
%然后分别将signal1和signal2变换到频域得到IQfreqDom1和
%IQfreqDom2，接着分别与经过傅里叶变换、复数共轭的C/A码相乘，
%得到convCodeIQ1和convCodeIQ2，最后进行傅里叶反变换，
%取模得到acqRes1和acqRes2
sinCarr=sin(frqBins(frqBinIndex)*phasePoints);
cosCarr=cos(frqBins(frqBinIndex)*phasePoints);
I1=sinCarr.*signal1;
I2=sinCarr.*signal2;
Q1=cosCarr.*signal1;
Q2=cosCarr.*signal2;
IQfreqDom1=fft(I1+j*Q1);
IQfreqDom2=fft(I2+j*Q2);
% 调试过程：注意I和Q的位置
convCodeIQ1=IQfreqDom1.*caCodeFreqDom;
convCodeIQ2=IQfreqDom2.*caCodeFreqDom;
acqRes1=abs(ifft(convCodeIQ1)) .^ 2;
acqRes2=abs(ifft(convCodeIQ2)) .^ 2;
%破案了，这里的问题，这里需要取模后平方？似乎与书上不太一致，但是阈值也会相应改变
%比较acqRes1和acqRes2中的最大值，并将结果赋给results
% results(frqBinIndex,:)=max(max(acqRes1,acqRes2));
if (max(acqRes1) > max(acqRes2))
results(frqBinIndex, :) = acqRes1;
else
results(frqBinIndex, :) = acqRes2;
end
%请注意：此处需要用if/else语句，如果用max(max())那么codePhase求出来是错的
end
[peakSize frequencyBinIndex]=max(max(results,[],2));
[peakSize codePhase]=max(max(results));
% samplesPerCodeChip = round(settings.samplingFreq / settings.codeFreqBasis);
samplesPerCodeChip = round(settings.samplingFreq / settings.codeFreqBasis);
excludeRangeIndex1 = codePhase - samplesPerCodeChip;
excludeRangeIndex2 = codePhase + samplesPerCodeChip;
%--- Correct C/A code phase exclude range if the range includes array
%boundaries 如果取值在边界附近，需要修正C/A码相位
if excludeRangeIndex1 < 2
codePhaseRange = excludeRangeIndex2 : ...
(samplesPerCode + excludeRangeIndex1);
elseif excludeRangeIndex2 >= samplesPerCode
codePhaseRange = (excludeRangeIndex2 - samplesPerCode) : ...
excludeRangeIndex1;
else
codePhaseRange = [1:excludeRangeIndex1, ...
excludeRangeIndex2 : samplesPerCode];
end
% if codePhase-samplesPerCodeChip<2
% codePhaseRange=codePhase+samplesPerCodeChip:...
% (samplesPerCode + codePhase - samplesPerCodeChip);
% elseif codePhase+samplesPerCodeChip>samplesPerCode
% codePhaseRange=(codePhase+samplesPerCodeChip - samplesPerCode):...
% codePhase-samplesPerCodeChip;
% else
% codePhaseRange=[1:codePhase-samplesPerCodeChip,...
% codePhase+samplesPerCodeChip:samplesPerCode];
% end
secondPeakSize=max(max(results(frequencyBinIndex,codePhaseRange)));
%检验peakSize/secondPeakSize是否达到门限条件，
%如果大于门限判别条件，就增加相干长度，准确判断频率值。
acqResults.peakMetric(PRN) = peakSize/secondPeakSize;
if (peakSize/secondPeakSize)>settings.acqThreshold
fprintf('%02d ', PRN);
%产生一组10ms的C/A码，具体过程：
% 1.使用generateCAcode（PRN）函数产生caCode，
% 2.确定10ms内，每个码片的序号codeValueIndex，大概是每个码片有38个采样样本，
% 因此前38个序号是0，第39-76个序号就是1，以此类推；
% 3.将其拼接成一个长C/A码longCaCode（10个周期）
caCode=generateCAcode(PRN);
% codeValueIndex=floor(settings.samplingFreq*...
% ((1:10*samplesPerCode))/settings.codeFreqBasis);
% 出错了，重写
codeValueIndex = floor((ts * (1:10*samplesPerCode)) / ...
(1/settings.codeFreqBasis));
longCaCode=caCode(rem(codeValueIndex,settings.codeLength)+1);
xCarrier=longCaCode.*signal0DC(codePhase:(codePhase+10*samplesPerCode-1));
%此处代码不理解
fftNumPts = 8*(2^(nextpow2(length(xCarrier))));
fftxc = abs(fft(xCarrier, fftNumPts));
uniqFftPts = ceil((fftNumPts + 1) / 2);
[fftMax, fftMaxIndex] = max(fftxc(5 : uniqFftPts-5));
fftFreqBins = (0 : uniqFftPts-1) * settings.samplingFreq/fftNumPts;
%输出捕获的频率值和相位值
acqResults.carrFreq(PRN) = fftFreqBins(fftMaxIndex);
acqResults.codePhase(PRN) = codePhase;
else
fprintf('. ');
end
end
fprintf(')\n');
plotAcquisition.m绘制捕获函数图形
绘制捕获函数的结果图，程序比较简单，实现了函数绘制了条形结果图，且没有显示未在清单上的卫星图。
主要内容：
1、输入参量acqResults，设置画布编号
2、将相关第一峰值和第二峰值之比使用条形图做出，acqResults.peakMetric，并设置标题“捕获结果”，横坐标“PRN编号”，纵坐标“比值”，设置坐标范围及图形格式
3、标记被捕获的信号为绿色，并添加图例
结果如图所示：
图5 捕获卫星与未捕获卫星的对比图
附上程序（提取码：1111）：
function plotAcquisition(acqResults)
figure(102)
bar(acqResults.peakMetric);
hold on
a=acqResults.peakMetric .* (acqResults.carrFreq > 0);
bar(a,'g');
legend('未被捕获的信号', '被捕获信号');
title("捕获结果");
xlabel('PRN编号');
ylabel('相关第一峰值和第二峰值之比');针对上文提出的疑问，在此做一次更新：
1. 对pwelch函数的理解问题：由于pwelch函数是利用Welch'平均功率图法进行的绘图。为了得到更好的功率谱曲线，Welch法允许数据段重叠，即在对x_N (n)分段时，可允许每一段的数据有部分的交叠。此外，在计算周期图之前对数据段进行加窗，每一段的数据窗口可以不是矩形窗口（默认为汉明窗），这样可以改善由于矩形窗所产生的谱失真。
通过查阅pwelch函数的解释，可以看出，该段语句使用[pxx,f] = pwelch(x, window, noverlap, f, fs)的形式，其中window默认选择汉明窗，pwelch函数将x分割成与window长度相等的重叠段，然后将每个信号段与window中指定的这一部分相乘。Noverlap是指重叠样本数，指定为小于窗口长度的正整数。fs采样率，采样率是单位时间内的采样数。如果时间的单位是秒，那么采样率的单位是Hz。因为转换为MHz，所以settings. samplingFreq需要除以106。
2. 在acquisition.m 捕获函数中，因为settings.acqSearchBand是指需要搜索的频率范围，而不是指多普勒频移，因此，需要除以2。
3. 对频率值进行精细化使用的是对FFT进行补零的方法，可以通过补零观察到更多的频点，因此可以找到更准确的频率点，但是这并不意味着补零能够提高真正的频谱分辨率。具体内容，可以查看论文《GPS软件接收机捕获方法研究》
扫一扫
413
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
