GPS信号的捕获(PMF+FFT方法)
目标
- 获取GPS信号的频率偏移和码相位
程序流程
-
读取/生成GPS信号
GPS的导航信号模型为 x ( n T s ) = A ∗ D ( n T s ) C ( n T s ) c o s ( 2 π ( f I F + f d ) n T s + θ ) x(nT_s) = A*D(nT_s)C(nT_s)cos(2\pi(f_{IF} + f_d)nT_s+\theta) x(nTs)=A∗D(nTs)C(nTs)cos(2π(fIF+fd)nTs+θ)
其中, D ( n T s ) D(nT_s) D(nTs) 为导航电文,此程序中只读取了导航电文中的一个bit(1ms); C ( n T s ) C(nT_s) C(nTs)为C/A码; f I F f_{IF} fIF为中频; f d f_d fd为多普勒频移; T s T_s Ts为抽样间隔, F s = 1 / T s F_s=1/T_s Fs=1/Ts为抽样频率
-
混频,得到I路信号和Q路信号
-
基础知识,三角函数和差化积、积化和差
s i n ( α ) ∗ c o s ( β ) = 1 / 2 ∗ ( s i n ( α + β ) + s i n ( α − β ) ) sin(\alpha)*cos(\beta) = 1/2 * (sin(\alpha + \beta) + sin(\alpha - \beta)) sin(α)∗cos(β)=1/2∗(sin(α+β)+sin(α−β))
c o s ( α ) ∗ c o s ( β ) = 1 / 2 ∗ ( c o s ( α + β ) + c o s ( α − β ) ) cos(\alpha)*cos(\beta) = 1/2 * (cos(\alpha + \beta) + cos(\alpha - \beta)) cos(α)∗cos(β)=1/2∗(cos(α+β)+cos(α−β))
s i n ( α ) ∗ s i n ( β ) = 1 / 2 ∗ ( c o s ( α + β ) − c o s ( α − β ) ) sin(\alpha)*sin(\beta) = 1/2 * (cos(\alpha + \beta) - cos(\alpha - \beta)) sin(α)∗sin(β)=1/2∗(cos(α+β)−cos(α−β))
设本地载波信号为 s i n ( 2 π f I F ∗ n T s ) sin(2 \pi f_{IF} * nT_s) sin(2πfIF∗nTs) 和 c o s ( 2 π f I F ∗ n T s ) cos(2 \pi f_{IF} * nT_s) cos(2πfIF∗nTs) ,导航信号与本地载波混频得到I路和Q路混频信号
I路
S I = x ( n T s ) ∗ c o s ( 2 π f I F ∗ n T s ) = A ∗ D ( n T s ) C ( n T s ) c o s ( 2 π ( f I F + f d ) n T s + θ ) ∗ c o s ( 2 π f I F ∗ n T s ) = A ∗ D ( n T s ) C ( n T s ) ∗ 1 / 2 ∗ ( c o s ( 2 π ( 2 f I F + f d ) n T s + θ ) + c o s ( 2 π f d + θ ) ) \begin{align} S_I & = x(nT_s)*cos(2 \pi f_{IF} * nT_s) \\ & = A*D(nT_s)C(nT_s)cos(2\pi(f_{IF} + f_d)nT_s+\theta) * cos(2 \pi f_{IF} * nT_s) \\ & = A*D(nT_s)C(nT_s) * 1/2 * (cos(2\pi(2f_{IF} + f_d)nT_s+\theta) + cos(2\pi f_d + \theta)) \\ \end{align} SI=x(nTs)∗cos(2πfIF∗nTs)=A∗D(nTs)C(nTs)cos(2π(fIF+fd)nTs+θ)∗cos(2πfIF∗nTs)=A∗D(nTs)C(nTs)∗1/2∗(cos(2π(2fIF+fd)nTs+θ)+cos(2πfd+θ))Q路
S Q = x ( n T s ) ∗ s i n ( 2 π f I F ∗ n T s ) = A ∗ D ( n T s ) C ( n T s ) c o s ( 2 π ( f I F + f d ) n T s + θ ) ∗ s i n ( 2 π f I F ∗ n T s ) = A ∗ D ( n T s ) C ( n T s ) ∗ 1 / 2 ∗ ( s i n ( 2 π ( 2 f I F + f d ) n T s + θ ) − s i n ( 2 π f d n T s + θ ) ) \begin{align} S_Q & = x(nT_s)*sin(2 \pi f_{IF} * nT_s) \\ & = A*D(nT_s)C(nT_s)cos(2\pi(f_{IF} + f_d)nT_s+\theta) * sin(2 \pi f_{IF} * nT_s) \\ & = A*D(nT_s)C(nT_s) * 1/2 * (sin(2\pi(2f_{IF} + f_d)nT_s + \theta) - sin(2\pi f_dnT_s + \theta)) \\ \end{align} SQ=x(nTs)∗sin(2πfIF∗nTs)=A∗D(nTs)C(nTs)cos(2π(fIF+fd)nTs+θ)∗sin(2πfIF∗nTs)=A∗D(nTs)C(nTs)∗1/2∗(sin(2π(2fIF+fd)nTs+θ)−sin(2πfdnTs+θ)) -
-
低通滤波
设 K = A ∗ D ( n T s ) C ( n T s ) ∗ 1 / 2 K = A*D(nT_s)C(nT_s) * 1/2 K=A∗D(nTs)C(nTs)∗1/2
则I信号 S I = K ∗ ( c o s ( 2 π ( 2 f I F + f d ) n T s + θ ) + c o s ( 2 π f d + θ ) ) S_I = K * (cos(2\pi(2f_{IF} + f_d)nT_s+\theta) + cos(2\pi f_d + \theta)) SI=K∗(cos(2π(2fIF+fd)nTs+θ)+cos(2πfd+θ))
Q路信号 S Q = K ∗ ( s i n ( 2 π ( 2 f I F + f d ) n T s + θ ) − s i n ( 2 π f d n T s + θ ) ) S_Q = K*(sin(2\pi(2f_{IF} + f_d)nT_s + \theta) - sin(2\pi f_dnT_s + \theta)) SQ=K∗(sin(2π(2fIF+fd)nTs+θ)−sin(2πfdnTs+θ))
可以看出I/Q两路信号都有两个主要的频率分量, 2 π ( 2 f I F + f d ) 2\pi(2f_{IF} + f_d) 2π(2fIF+fd) 和 2 π f d 2\pi f_d 2πfd ,分别对应GPS导航信号与本地载波的和频与差频。如果这时滤掉
和频
保留差频
,然后找到差频的频率,就是多普勒频移。 -
降采样,减少不必要的计算
在经过低通滤波之后,信号中保留着
多普勒频移
、伪码
和数据
,其中速率最高的伪码的码率为1023kbps,根据奈奎斯特采样定理,此时数据的采样率只要高于1023*2就可以保证数据的有效性。为了尽可能的降低运算量,可以进行降采样,把采样率从原来的26MSPS降低到2046SPS。伪代码如下:
/** fs -- 原始采样率 * fs_re -- 降采样后的采样率 * N -- 原始数据的采样点数 * cnt -- 重采样后数据的下标 **/ temp = -1; for(int i = 0; i < N; i++) { cnt = 向下取整(i * fs_re / fs); if(temp == cnt - 1) { data_re(cnt) = data(i); temp = cnt; } }
-
相干积分
经过低通滤波和降采样后的数据里面主要有两种成分,伪码和多普勒频移。因为伪码的频率再在频谱上均匀分布,只有去掉伪码才能将多普勒频移找出来。伪码调制再载波信号中,是一种双极性不归零的信号(GPS导航电文也是一种双极性不归零的信号)。要去掉伪码,就需要利用双极性不归零信号的性质。
在BPSK调制中,数字信号是双极性不归零的信号。比如
1 0 1 0 1 1 0 0 需要转化为
1 -1 1 -1 1 1 -1 -1 再进行调制。
在捕获GPS信号时,需要本地生成一个标准的伪码信号,也需要是双极性不归零的信号。假设伪码就是上文中的10101100,在伪码对齐的情况下,两个相同的伪码相乘,$ 1*1=1 $ 、 $ -1 * -1 = 1 $。那么
1 -1 1 -1 1 1 -1 -1 和
1 -1 1 -1 1 1 -1 -1 对应位相乘的结果就是
1 1 1 1 1 1 1 1 伪码全部变成1就代表全部变成了直流分量,然后在做FFT或者分组求和后再做FFT就可以得到多普勒频移。再实际的操作中,需要将本伪码不断的循环右移,再相乘后做FFT,然后记录FFT的结果。将记录的结果绘制为立体图如线图所示
其中尖峰所在的x坐标和y坐标对应的坐标就是码相位和多普勒频移,图中使用256点FFT。 -
找到最高点,计算频率偏移和相位偏移
遍历记录的FFT结果,找到最高点,计算对应的频率和码相位,就是多普勒频移和码相位。
说明
- 此文描述都是我个人理解,如有错误请大佬指正
- Matlab源码和使用的数据将会上传到Gitee和Github。
Gitee: https://gitee.com/chentianya000/gps_-pmf-fft.git