总体介绍
为了加深对gps软件接收机的理解，首先阅读了基于matlab语言编写的gps软件接收机程序，最近接触到了python语言，所以现在拟用python语言复现一部分，同时在这里记录自己的学习过程，目前实现了从文件中读取数据、ca码生成，以及捕获方法的实现，接下来主要介绍这三方面的内容。由于是python语言初学者，所以代码风格这方面还有很多欠缺。
settings到self
基于matlab实现的gps软件接收机用一个settings结构体定义了接收机自身的属性，包括采样频率（samplingfreq）伪码码率（codeFreqBasis）伪码长度（codelength）等等属性，python面向对象编程，所以我将接收机作为一个对象，用self定义的对象所具有的属性，目前只定义了实现捕获方法所需要的几个属性，（后续会逐渐更新）
def __init__(self):
self.filepath = 'gnsa14.bin' # 文件路径
self.skipNumberOfBytes = 0 # 移动处理的节点
self.samplingFreq = 16.3676e6 # 采样频率
self.codeFreqBasis = 1.023e6 # 伪码码率
self.codeLength = 1023 # 伪码长度
self.dataType = 'int8' # 数据类型
self.acqsearchband = 14 # (khz) 最大多普勒频移的估算过程
self.samplesPerCode = round(
self.samplingFreq / (self.codeFreqBasis / self.codeLength)) # 一个ca码周期（1023个码片）或者（1ms）的采样点数
self.acqsatellitelist = [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24,
25, 26, 27, 28, 29, 30, 31, 32]
self.IF = 38400 # 中频频率
self.acqThreshold = 2.5 # 阈值上面为编写后续程序所用到的一些属性，命名可能有不规范之处。
读取文件中的数据probodata方法
先放代码和程序实现对比结果，代码如下：
def probodata(self): # 读取数据方法
# 函数拟实现的功能为读取文件数据并做出时域、频域以及直方图
with open(self.filepath, 'rb') as f:
f.seek(self.skipNumberOfBytes,
os.SEEK_CUR) # seek方法用于移动文件指针到指定位置，当前语句表示移动文件指针到settings.skipNumberOfBytes处（seek_cur表示从当前位置开始移动）
samplesPerCode = round(
self.samplingFreq / (self.codeFreqBasis / self.codeLength)) # 计算每个码周期的采样数
data = np.fromfile(f, dtype=self.dataType, count=11 * samplesPerCode)
count = len(data)
try:
if count < 11 * samplesPerCode:
# The file is too short
raise ValueError("Could not read enough data from the data file.")
except ValueError as e:
print(e)
t = np.arange(0, 1e-3, 1 / self.samplingFreq)
plt.figure(1)
plt.subplot(2, 2, 1)
plt.plot(t[0:round(samplesPerCode / 50)], data[0:round(samplesPerCode / 50)])
plt.axis('tight')
plt.grid(True)
plt.title('Time domain')
plt.xlabel('Time [s]')
plt.ylabel('Amplitude')
plt.subplot(2, 2, 2)
f, pxx = welch(data - np.mean(data), fs=self.samplingFreq / 1e6, nperseg=1024, noverlap=256, nfft=2048)
plt.semilogy(f, pxx)
plt.xlabel('Frequency [MHz]')
plt.ylabel('magtitude')
plt.title('Frequency domain')
plt.subplot(2, 2, 3)
plt.hist(data, bins=np.arange(-10, 11))
dmax = max(abs(data)) + 1
plt.axis('tight')
adata = plt.axis()
plt.axis([-dmax, dmax, adata[2], adata[3]])
plt.grid(True)
plt.title('Histogram')
plt.xlabel('Bin')
plt.ylabel('Numberinbin')
# 显示图形
plt.show()matlab程序实现结果和以上程序实现结果对比如下：
用python实现的程序在画直方图的时候区间划分有些问题，没有matlab程序更加直观，后续会尝试改进源程序，频域图和时域图形基本一致，主要遇到的问题有以下几个方面，首先是分析功率谱密度方法welch（）的使用，使用的是python内置方法，它的主要输入参数有以下几个：（从科大讯飞星火大语言模型得到的回答：）
welch函数的参数主要包括以下几种：
- x：这是输入的时间序列数据，通常是一维数组形式。
- fs：表示采样频率，是一个浮点数，默认值为1.0。这个参数用于指示每个时间序列的采样率。
- window：这个参数是窗口函数，可以是一个字符串，如'hamming'、'hanning'、'bartlett'、'blackman'、'flattop'等，也可以是自定义的窗口函数数组。窗口函数主要用于减少频谱泄漏和由于信号不连续性或周期性引起的误差。
- nperseg：表示每个段的长度，是一个整数。这个参数决定了每个窗口处理的数据长度。
- noverlap：表示窗口之间的重叠长度。通常取值为33%到50%。重叠得越多，得到的图像越平滑；反之，则更粗糙。
- NFFT：表示FFT（快速傅里叶变换）数据点的个数，它决定了频域分辨率。通常设置为大于每一段点数的最小2次幂，以获得最高的频域分辨率。
这些参数最后选择2的n的幂效果较佳，nperseg参数必须大于noverlap，还查询了npseseg参数的选取原则：
nperseg 参数是用于指定每个段的长度，在计算功率谱密度（PSD）时使用。这个参数对于频谱分析的精度和分辨率有重要影响。
-
长度选择：
nperseg的值通常取决于信号的特性和你想要捕捉的频率成分。- 较短的
nperseg值会提供较高的频率分辨率，但可能导致估计的方差增加，因为每个段的数据点较少。 - 较长的
nperseg值会提供更平滑的功率谱密度估计，但频率分辨率会降低。
-
与采样率的关系：
nperseg的选择还应该考虑采样率 (fs)。例如，如果你的信号是 1000 Hz 采样的，而你想要分辨 10 Hz 的成分，你需要至少大约 50 个样本点来清晰分辨该频率成分（根据 Nyquist 准则）。
-
常用规则：
- 一个常用的经验公式是将
nperseg设置为你想要检测的最小频率的两倍的倒数，再乘以采样率。例如，如果你想检测最低为 10 Hz 的频率，那么nperseg应至少为1/（10 * 2） = 0.05秒的数据，如果采样率为 1000 Hz，则nperseg应至少为0.05 * 1000 = 50个样本点。
- 一个常用的经验公式是将
- 另外两个问题是np.fromfile方法的使用和datatype的选取，matlab源程序使用的是’schar‘或者int8，这里使用的是int8类型，由于遇到的一个小问题就是plt.show和plt.show(),加不加括号，是两种完全不同的效果，一开始没加括号导致程序直接运行结束，调试了很长时间，发现是这个括号没加，
ca码生成涉及到的两个方法：generateCaCode和maketables
首先介绍generateCaCode，源程序如下：
def generateCaCode(self, PRN):
##该函数的作用是产生ca码
##返回值为一维ca码数组，长度为1023
g2_yiwei = np.array([5, 6, 7, 8, 17, 18, 139, 140, 141, 251,
252, 254, 255, 256, 257, 258, 469, 470, 471, 472,
473, 474, 509, 512, 513, 514, 515, 516, 859, 860,
861, 862]) # 相位选择器对应于ca码的平移
##产生g1码，方法为p24图2.9
#
g1 = np.zeros(1023, dtype=int)
reg = -1 * np.ones(10)
for j in range(0, 1023):
g1[j] = reg[9]
savebit = reg[2] * reg[9]
reg[1: 10] = reg[0: 9]
reg[0] = savebit
# 产生g2码
n = g2_yiwei[PRN - 1]
g3 = np.zeros(1023, dtype=int)
# print(type(g3))
reg = -1 * np.ones(10)
# print()
for i in range(0, 1023):
g3[i] = reg[9]
savebit = reg[1] * reg[2] * reg[5] * reg[7] * reg[8] * reg[9]
reg[1: 10] = reg[0: 9]
reg[0] = savebit
g3 = np.array(list(g3[1023 - n:]) + list(g3[0:1023 - n]))
cacode = -1 * g1 * g3
return cacode程序实现原理如下图：
将ca码发生器结构编程实现，首先对于异或运算，程序是采用相乘的方法间接实现的：实现原理如下：
0设置为1，1设置为（-1）相乘实现异或运算，最后可以通过乘-1转换过来。
第二个问题是G2i码中相位选择器的实现，相位选择器实现有两种方式，一个是将原ca码平移得到，代码中g2-yiwei列表中的元素代码要移动的位数，例如，对于卫星编号PRN=1发射的ca码，需要将原ca码移动5位得到，另一种方法是直接进行异或运算，例如对于卫星编号PRN=1，它对应的G2i码是寄存器2和寄存器6经过异或运算得到（见下图）的，具体的ca码生成表如下：
下面展示几个运行结果：
代码生成了PRN=1、6、18、27的前十位ca码，其中1代表1，-1代表0.，转变位对应的8进制码分别为：1440、1455、1467和1770（最后一个被水印遮挡了）可以看出与上图的前十码是吻合的。
接下来时maketables方法的实现：源代码如下：
def makecatable(self):
# 该方法的作用是产生ca码
cacodestable = np.zeros((32, self.samplesPerCode))
# 每个ca码周期采样16368个点，所以每个码片周期采样点数为16368/1023 = 16个点 所以有generatecode产生的ca码每个数字要重复16此
ts = 1 / self.samplingFreq
tc = 1 / self.codeFreqBasis
for PRN in range(1, 33):
# 对于给定的卫星编号PRN产生相应的ca码
cacode = self.generateCaCode(PRN)
codevalueindex = np.floor(np.arange(0, self.samplesPerCode, 1) / tc * ts)
# print(codevalueindex[16200:16368])
for i in codevalueindex:
cacodestable[PRN - 1, int(i * 16):int((i + 1) * 16)] = cacode[int(i)]
return cacodestable该方法实现的功能是生成一个包含32颗gps卫星的ca码表，输出是一个二维数组，共32行，但此ca码表是经过采样的ca码表，每个码周期（1023个码片）采样点数位samplespercode，计算方法在第二部分（settings到self）有注释，换算到每个码片采样16个点，这部分因为运行结果太长，故结果不做展示。（ca码表的生成算法很重要，直接影响后续的捕获结果，按照上面的算法，可以实现对第一份接收信号的捕获，但是对第二份数据的捕获无能为力，算法有些瑕疵，以下是改进的ca码表生成算法，
def makecatable(self):
# 该方法的作用是产生ca码
# cacodestable = np.zeros((32, 1))
cacodestables = []
# 每个ca码周期采样16368个点，所以每个码片周期采样点数为16368/1023 = 16个点 所以有generatecode产生的ca码每个数字要重复16此
ts = 1 / self.samplingFreq
tc = 1 / self.codeFreqBasis
for PRN in range(1, 33):
# 对于给定的卫星编号PRN产生相应的ca码
cacode = self.generateCaCode(PRN)
codevalueindex = np.ceil(ts*np.arange(0, self.samplesPerCode, 1) / tc)
# print(codevalueindex[16200:16368])
cacodestable = []
for i in codevalueindex:
cacodestable.append(cacode[int(i-1)])
cacodestables.append(cacodestable)
tables = np.array(cacodestables)生成的ca码表本质上是一个二维数组，行数是gps卫星数（32），列数是一个码周期采样点数，ca码的码长是1023个码片，每个码片都需要采样一部分的点（大约为tc/ts个），在制作ca码表的时候，一个码片对应ca码表的多个点，在一个码片持续时间内，采样点的幅值都是一样的，算法实现时，需要先生成一个索引列表，索引列表中每tc/ts个的元素是一样的值，然后利用索引在原ca码表中生成相应的值。
捕获方法acquisition的实现
首先是实现效果图对比：
下面的两张图是另一份数据的捕获结果图，与matlab的捕获结果是一样的
程序源代码如下：
def acquisition(self, longsignal):
# 输入两个参量，长度为11ms的中频信号，和类的一些属性
# 返回一个字典，字典中有两个键，分别是检测信号的码相位和频率
# 创建两个1ms的相关数据矩阵（signal1和signal2）
# 和一个零直流数据
signal1 = longsignal[0:self.samplesPerCode]
signal2 = longsignal[self.samplesPerCode:2 * self.samplesPerCode]
signal0DC = longsignal - np.mean(longsignal)
longcacode = np.zeros(10*self.samplesPerCode)
ts = 1 / self.samplingFreq # 计算采样时间ts
tc = 1 / self.codeFreqBasis
phasepoint = np.arange(self.samplesPerCode) * 2 * np.pi * ts # self.samplesPerCode = 16368
numberoffrqbins = int(math.floor(self.acqsearchband) * 2 + 1) # 计算需要搜索的中心频点数目 gps原理与接收机设计p353
cacodestable = self.makecatable() # 调用产生ca码的方法makecatable(self),产生ca码
results = np.zeros((numberoffrqbins, self.samplesPerCode)) # 构建一个二维矩阵results，并赋初值，包含所有搜获的频率带宽
frqbins = np.zeros((1, numberoffrqbins))
acqresult = {'carrfreq': np.zeros((1, 32)), 'codephase': np.zeros((1, 32)),
'peakmetric': np.zeros((1, 32))} # 初始化捕获结果
# #并行码相位搜索捕获算法
print('(')
for PRN in self.acqsatellitelist:
caCodeFreqDom = np.conj(np.fft.fft(cacodestable[PRN-1, :])) # 对ca码求傅里叶变换，在求共轭
# 对整个频带进行相关性分析
for frebinindex in range(1, numberoffrqbins + 1):
frqbins[0, frebinindex - 1] = self.IF - (self.acqsearchband/2) * 1000 + \
0.5e3*(frebinindex-1) # 等间距产生载波（频率间隔为0.5khz）
# 产生本地正弦和余弦载波
sincarr = np.sin(phasepoint * frqbins[0, frebinindex - 1])
coscarr = np.cos(phasepoint * frqbins[0, frebinindex - 1])
# 从信号中分离载波
I1 = sincarr * signal1
Q1 = coscarr * signal1
I2 = sincarr * signal2
Q2 = coscarr * signal2
# 将基带信号转化到频域（通过傅里叶变换） # warning matlab 和 python fft函数对比
IQfreqdom1 = np.fft.fft(I1+Q1*1j)
IQfreqdom2 = np.fft.fft(I2+Q2*1j)
# 频域相乘
IQ1convcode = IQfreqdom1 * caCodeFreqDom
IQ2convcode = IQfreqdom2 * caCodeFreqDom
# 执行fft反变换
acqRes1 = np.abs(np.fft.ifft(IQ1convcode)) ** 2
acqRes2 = np.abs(np.fft.ifft(IQ2convcode)) ** 2
# 求解同一频率单元哪一毫秒的相关值最大，并将其保存
if np.max(acqRes1) > np.max(acqRes2):
results[frebinindex - 1, :] = acqRes1
else:
results[frebinindex - 1, :] = acqRes2
# for frebinindex in range(1, numberoffrqbins + 1):循环结束
# 在结果中寻找相关峰值 对下面的两个求最大值的方法有疑问？我觉得可以直接求出矩阵中所有元素的最大值 其行索引和列索引就是下面代码完成的任务
# 获取每行的最大值及其索引
row_max_values = np.max(results, axis=1) # 返回值row_max_values是个列向量，元素为每一行的最大值
peaksize = np.max(row_max_values)
frequencybinindex = np.argmax(row_max_values) # 求列向量row_max_values中元素最大值的索引
# 获取码相位
col_max_values = np.max(results, axis=0) # 返回值col_max_values是个hang向量，元素为每一lie的最大值
peaksize = np.max(col_max_values)
codephase = np.argmax(col_max_values) # 求列向量row_max_values中元素最大值的索引
samplesPerCodeChip = round(self.samplingFreq / self.codeFreqBasis)
excludeRangeIndex1 = codephase - samplesPerCodeChip
excludeRangeIndex2 = codephase + samplesPerCodeChip
if excludeRangeIndex1 < 2:
codephaserange = np.arange(excludeRangeIndex2, self.samplesPerCode + excludeRangeIndex1+1)
elif excludeRangeIndex2 >= self.samplesPerCode:
codephaserange = np.arange(excludeRangeIndex2-self.samplesPerCode, excludeRangeIndex1+1)
else:
codephaserange = np.concatenate((np.arange(0,excludeRangeIndex1+1) , np.arange(excludeRangeIndex2 , self.samplesPerCode))).astype(int)
secondpeaksize = results[frequencybinindex, codephaserange[0]] # 用于寻找第二峰值，
for i in codephaserange: # 接上行注释，codephaserange[0]后面是否用减去1？？？
if results[frequencybinindex, i] > secondpeaksize:
secondpeaksize = results[frequencybinindex,i]
# 检查是否过门限 储存结果
print(peaksize, secondpeaksize)
acqresult['peakmetric'][0, PRN - 1] = peaksize / secondpeaksize
if peaksize / secondpeaksize > self.acqThreshold:
print(f'{PRN}')
caCode = self.generateCaCode(PRN)
codevalueindex = np.floor(np.arange(0, 10*self.samplesPerCode, 1) / tc * ts)
# codevalueindex[-1] = 1022
for i in codevalueindex:
# cacodestable[PRN - 1, int(i * 38):int((i + 1) * 38)] = cacode[int(i)%1023]
longcacode[int(i * 38):int((i + 1) * 38)] = caCode[int(i)%1023]
# 剥离ca码，采用相乘的方式
xcarr = longcacode*signal0DC[codephase:codephase+10*self.samplesPerCode]
# 通过傅里叶变换求载波频率
fftnumpoint = int(8*nextpow2(len(xcarr)))
# print(type(fftnumpoint))
fftxc = abs(np.fft.fft(xcarr, fftnumpoint));
uniqFftPts = math.ceil((fftnumpoint + 1) / 2);
fftmaxindex = np.argmax(fftxc[5:uniqFftPts-4])
fftfreqbins = np.arange(0, fftnumpoint/2, 1) * (self.samplingFreq / fftnumpoint)
carrfreq = fftfreqbins[fftmaxindex]
acqresult['carrfreq'][0, PRN - 1] = carrfreq
acqresult['codephase'][0, PRN - 1] = codephase
else:
print('-')
print(')\n')
return acqresult
# 绘制捕获结果原matlab程序返回的是一个结构体，包含三个变量，分别为估算的载波频率（carrfreq），码相位（codephase）以及第一第二相关峰比值（peakmetric），本程序返回的是一个字典acqresult，字典的键包含上面三个元素。这段程序的编写主要参考了文章（从零编写基于MATLAB的GNSS_SDR程序（GNSS软件接收机）——学习记录(1)_gnss-sdr-CSDN博客），这篇文章系统介绍了捕获算法的原理，这个博主系统介绍了基于matlab的gps软件接收机程序并给出了详细的中文注释。
总结
首先是用python处理矩阵远没有用matlab处理矩阵的各种运算方便，特别是ca码表的生成方法，其次python对于矩阵的运算需要重点主义矩阵维度的转换，特别是在python中一维数组和和二维数组但是只有一行的情形，运算时很容易报错。数组索引方面从0开始和从1开始有很大的不同，需要非常注意。
其次是画图程序使用还不熟练，特别是对生成的图片添加图例和注释时，出现了一些问题，目前生成的图片只能添加英文注释，中文暂且加不了，还问解决/
还有使用了类的方法和字典的方法来代替matlab程序中的结构体，对面向对象编程有了更进一步的理解，编程过程中，使用了辅助工具星火大语言模型，确实提高了编写效率。
最后就是程序的调试部分，目前我只是使用最笨的调试方式，使用print方式查看某些部分的输出，如果输出的数据很常，很多，这种就很低效，matlab中有变量区，方便了很多，后续需要学习调试方法，特别是对这种代码量很大的程序调试很有必要。
扫一扫
447
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
