M序列的特征多项式与生成电路[原创cnblogs.com/helesheng]

M序列生成器是基于线性反馈移位寄存器（LFSR）的伪随机序列产生电路，通过移位寄存器与异或反馈结构，生成周期最长、统计特性接近白噪声的二进制序列，广泛用于通信扩频、数据加密、硬件测试中的伪随机信号源。

一、问题和背景

第十届全国大学生电子设计竞赛F题还要求实现过特征多项式为\(f(x)=1+x^2+x^3+x^4+x^8\)的M序列，最近想设计一个考题，要求由已知的M序列反推生成这个M序列的特征多项式。最初的想法是让学生使用经典的Berlekamp-Massey 算法（简称B-M 算法，也译作"伯利坎普-梅西算法"）来实现。我自己做了一下，发现怎么也调试不对，用AI写的MATLAB程序也总做不对。仔细检查后才发现，问题发生在反推特征多项式之前的步骤——特征多项式的实现上。这一步就错，后续的B-M算法自然怎么也无法得到正确的结果。

先做自我检讨，我早年参与编写的一本书《基于FPGA技术的工程应用与实践》（9787030574299，科学出版社2019）中我负责的最后一章中的关于M序列生成电路就是错误的，当时我把特征多项式中的常数项1也理解成了需要参与反馈的项（关于这一点的数学原理，请参见下文）。但是由于电子设计竞赛评分标准中并未要求验证生成序列与所要求的特征多项式是否匹配，导致当时未发现这个错误。另外，网络上https://blog.51cto.com/u_14152/8173741的理解也是错误的。这篇网文的错误在于把特征多项式中次数和移位寄存器编号的对应关系搞反了。特征多项式中出现次数为i的项，意味着i个时钟节拍之前送入移位寄存的值将参与异或计算产生反馈值。

（以下内容欢迎转载，请注明出处： https://www.cnblogs.com/helesheng）

二、正确的对应关系

直接上第十届电赛那个题作为实例：

注意事项：

1、特征多项式中x的次数就是延迟的时钟节拍数（上图没有画出寄存器的同步时钟，不是没有，而是为了突出重点省略了），n次特征多项式，需要n个寄存器生成（虽然n次本原多项式有n+1个系数）。特征多项式中的“+”运算，表示的是\(GF(2^N)\)域中的第一种群运算“异或”，并非加法。

2、n次特征多项式的n次系数和0次系数一定为1。n次系数表示M序列生成器共保存了n次之前的输入，寄存器一定有n个；0次系数表示的是移位寄存的输入（如上图），因此系数也必须为1（我觉得，把代表输入的常数项放到等式的左边更确切，——表明反馈的输入是怎样计算出来的，但实际上异或运算的逆运算也是异或运算，这个代表输入项的1就是从等式的左边搬过来的，只不过从从左边搬过来后运算变为逆运算还是异或“+”而已）。

3、“本原多项式”的本质是使得M序列具有\(2^8-1=255\)个状态（也就是生成的伪随机序列循环周期最长）的多项式，从循环群理论解释，本原多项式就是和多项式（\(f(x)=1+x^1+x^2+x^3+x^4+x^5+x^6+x^7+x^8\)）“互质”的多项式（因为只有这样它才是循环群的生成元）。但不是循环群生成元的多项式也可以生成伪随机序列，只是它们生成的都是循环的子群，状态数小于255，且可以整除。这样的多项式不能叫做“本原多项式”，但可以更一般的称为所生成序列的“特征多项式”。

4、用硬件描述这个M序列生成器时，按一般习惯，用Verilog HDL如下：

// 8 级移位寄存器，上电初值非全零（确保自启动）  
reg [7:0] shift_reg = 8'b0000_0001;
// 输出最高位
assign seq_out = shift_reg[7];
assign feedback = shift_reg[1] ^ shift_reg[2] ^ shift_reg[3] ^ shift_reg[7];
always @(posedge ext_clk) begin
        shift_reg[7:0] <= {shift_reg[6:0], feedback};
end

其中shift_reg后面的编号，由于移位寄存器存在一个节拍延迟。因此如shift_reg[0]实际上是1个节拍延迟的输入，写程序时，不能直接用特征多项式中x的次数。与上图中相同的特征多项式\(f(x)=1+x^2+x^3+x^4+x^8\)，应该描述为feedback = shift_reg[1] ^ shift_reg[2] ^ shift_reg[3] ^ shift_reg[7];

三、Berlekamp-Massey算法的实现

最后附上Berlekamp-Massey算法的MATLAB实现，输入是测得的M序列输出，输出是可以生成该输出的最小特征多项式。一般规律是，猜测N阶（次）特征多项式，需要提供2N位的M序列输出。

function [C, L, sufficient] = BMA(s)
% BMA  用 Berlekamp-Massey 算法求二进制序列的极小特征多项式
%
% 输入:  s —— 长度为 N 的行向量，元素取 0/1（GF(2)）
% 输出:  C —— 长度为 L+1 的行向量 [1 c1 c2 ... cL]，即 C(x) = 1 + c1*x + ... + cL*x^L
%        L —— 线性复杂度（LFSR 级数）
%        sufficient —— 数据充分性标志:  true  = 2L <= N，结果唯一确定
%                                       false = 2L >  N，结果欠定（数据不足）
%
% 递推约定: s(n) = c1*s(n-1) + ... + cL*s(n-L)  (mod 2, n > L)
%
% 示例:
%   s = [1 0 0 0 1 0 1 0 1 1 0 0 1 1 0 0];
%   [C, L, ok] = BMA(s)   % -> C = [1 0 1 1 1 0 0 0 1], L = 8, ok = true

N = length(s);
s = s(:).';                              % 统一为行向量
assert(all(s == 0 | s == 1), '输入序列必须只含 0 和 1');

C = zeros(1, N+1);  C(1) = 1;            % 当前连接多项式 C(x) = 1
B = zeros(1, N+1);  B(1) = 1;            % 失配存档多项式   B(x) = 1
L = 0;                                   % 当前线性复杂度
b = 1;                                   % 存档多项式首次生效时的失配值
m = 1;                                   % 距上次失配的步数（修正项移位量）

for n = 1:N
    % ---- ① 失配值 d = s(n) + sum c_i * s(n-i)  (mod 2) ----
    d = s(n);
    for i = 1:L
        if C(i+1)
            d = bitxor(d, s(n-i));
        end
    end

    if d == 0
        % ---- ② 无失配 ----
        m = m + 1;
    elseif 2*L <= n-1
        % ---- ③ 失配且允许长度扩展: C <- C + x^m*B ----
        T = C;
        C = xor(C, shiftPoly(B, m, N+1));
        L = n - L;
        B = T;
        b = d;
        m = 1;
    else
        % ---- ④ 失配但数据不足，仅微调: C <- C + (d/b)*x^m*B ----
        if b == 1
            C = xor(C, shiftPoly(B, m, N+1));
        end
        m = m + 1;
    end
end

C = C(1:L+1);                            % 截断到 x^L
sufficient = (2*L <= N);                 % 数据充分性判定

if ~sufficient
    warning('数据不足: 2L=%d > N=%d，当前多项式欠定，请补采数据至 %d 位以上。', ...
            2*L, N, 2*L);
end
end

% ---- 辅助函数: 计算 x^m * B 并截断到指定长度 ----
function p = shiftPoly(B, m, len)
p = zeros(1, len);
p(1+m:min(1+m+find(B,1,'last')-1, len)) = B(1:min(1+find(B,1,'last')-1, len-m));
end