1. 从零开始理解m序列与Gold序列的“前世今生”如果你正在学习通信工程或者对CDMA、卫星导航这些技术背后的“密码”感到好奇那你一定绕不开两个名字m序列和Gold序列。它们不是某种神秘代码而是现代无线通信系统中用来区分不同用户、抵抗干扰的“扩频码”家族里的核心成员。今天我就带你用MATLAB这个强大的工具亲手生成它们并像侦探一样通过分析它们的“指纹”——也就是自相关和互相关特性来验证它们是否能组成一对合格的Gold序列搭档。简单来说你可以把m序列想象成一种由特定规则反馈抽头生成的、周期非常长的0/1序列。它有个很酷的特性自相关特性极好。什么叫自相关就是把这个序列和它自己错开一点位置再对比只有在完全对齐时相似度才爆表峰值很高一旦错开相似度就立刻降到很低。这个特性让接收机在嘈杂的环境里能精准地锁定信号就像在喧闹的派对上你能立刻听出自己朋友的声音一样。但m序列有个“小毛病”不同m序列之间的互相关特性可能不太好。互相关就是拿一个序列和另一个不同的序列去对比相似度。如果两个序列的互相关值总是很高那在通信系统里就麻烦了用户A的信号很容易被误认为是用户B的造成串扰。这时候Gold序列就登场了。Gold序列不是凭空造出来的它是由两个特定的、周期相同的m序列通过模二加简单理解就是异或运算生成的一整个序列家族。这个家族的成员们不仅继承了父辈m序列良好的自相关特性更重要的是家族成员之间的互相关值被严格限制在很低的水平非常适合让大量用户同时通信而互不干扰。那么问题来了随便找两个m序列就能生出优质的Gold序列家族吗当然不是这对“父母”需要满足严格的“婚配条件”。我们的实战目标就是验证给定的两个m序列比如反馈抽头为[6,1]和[6,5,2,1]的是否是一对“佳偶”。验证方法就是计算并分析它们的互相关函数。如果它们的互相关函数值只有少数几种且绝对值都很小理论上有明确的上限那它们就具备了构造Gold序列的资格。纸上谈兵不如亲手一试接下来我们就打开MATLAB一步步实现这个过程。2. 实战第一步在MATLAB中生成你的第一个m序列理论懂了手开始痒了我们这就开始敲代码。生成m序列的核心是一个叫做线性反馈移位寄存器的模型。你可以把它想象成一条有6个对应周期63即2^6-1格子的传送带每个格子放一个比特0或1。每次动作最右边的格子里的数被输出同时根据预设的“反馈抽头”规则从某几个格子里取出数做计算模二加得到的结果塞到最左边的格子里所有格子里的数整体向右移动一位。在MATLAB里我们首先需要把“反馈抽头”这个抽象概念变成代码。比如抽头[6,1]意味着参与反馈计算的是第6个寄存器和第1个寄存器通常从左往右编号为1到m。在编程时我们用一个长度为m的向量数组来表示抽头系数对应位置为1表示该寄存器参与反馈为0则不参与。对于[6,1]m6对应的向量就是[1, 0, 0, 0, 0, 1]。注意这里第1位对应最高位寄存器1第6位对应最低位寄存器6顺序别搞反了。下面是我写的一个健壮且易理解的mseq函数它比原始代码更清晰加了详细注释方便你每一步都看懂function pn_sequence mseq(feedback_taps) % 生成最大长度序列m序列 % 输入参数 % feedback_taps: 一个行向量表示反馈抽头系数。例如[1,0,0,0,0,1]对应抽头[6,1] % 输出参数 % pn_sequence: 生成的m序列逻辑0/1序列长度为 2^m - 1 m length(feedback_taps); % 确定移位寄存器的级数 N 2^m - 1; % m序列的周期最大长度 % 初始化寄存器状态通常初始化为除全零外的任何状态这里用[0 ... 0 1] registers [zeros(1, m-1), 1]; % 预分配输出序列内存提升效率 pn_sequence zeros(1, N); for i 1:N % 输出当前最右边寄存器最低位的值 pn_sequence(i) registers(end); % 计算反馈值抽头系数向量与寄存器状态向量点乘然后模2加等价于异或 % 注意MATLAB中mod(sum(a.*b), 2) 实现了模2加法 feedback_bit mod(sum(feedback_taps .* registers), 2); % 寄存器右移一位并将反馈值插入最左端最高位 registers [feedback_bit, registers(1:end-1)]; end end使用这个函数超级简单。在MATLAB命令窗口或你的脚本里直接调用即可生成序列% 定义反馈抽头系数向量 tap1 [1, 0, 0, 0, 0, 1]; % 对应抽头 [6, 1] tap2 [1, 1, 0, 0, 1, 1]; % 对应抽头 [6, 5, 2, 1] tap3 [1, 0, 0, 1, 1, 1]; % 对应抽头 [6, 5, 4, 1] 注意原始文章此处描述可能有误我们以代码为准 % 生成m序列 pn_seq1 mseq(tap1); % 序列1 pn_seq2 mseq(tap2); % 序列2 pn_seq3 mseq(tap3); % 序列3 % 查看前10个码片 disp(序列1前10个码片); disp(pn_seq1(1:10));生成的pn_seq1等是一个由0和1组成的、长度为63的序列。但等一下在计算相关性时我们通常希望序列的值是1和-1因为这样计算相关函数更方便对应BPSK调制。所以我们还需要一个简单的转换% 将0/1序列转换为1/-1序列 y1 2 * pn_seq1 - 1; % 0 - -1, 1 - 1 y2 2 * pn_seq2 - 1; y3 2 * pn_seq3 - 1;好了现在y1,y2,y3就是我们后续进行相关性分析的“原材料”。你可以试着画一下其中一个序列的片段看看它看起来是不是像一串随机的正负一脉冲这就是伪随机序列的特点。3. 核心武器编写相关性计算函数并深入解读序列有了接下来就是检验它们特性的关键步骤——计算相关函数。这是整个实验的灵魂。我强烈建议你不要直接调用MATLAB内置的xcorr函数了事自己动手写一个才能真正理解周期相关是怎么算的这对理解通信原理至关重要。自相关函数描述的是一个序列和它自身经过不同时延τ后的相似程度。对于周期为N的序列我们计算所有可能循环移位后的内积。互相关函数则是计算两个不同序列在不同相对时延下的相似程度。理想情况下我们希望自相关函数在零时延处有一个尖锐的高峰方便同步而在其他时延处值尽可能低抗多径干扰希望互相关函数在所有时延处的绝对值都尽可能小减少用户间干扰。下面是我优化后的pncorr函数。它同时计算了循环左移和循环右移并将结果拼接成一个完整的、中心为零时延的相关函数序列长度是2N-1非常便于绘图观察。function corr_result pncorr(seq_a, seq_b, N) % 计算两个周期序列的周期互相关函数若输入相同序列则为自相关 % 输入参数 % seq_a: 第一个序列1/-1格式 % seq_b: 第二个序列1/-1格式必须与seq_a等长 % N: 序列的周期长度 % 输出参数 % corr_result: 长度为2*N-1的相关函数序列中心点第N个元素对应零时延 % 初始化存储左移和右移相关结果的数组 left_shift_corr zeros(1, N-1); % 存储τ1,2,...,N-1seq_a左移 right_shift_corr zeros(1, N); % 存储τ0,-1,...,-(N-1)seq_a右移 % 计算零时延互相关τ0 right_shift_corr(1) sum(seq_a .* seq_b); % 计算seq_a循环左移τ位τ为正时与seq_b的相关 for tau 1:N-1 shifted_a [seq_a(tau1:end), seq_a(1:tau)]; % 循环左移tau位 left_shift_corr(tau) sum(shifted_a .* seq_b); end % 计算seq_a循环右移τ位τ为负这里用绝对值时与seq_b的相关 for tau 1:N-1 shifted_a [seq_a(N-tau1:end), seq_a(1:N-tau)]; % 循环右移tau位 right_shift_corr(tau1) sum(shifted_a .* seq_b); end % 组合结果左移部分τ正 右移部分τ零和负 % 注意顺序left_shift_corr是τ N-1, N-2, ..., 1 % 我们需要将其反转使其对应τ 1, 2, ..., N-1 corr_result [fliplr(left_shift_corr), right_shift_corr]; end让我解释一下这个函数的精妙之处。它输出的corr_result是一个长度为125对于N63的向量。中间第63个点索引63对应零时延τ0。向左走索引减小对应τ为正seq_a相对左移向右走索引增大对应τ为负seq_a相对右移。这样画图时横坐标从-62到62图像关于零时延对称对于自相关函数而言非常直观。为什么要自己写而不直接用xcorr因为xcorr默认计算的是线性相关对于周期序列需要指定biased或unbiased选项并且结果的范围和归一化方式需要小心处理。自己写的函数概念更清晰更能体现周期相关的本质——循环移位后点积。你可以用这个函数计算自相关看看Rxx pncorr(y1, y1, 63);你会发现零时延的值是63序列长度因为每个1或-1自乘都是163个1相加就是63这就是峰值。4. 可视化对决绘制并解读相关函数图谱计算出的数据是冰冷的图形才是直观的判决书。我们把四个关键的相关函数画在一起进行对比分析。下面的脚本将生成一个2x2的子图让你一眼看清所有特性。clc; clear; close all; % 清空环境 N 63; % 生成序列使用前面定义的函数和抽头 tap1 [1,0,0,0,0,1]; tap2 [1,1,0,0,1,1]; tap3 [1,0,0,1,1,1]; pn1 mseq(tap1); pn2 mseq(tap2); pn3 mseq(tap3); y1 2*pn1 - 1; y2 2*pn2 - 1; y3 2*pn3 - 1; % 计算相关函数 R11 pncorr(y1, y1, N); % 序列1自相关 R22 pncorr(y2, y2, N); % 序列2自相关 R12 pncorr(y1, y2, N); % 序列1与序列2互相关 R23 pncorr(y2, y3, N); % 序列2与序列3互相关 % 设置时延轴 tau -(N-1):(N-1); % 从-62到62 % 开始绘图 figure(Position, [100, 100, 1000, 700]); % 设置图形窗口大小 % 子图1序列1抽头[6,1]的自相关函数 subplot(2,2,1); stem(tau, R11, b, Marker, none, LineWidth, 1.2); xlabel(时延 \tau); ylabel(R_{11}(\tau)); title((a) m序列 [6,1] 的自相关函数); grid on; axis tight; hold on; plot([0,0], ylim, k--, LineWidth, 0.5); % 在零时延处画虚线 % 突出峰值 plot(0, R11(N), ro, MarkerSize, 8, MarkerFaceColor, r); legend(自相关, 零时延, 峰值点, Location, best); % 子图2序列2抽头[6,5,2,1]的自相关函数 subplot(2,2,2); stem(tau, R22, r, Marker, none, LineWidth, 1.2); xlabel(时延 \tau); ylabel(R_{22}(\tau)); title((b) m序列 [6,5,2,1] 的自相关函数); grid on; axis tight; hold on; plot([0,0], ylim, k--, LineWidth, 0.5); plot(0, R22(N), ro, MarkerSize, 8, MarkerFaceColor, r); % 子图3序列1与序列2的互相关函数 subplot(2,2,3); stem(tau, R12, g, Marker, none, LineWidth, 1.2); xlabel(时延 \tau); ylabel(R_{12}(\tau)); title((c) [6,1] 与 [6,5,2,1] 的互相关函数); grid on; axis([-N1, N-1, -20, 20]); % 固定y轴范围便于比较 hold on; plot(xlim, [0,0], k-, LineWidth, 0.5); % 画零线 % 标注最大值和最小值 [max_val, max_idx] max(R12); [min_val, min_idx] min(R12); plot(tau(max_idx), max_val, ^, Color, [0.8,0,0], MarkerSize, 10, MarkerFaceColor, [0.8,0,0]); plot(tau(min_idx), min_val, v, Color, [0,0.5,0], MarkerSize, 10, MarkerFaceColor, [0,0.5,0]); legend(互相关, 零线, sprintf(最大值: %.0f, max_val), sprintf(最小值: %.0f, min_val), Location, best); % 子图4序列2与序列3的互相关函数 subplot(2,2,4); stem(tau, R23, m, Marker, none, LineWidth, 1.2); xlabel(时延 \tau); ylabel(R_{23}(\tau)); title((d) [6,5,2,1] 与 [6,5,4,1] 的互相关函数); grid on; axis([-N1, N-1, -20, 20]); hold on; plot(xlim, [0,0], k-, LineWidth, 0.5); [max_val2, max_idx2] max(R23); [min_val2, min_idx2] min(R23); plot(tau(max_idx2), max_val2, ^, Color, [0.8,0,0], MarkerSize, 10, MarkerFaceColor, [0.8,0,0]); plot(tau(min_idx2), min_val2, v, Color, [0,0.5,0], MarkerSize, 10, MarkerFaceColor, [0,0.5,0]); legend(互相关, 零线, sprintf(最大值: %.0f, max_val2), sprintf(最小值: %.0f, min_val2), Location, best); % 整体标题 sgtitle(m序列与Gold序列候选对相关性分析对比 (N63), FontSize, 14, FontWeight, bold);运行这段代码你会得到一张信息量丰富的图。现在我们来当一回“图谱医生”诊断每一幅图图(a)和图(b)这是两个m序列的自相关函数。你应该能看到一个非常典型的“冲击函数”形状在时延τ0处有一个高达63的尖峰而在所有其他时延位置τ±1, ±2, ...相关值都是一个很小的负数-1。这个“主瓣尖锐旁瓣平坦且低”的特性正是m序列作为优秀同步码和扰码的资本。接收机通过滑动相关寻找这个峰值就能精确确定信号的起始时刻。图(c)这是重头戏候选Gold序列对[6,1]和[6,5,2,1]的互相关函数。仔细观察它的纵坐标值。你会发现它不像自相关那样只有一个高峰而是在多个时延上起伏。但关键点在于这些起伏的值被限制在有限的几个电平上。通过unique(R12)命令你可以看到它可能只取-17, -1, 15这几个值具体值可能因初始相位略有不同。更重要的是其最大绝对值本例中为17远小于自相关的峰值63。根据Gold序列理论对于周期N63的序列优选对互相关的最大绝对值上限约为12^((m2)/2) ≈ 12^4 17。我们的计算结果符合这个理论值这是一个强烈的信号表明这对m序列是“优选对”。图(d)作为对比我们计算了序列[6,5,2,1]与另一个序列[6,5,4,1]的互相关。你会发现它的波形看起来更“杂乱”最大值和最小值的绝对值很可能更大可能超过20甚至更高。这直观地说明了不是任意两个m序列都能组成低互相关的“优选对”。选择错误的配对会导致在CDMA系统中产生严重的多址干扰。5. 工程意义与深入探索从验证到应用通过上面的计算和绘图我们已经从数据上验证了抽头为[6,1]和[6,5,2,1]的m序列是一对合格的Gold序列生成器“优选对”。但这在工程上到底意味着什么让我结合CDMA系统的例子再深入聊聊。在CDMA手机网络或GPS中每个用户或每颗卫星都被分配一个独特的扩频码。Gold序列家族就是从一对优选对m序列派生出来的庞大码组。假设我们有了优选对G1和G2那么通过将G2相对于G1进行不同的循环移位然后模二加就能生成一系列新的序列Gold(i) G1 ⊕ shift(G2, i)其中i取0到N-1。这样我们就能得到N2个序列包括G1、G2和N个Gold序列。这些序列家族内的互相关特性都继承了其“父母”优选对的优良基因即互相关值被限制在很低的水平。我们可以用MATLAB简单演示一下生成一个Gold序列家族并抽查其中几个序列间的互相关% 假设y1, y2已经是1/-1格式的优选对序列 G1 y1; % 对应[6,1] G2 y2; % 对应[6,5,2,1] % 生成Gold序列家族这里生成前5个作为示例 GoldFamily zeros(5, N); for i 0:4 shifted_G2 circshift(G2, i); % 循环移位G2 GoldFamily(i1, :) G1 .* shifted_G2; % 注意对于1/-1序列模二加等价于逐元素相乘 % 因为 (1) XOR (1) 0 - 映射为 -1? 这里需要小心。 % 更严谨的做法是先用回0/1格式做XOR再转成1/-1。 % 以下为更清晰的生成方法 pn_G1 (G11)/2; % 1/-1 转 0/1 pn_G2_shift (circshift(G2, i)1)/2; gold_bits xor(pn_G1, pn_G2_shift); % 模二加 GoldFamily(i1, :) 2 * gold_bits - 1; % 转回1/-1 end % 计算家族内序列0和序列4的互相关 R_gold pncorr(GoldFamily(1,:), GoldFamily(5,:), N); figure; stem(tau, R_gold); title(Gold家族内两个序列的互相关函数示例); xlabel(时延 \tau); ylabel(R(\tau)); grid on; axis([-N1, N-1, -20, 20]);你会发现这样生成的序列间的互相关其最大值仍然被限制在较低的水平。这就是Gold序列在工程上如此强大的原因用少量的硬件两个线性反馈移位寄存器就能产生大量互相关性良好的伪随机码极大地节省了资源。在实际项目中选择哪一对m序列作为优选对是有表格可查的称为“优选对表”它基于本原多项式的理论。我们的实验验证了表格中的一对。通过这个实战你不仅学会了用MATLAB生成和分析序列更重要的是你理解了背后“为什么”——为什么Gold序列能用于CDMA为什么互相关特性如此关键。下次当你听到“扩频通信”、“码分多址”这些词时你脑子里浮现的就不再是抽象概念而是一幅幅由MATLAB绘制的、跳动着的相关函数图以及那对精心挑选的、共同孕育出一个低干扰码家族的“优选对”m序列。