GPS信号生成：C/A码序列生成【MATLAB实现】

GPS C/A码序列生成【MATLAB实现】

在本文中，将简要介绍GPS C/A码及其生成原理，并且用MATLAB代码实现。

GPS信号与C/A码

GPS的信号主要有三类：载波（carrier）、测距码（也可以说是伪随机噪声码，Pseudo Random Noise，PRN）和导航电文（Navigation message）。GPS卫星在实现定位时就使用以上三种信号进行相应的伪距（Pseudorange）测量或载波相位（Carrier phase）测量。
本文介绍实现伪距测量的一种信号：C/A码（Coarse/Acquisition Code）是用于进行粗略测距和捕获P码的粗码，也称捕获码。周期T为1毫秒，一个周期含有码元即码长N为2^10 - 1=1023，每个码元持续的时间即码元周期t=1 ms/1023=0.977517微秒，相应的码元宽度为293.05米。C/A码是一种公开的明码，可供全球用户免费使用。但C/A码在L1载波上，所以无法精确地消除电离层延迟，测距精度一般为±（2~3）米。C/A码和P码（Precise code）同属PRN码。但是C/A码为民用，P码通常为军用，所以我们在学习和实验时基本上使用的是C/A码。
C/A码有1023个元素（chips），这里面含有512个1与511个0，他们的排列是看起来仿佛随机的，但是却是完全可决定的，因此叫做伪随机噪声序列。

GPS C/A码序列生成原理

GPS C/A码的生成通常需要经过几个步骤：如下图，窄带宽基带信号（baseband signal）经 PRN 码调制后，基带信号的频谱明显变宽。也就是说，信号从narrow bandwidth变成了wide bandwidth。基带信号经过扩频处理后，使用二进制相移键控（BPSK）调制方法对扩频信号进行调制。最后，BPSK 扩频信号与在L1频段1575.42 MHz 的载波信号一起传输。
GPS C/A码信号生成
GPS的L1频段的码率是50 bps（bit per second）。在这50 bps中，信号的频谱宽度是-50 Hz到50 Hz。另外，C/A码的码率是1.023 MHz，也就是1 ms内可以得到1023 chips的C/A码。
C/A码的生成遵循Gold code原理，这里我放了1967年的论文原文。在C/A码生成时，需要G1和G2两个线性反馈移位寄存器（Linear feedback shift register，LFSR）。每个寄存器都有10位，所以按照Gold code原理生成的序列长度是2^10-1=1023 chips。寄存器G1和G2按照异或（Exclusive or，XOR）的方式生成C/A码序列，如下图所示。一般来说，移位寄存器是一组一位存储单元（存储单元）。在时钟脉冲作用下，每个寄存器的内容将向右移动一位。输出是最后一个单元（最左边的单元或寄存器）的内容。
C/A码序列生成
需要注意的是，G1和G2分别有自己的反馈准则：

对于 G1，在每个时钟脉冲下，寄存器位置 3 和 10 的进行异或运算，并将结果移至寄存器位置 1。其他寄存器向右移动一个位置，最后寄存器位置 10 的内容成为 G1 的输出G1。对于 G2，每个时钟脉冲都会对寄存器位置 2、3、6、7、8、9、10 的进行异或运算，并将结果移至寄存器位置 1。同样，所有其他寄存器的内容也会向右移动一位。然后，选择两个selector，进行异或运算，再offset移位运算，成为G2i。最后，将两个寄存器的输出G1 和 G2i 进行异或，得到每个epoch脉冲上产生的 PRN 码位。
每个PRN对应的selector和offset的选择参考GPS Interface Specification IS-GPS-200 的表，如下图：
在这里插入图片描述
像这样重复运算1023次之后，我们就得到1 ms的C/A码序列。

MATLAB代码实现

附上生成C/A码序列的MATLAB代码，输入PRN号，得到相应的C/A码。

function code = CA_generation(prn)
% Generate the C/A code for a given PRN (1-32)
% Input:
%   prn - PRN number (1-32)
% Output:
%   code - 1023-length C/A code sequence

% Initial states for Linear feedback shift register (LFSR) of length G1 and G2
G1 = ones(1, 10); 
G2 = ones(1, 10);

% PRN taps (satellite-specific tap positions)
prn_taps_table = [
    2, 6;   % PRN1
    3, 7;   % PRN2
    4, 8;   % PRN3
    5, 9;   % PRN4
    1, 9;   % PRN5
    2, 10;  % PRN6
    1, 8;   % PRN7
    2, 9;   % PRN8
    3, 10;  % PRN9
    2, 3;   % PRN10
    3, 4;   % PRN11
    5, 6;   % PRN12
    6, 7;   % PRN13
    7, 8;   % PRN14
    8, 9;   % PRN15
    9, 10;  % PRN16
    1, 4;   % PRN17
    2, 5;   % PRN18
    3, 6;   % PRN19
    4, 7;   % PRN20
    5, 8;   % PRN21
    6, 9;   % PRN22
    1, 3;   % PRN23
    4, 6;   % PRN24
    5, 7;   % PRN25
    6, 8;   % PRN26
    7, 9;   % PRN27
    8, 10;  % PRN28
    1, 6;   % PRN29
    2, 7;   % PRN30
    3, 8;   % PRN31
    4, 9;   % PRN32
];

offset_table=[
    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;
];

% Validate PRN input
if prn < 1 || prn > 32
    error('PRN must be between 1 and 32.');
end

% Get the taps for the given PRN
selector = prn_taps_table(prn, :);
select1 = selector(1);
select2 = selector(2);

% Get the offset
offset=offset_table(prn,:);

% Generate the C/A code sequence
code = zeros(1, 1023);

for i = 1:1023
    % Generate the code bit (modulo-2 addition of selected G2 taps)
    G1_output = G1(10);
    G2_output = xor(G2(select1), G2(select2)); % Adjust G2 output by offset

    % Shift G1
    feedback1 = xor(G1(3), G1(10));
    G1 = [feedback1, G1(1:end-1)];
    
    % Shift G2
    feedback2 = xor(G2(2), xor(G2(3), xor(G2(6), xor(G2(8), xor(G2(9), G2(10))))));
    G2 = [feedback2, G2(1:end-1)];

    % Apply offset to G2
    G2_shifted = circshift(G2_output, [0, offset]);

    code(i) = xor(G1_output, G2_shifted); % Generate C/A code
end

% Convert 0 to -1 for numerical representation
code = 1 - 2 * code; % Change logic value 0 to +1 and 1 to -1
end