MATLAB m 序列生成函数 & 相关函数

观摩了 Simon Haykin 《Communication Systems 4th Edition》第七章课后的代码，感觉有点谜人，做了注释及部分修改。

PN 序列生成函数

function x = PNseq(p)
% 根据多项式 p 生成 PN 序列的线性移位寄存器
N = length(p) - 1; % 不包括 0 次项
p = fliplr(p);
X = [1 zeros(1, N - 1)]; % 初始状态
n = 1; % 输出的周期数
% 参考程序输出序列是反向序列，这里按教材顺序输出
for i = 1 : n * (2 ^ N - 1)
    x(i) = X(N); % 输出
    X = [p(1) * rem(sum(p(2 : N + 1) .* X(1 : N)), 2) X(1 : N - 1)];
end   
end

相关函数

function [corrf] = pn_corr(u, v, N)
% 求 u, v 相关函数，不做归一化处理
% 处理循环位移的 N 种情况
for m = 0 : N - 1
    shifted_u = [u(m + 1 : N) u(1 : m)];
    corr(m + 1) = sum(v .* shifted_u);
end
% 负方向的相关函数是对称的
corrf = [corr(2 : N) corr];
end

函数使用

考虑周期为 N = 63 的两个 m 序列，反馈抽头为 [6, 5, 2, 1], [6, 5, 4, 1]，求互相关函数。

% 反馈抽头对应的多项式 x^6+x^5 +x^2+x+1 & x^6+x^5 +x^4+x+1
pol1 = [1 1 0 0 1 1 1];
pol2 = [1 1 1 0 0 1 1];
% 码片周期
N = 63;
% 生成 m 序列
pnseq1 = PNseq(pol1);
pnseq2 = PNseq(pol2);
% 将 01 序列映射为电平表示
u = 2 * pnseq1 - 1;
v = 2 * pnseq2 - 1;
% 求互相关函数
[c_corrf] = pn_corr(u, v, N);
plot(-62 : 62, c_corrf(1 : 125)); axis([-62, 62, -20, 20]);
xlabel('Delay \tau'); ylabel('Cross function R_{c}(\tau)');