Matlab中M序列循环移位实现与自相关验证 1. 项目概述与背景最近在折腾一个扩频通信相关的项目核心需求是生成一组性能优良的扩频码。对于入门或者对实时性要求不那么极致的场景M序列最大长度序列是个非常经典且实用的选择。它具备良好的自相关和互相关特性特别适合用来做码分多址或者抗干扰。我的思路是先生成一个本原M序列然后通过循环移位的方式派生出多组不同的码序列这样就能用一套简单的硬件或算法支撑起多个用户的扩频需求。在Matlab里实现这个想法时我发现“移位”这个看似简单的操作其实藏着不少细节。Matlab提供了circshift和bitshift这两个函数但它们的适用场景和效果完全不同。直接用bitshift做循环移位会掉数据而用circshift处理二进制向量又涉及到数据类型的转换。我花了一些时间摸索最终整理出了一套从生成、移位到性能验证自相关计算的完整流程。这篇文章我就把这些踩过的坑和验证过的方法详细地分享出来特别是如何高效地实现比特级的循环移位以及如何直观地验证生成的M序列是否具备我们期望的尖锐自相关特性。2. M序列基础与循环移位原理2.1 什么是M序列及其在扩频通信中的价值M序列全称最大长度线性反馈移位寄存器序列是由一个n级线性反馈移位寄存器在特定反馈逻辑下产生的周期为2^n - 1的伪随机二进制序列。它之所以在扩频通信里备受青睐核心在于它的几个关键特性首先是尖锐的自相关特性在零时延处自相关值达到最大序列长度而在其他时延处自相关值迅速跌落到一个很小的负值-1这个特性对于接收端的同步捕获至关重要其次是良好的互相关特性虽然同族M序列之间的互相关值不是绝对理想但通过优选对或者Gold序列等方式可以改善对于由同一个M序列循环移位产生的不同序列它们在完整周期内的互相关性能也是有理论保障的。在工程上我们常常用一个本原多项式来定义M序列的生成逻辑。例如对于一个5级移位寄存器n5使用本原多项式x^5 x^2 1可以产生一个周期为31的M序列。在Matlab中我们可以用prbs函数或自己编写LFSR逻辑来生成它。得到这个“母序列”后循环移位的价值就体现出来了它的每一个循环移位版本本身也是一个周期相同的M序列并且与原始序列及其他移位序列之间在完整周期内具有固定的、可计算的互相关值。这意味着我只需要存储或生成一个序列就能通过简单的移位操作得到多个可用的扩频码极大地节省了存储资源或计算开销在FPGA或嵌入式系统中尤其有用。2.2 Matlab中的移位操作circshift与bitshift深度解析Matlab提供了两种主流的移位函数但它们的操作对象和逻辑有本质区别用错了地方结果会南辕北辙。circshift面向数组的循环移位这个函数操作的对象是数组向量、矩阵的元素。它的逻辑是将数组的边界视为首尾相连的环移出边界的元素会从另一侧补回来。正如我例子中所示x [1,2,3;4,5,6;7,8,9]; y circshift(x, [1, -1]);这里[1, -1]表示整个矩阵向下循环移动1行同时向左循环移动1列。元素7移动到了(1,2)的位置元素3移动到了(2,3)的位置。它不关心元素本身代表什么是数字、字符还是别的只负责改变元素在容器中的位置。当我们把一个二进制序列表示成一个由字符0和1组成的向量时circshift就可以对这个向量进行循环移位。bitshift面向整数的比特移位这个函数操作的对象是整数的二进制位。它的逻辑是经典的逻辑移位或算术移位默认是逻辑移位移出边界的比特会被直接丢弃空出的位补零。x 131; % 二进制 10000011 y bitshift(x, 2); % 左移2位结果为 1000001100但Matlab默认用双精度存储显示为524 z bitshift(x, 2, 8); % 左移2位并只保留低8位。10000011左移2位是1000001100取低8位是00001100即12。关键问题来了bitshift不支持循环移位。左移时高位移出就丢了右移时低位移出也丢了。这对于需要保持序列完整性的M序列生成来说是致命的。因此如果我们直接对一个代表M序列的十进制整数使用bitshift得到的将不再是同一个M序列的循环移位而是被截断或补零后的另一个数。注意bitshift的第三个参数n非常有用它指定了将输入数据视为一个n位无符号整数进行移位和截断。这在处理固定位宽数据如FPGA中的寄存器时能精确模拟硬件行为。那么如何实现比特级的循环移位呢核心思路是将整数转换为一个显式的比特数组向量然后对这个数组进行循环移位最后再转换回整数。这样我们就将“比特循环移位”这个bitshift做不到的功能转化成了circshift可以轻松处理的数组操作。3. 循环移位M序列的生成与实现3.1 从十进制到二进制向量的转换策略既然要对比特进行操作第一步就是把代表M序列的十进制数转换成一个由0和1组成的向量。这里有几个方法各有优劣方法一使用dec2bin和字符数组这是我最开始采用的方法也是比较直观的一种。x_decimal 3139384450; % 假设这是一个32位M序列对应的十进制数 x_bin_str dec2bin(x_decimal, 32); % 转换为32位宽的二进制字符串例如 10111100101100101001100100000010 x_bit_vector x_bin_str - 0; % 巧妙利用ASCII码差将字符数组转换为数值向量dec2bin(x, n)会生成一个长度为n的字符串非常适合固定位宽的场景。后续的- 0是一个经典技巧因为字符0和1的ASCII码是连续的相减后直接得到数值0和1。方法二使用bitget函数更推荐这是更专业、更直接的方法尤其适合处理固定位宽的硬件数据。x_decimal 3139384450; % bitget 从最低位LSB开始获取指定位的值。为了得到从最高位MSB到最低位的向量需要反转索引。 x_bit_vector bitget(x_decimal, 32:-1:1); % 获取第32位到第1位的值组成向量bitget直接操作整数的比特位无需经过字符串转换效率更高语义也更清晰。32:-1:1这个索引顺序确保了生成的向量x_bit_vector是[MSB, ..., LSB]的顺序这与我们通常书写二进制数的习惯一致。3.2 实现高效的比特循环移位得到了比特向量x_bit_vector后实现循环移位就非常简单了直接使用circshift即可。% 假设 x_bit_vector 是一个行向量 [1, 0, 1, 1, ...] shift_amount 5; % 循环左移5位 shifted_bit_vector circshift(x_bit_vector, [0, -shift_amount]);这里[0, -shift_amount]表示在行方向第一维不移位在列方向第二维向左移动shift_amount位负数表示左移。circshift会自动将移出左边的比特从右边补回来完美实现了循环移位。如果需要循环右移则将参数改为[0, shift_amount]正数表示右移。3.3 生成完整移位序列集并导出在实际项目中我通常需要生成一整套例如32个循环移位后的序列并将它们以某种格式导出供其他模块如FPGA使用。下面是一个完整的示例流程%% 1. 定义参数与初始序列 sequence_length 32; % M序列位宽 m_sequence_decimal 3139384450; % 种子序列需确保是有效的2^n-1周期序列此处为示例值 % 生成初始比特向量 (MSB - LSB) initial_bits bitget(m_sequence_decimal, sequence_length:-1:1); %% 2. 生成所有循环移位序列 num_sequences sequence_length; % 生成与位宽相同数量的序列 all_sequences_bits zeros(num_sequences, sequence_length); % 预分配矩阵每行一个序列 all_sequences_bits(1, :) initial_bits; % 第一行是原始序列 for i 2:num_sequences % 每次在前一个序列的基础上左移1位 all_sequences_bits(i, :) circshift(all_sequences_bits(i-1, :), [0, -1]); end %% 3. 转换为十进制数可选便于查看或某些接口使用 all_sequences_decimal zeros(num_sequences, 1); for i 1:num_sequences % 将比特向量转换为十进制数 % 方法将二进制向量转换为字符串再用 bin2dec 转换。注意 bin2dec 接收字符向量。 bin_str sprintf(%d, all_sequences_bits(i, :)); all_sequences_decimal(i) bin2dec(bin_str); end %% 4. 导出数据例如导出为.mif文件供FPGA使用 % 假设我们需要将32个32位的十进制数写入.mif文件 % .mif (Memory Initialization File) 是Altera/Intel FPGA常用的一种内存初始化格式。 data_width 32; data_depth num_sequences; hex_data dec2hex(all_sequences_decimal, ceil(data_width/4)); % 转换为16进制字符串 fid fopen(m_sequences.mif, w); fprintf(fid, -- M序列循环移位组深度%d位宽%d\n, data_depth, data_width); fprintf(fid, WIDTH%d;\n, data_width); fprintf(fid, DEPTH%d;\n, data_depth); fprintf(fid, ADDRESS_RADIXDEC;\n); fprintf(fid, DATA_RADIXHEX;\n); fprintf(fid, CONTENT BEGIN\n); for addr 0:data_depth-1 fprintf(fid, %d : %s;\n, addr, hex_data(addr1, :)); end fprintf(fid, END;\n); fclose(fid); disp(M序列.mif文件生成完毕。);实操心得在将比特向量转换回十进制时bin2dec函数接收的是字符向量。使用sprintf(%d, bit_vector)可以快速将数值向量[1,0,1]转换成字符向量101比用num2str后再去除空格要更简洁高效。另外在生成用于FPGA的.mif文件时使用十六进制HEX格式比二进制BIN更紧凑可读性也更好。4. M序列自相关特性的计算与验证生成了序列我们必须要验证它是否具备M序列应有的尖锐自相关特性。这是衡量序列性能、确认生成过程是否正确无误的关键一步。4.1 使用xcorr函数计算自相关Matlab的信号处理工具箱提供了强大的xcorr函数用于计算互相关和自相关。对于二进制序列我们可以直接使用。%% 计算单个M序列的自相关函数 % 使用之前生成的 initial_bits (它是一个由0和1组成的向量) % 注意xcorr 默认计算的是‘非归一化’的原始相关值。 auto_corr_result xcorr(initial_bits); % 为了更直观我们通常计算归一化的自相关函数 % 归一化自相关 xcorr(sequence) / (序列长度) normalized_auto_corr auto_corr_result / length(initial_bits); %% 绘制自相关函数图 lag -(length(initial_bits)-1):(length(initial_bits)-1); % 时延滞后轴 figure; stem(lag, normalized_auto_corr, filled, MarkerSize, 4); xlabel(Lag (chip)); ylabel(Normalized Correlation); title(自相关函数 of M-sequence); grid on;对于一个理想的、周期为L 2^n - 1的M序列其周期自相关函数理论值应为在零时延lag0时相关值为L归一化后为1。在其他任何非零时延lag ≠ 0, ±L, ±2L, ...时相关值为-1归一化后为-1/L。我们绘制的图应该清晰地显示出在零点有一个尖锐的峰值而在其他位置的值都非常接近-1/L对于一个32位长的序列L31-1/L ≈ -0.0323形成一个低矮的“基底噪声”。4.2 验证循环移位序列的自相关与互相关生成了多组循环移位序列后我们还需要验证两件事1每个移位序列自身的自相关特性是否保持2不同移位序列之间的互相关特性如何。%% 验证所有循环移位序列的自相关特性是否一致 figure; hold on; for i 1:min(5, num_sequences) % 为了图面清晰只画前5个 corr_val xcorr(all_sequences_bits(i, :)); normalized_corr corr_val / sequence_length; plot(lag, normalized_corr, DisplayName, [Shift , num2str(i-1)]); end hold off; xlabel(Lag); ylabel(Normalized Correlation); title(自相关函数 of Different Cyclic Shifts); legend show; grid on; % 观察所有曲线是否基本重合仅在峰值位置发生平移。 %% 计算并观察互相关示例计算第一个序列和第三个序列的互相关 seq1 all_sequences_bits(1, :); seq3 all_sequences_bits(3, :); cross_corr_result xcorr(seq1, seq3); normalized_cross_corr cross_corr_result / sequence_length; figure; stem(lag, normalized_cross_corr, filled, MarkerSize, 4); xlabel(Lag); ylabel(Normalized Cross-Correlation); title(互相关函数 between Seq(0) and Seq(2)); grid on; ylim([-0.2, 0.2]); % 限制y轴范围更清晰地观察互相关值对于由同一个M序列循环移位产生的两个不同序列它们的互相关函数值在整个周期内也是三值特性通常为几个特定的值但其绝对值远小于自相关的峰值。在图中我们应该看到互相关值在一个较小的范围内波动没有像自相关那样突出的单一峰值。这是码分多址系统能够区分不同用户的基础。4.3 自相关计算中的细节与陷阱序列长度与周期我们上面计算的是非周期的自相关/互相关即只对给定的一段有限长序列进行计算。而M序列的优良特性是定义在其完整周期上的。如果输入xcorr的序列长度不是一个完整周期计算结果会偏离理论值。因此最好确保用于计算的initial_bits长度等于2^n - 1。如果我们的sequence_length是322^5那它可能不是一个标准的M序列长度31才是需要检查你的生成多项式。二进制表示xcorr函数将输入序列视为实数序列。对于二进制序列[0, 1]其自相关特性与将其映射为[-1, 1]的双极性序列[-1, 1]是不同的。在通信理论中我们通常使用双极性表示因为这样零时延的自相关值等于序列长度L非零时延为-1归一化后就是1和-1/L对比更鲜明。你可以通过转换来计算bipolar_sequence 2 * initial_bits - 1; % 将0/1映射为-1/1 auto_corr_bipolar xcorr(bipolar_sequence) / length(bipolar_sequence);使用双极性序列计算得到的图形在零时延的峰值会更“干净”旁瓣非零时延的相关值会更平整地处于-1/L的水平线上。xcorr的缩放选项xcorr(a, ‘normalized’)或xcorr(a, ‘coeff’)可以直接得到归一化的相关结果其零时延的系数为1。这对于快速查看形状很方便。但要注意它的归一化方式是除以sqrt(sum(a.*a) * sum(b.*b))的估计值对于计算互相关时理解具体数值含义有时自己手动归一化更清晰。5. 常见问题、调试技巧与工程化考量5.1 为什么我生成的序列自相关特性不理想这是调试阶段最常见的问题。可能的原因和排查步骤如下种子值错误你使用的十进制数m_sequence_decimal可能根本不是一个有效的M序列。M序列是由线性反馈移位寄存器产生的不是任意一个32位数都行。解决方案使用标准的本原多项式来生成序列。例如用通信工具箱的prbs函数或者自己编写一个LFSR生成函数。% 示例使用通信工具箱生成周期为31的M序列 if exist(prbs, file) p prbs(5, 1); % 5阶初始种子为1非全零 % prbs输出可能是双极性或单极性查看文档确认。假设输出0/1。 m_seq_bits p(1:31); % 取一个完整周期 % 然后可以将这个0/1序列转换成一个十进制数作为你的种子 seed_decimal bin2dec(sprintf(%d, m_seq_bits)); else % 手动实现一个5级LFSR反馈多项式 x^5 x^2 1 register [1 0 0 0 0]; % 初始状态非全零 m_seq_bits zeros(1, 31); for i 1:31 m_seq_bits(i) register(5); % 输出最高位 feedback xor(register(5), register(2)); % 根据多项式计算反馈 register [feedback, register(1:4)]; % 右移反馈输入最低位 end end序列长度非周期你生成的比特向量长度可能不是2^n - 1。比如你生成了32位但标准的5阶M序列周期是31位。多出来的一位破坏了循环结构。解决方案检查你的序列长度是否符合2^n - 1。在循环移位时确保移位位数是对序列长度取模运算circshift(seq, [0, -mod(shift, length(seq))])。数据类型混淆在循环移位和转换过程中数据类型可能在数值double、逻辑logical和字符char之间意外转换导致计算错误。解决方案在关键步骤后使用class()函数检查变量类型并使用whos查看变量维度。确保进行xcorr计算前输入是一个数值向量double或logical。5.2 如何将Matlab生成的序列用于FPGA或嵌入式系统格式转换如前面所示生成.mif或.coeXilinx FPGA使用文件是最直接的方式。确保导出的数据位宽与FPGA中ROM或移位寄存器的位宽一致。资源优化在FPGA中有时不需要存储所有移位序列。更常见的做法是只存储一个本原M序列在硬件上用一个线性反馈移位寄存器实时生成并通过使能信号控制其循环移位状态。这比存储一个大ROM更节省资源。同步与初始化无论在Matlab中仿真还是在硬件中实现LFSR的初始状态种子不能是全零否则寄存器会一直输出零。必须确保有一个非零的初始化过程。测试向量生成可以将Matlab计算出的自相关函数理论值特别是峰值和旁瓣值也导出作为FPGA仿真时的参考标准用于验证硬件相关器计算的结果是否正确。5.3 性能与扩展思考计算效率当需要生成大量、极长周期的M序列并进行相关运算时循环移位和xcorr可能会成为性能瓶颈。对于超长序列可以考虑使用频域方法计算相关利用FFT和卷积定理即ifft( fft(seq1) .* conj(fft(seq2)) )这在序列长度很大时速度远快于时域的xcorr。Gold序列与优选对如果项目对互相关性要求极高单一的M序列族可能不够。可以考虑使用两个特定本原多项式生成的M序列进行模二加产生Gold序列。Gold序列族拥有更多可用的序列且最大互相关值有理论上限。在Matlab中你可以生成两个M序列然后通过按位异或xor来生成Gold序列。多进制扩频码M序列是二进制的。在某些应用场景下可能需要多进制如四进制、复数的扩频码以获得更好的性能。其生成和循环移位的原理类似但操作对象从比特向量变成了符号向量circshift依然适用。整个流程走下来从理解M序列的原理到在Matlab中克服比特移位的障碍再到验证其核心的相关特性最后考虑如何落地到硬件是一个完整的通信算法链路仿真与验证过程。最关键的是不要被bitshift的名字迷惑对于循环移位circshift配合清晰的比特向量表示才是更通用、更可靠的工具。把序列用bitget变成向量一切操作就都变得直观而简单了。