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级移位寄存器（n=5），使用本原多项式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_RADIX=DEC;\n'); fprintf(fid, 'DATA_RADIX=HEX;\n'); fprintf(fid, 'CONTENT BEGIN\n'); for addr = 0:data_depth-1 fprintf(fid, ' %d : %s;\n', addr, hex_data(addr+1, :)); 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序列，其周期自相关函数理论值应为：
- 在零时延（lag=0）时，相关值为
L（归一化后为1）。 - 在其他任何非零时延（lag ≠ 0, ±L, ±2L, ...）时，相关值为
-1（归一化后为-1/L）。
我们绘制的图应该清晰地显示出在零点有一个尖锐的峰值，而在其他位置的值都非常接近-1/L（对于一个32位长的序列，L=31，-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是32（2^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或.coe（Xilinx 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变成向量，一切操作就都变得直观而简单了。