M序列
M 序列是一种伪随机序列,具有很好的伪噪声特性,常用于信道噪声测试和保密通信。不过 M 序列还有一个用途,也就是本文所介绍的——通过 M 序列测量频率响应。在讨论这个问题之前,我们先介绍 M 序列的特征与生成方法。
M 序列通过 N 阶线性反馈移存器生成,其周期为 2 N − 1 2^N-1 2N−1 。
如上图所示,反馈到最左端的信号为
a
n
=
∑
i
=
1
n
c
i
a
n
−
i
(
m
o
d
2
)
a_n=\sum_{i=1}^{n}c_ia_{n-i}\ (mod\ 2)
an=i=1∑ncian−i (mod 2)
经过
k
k
k 次移位后的输出为
a
k
=
∑
i
=
1
n
c
i
a
k
−
i
(
m
o
d
2
)
a_k=\sum_{i=1}^{n}c_ia_{k-i}\ (mod\ 2)
ak=i=1∑nciak−i (mod 2)
用特征多项式表示反馈的连接状态
f
(
x
)
=
c
n
x
n
+
c
n
−
1
x
n
−
1
+
⋯
+
c
1
x
+
c
0
=
∑
i
=
0
n
c
i
x
i
f(x)=c_nx_n+c_{n-1}x^{n-1}+\dots+c_1x+c_0=\sum_{i=0}^nc_ix^i
f(x)=cnxn+cn−1xn−1+⋯+c1x+c0=i=0∑ncixi
由于连接状态中,
c
0
、
c
n
c_0、c_n
c0、cn 均为 1,因此特征多项式必然包含
x
n
+
1
x^n+1
xn+1 ,特征多项式的次数 n 即是移存器的级数。
输出序列用多项式表示(母函数)
G
(
x
)
=
a
0
+
a
1
x
+
a
2
x
2
+
⋯
=
∑
i
=
0
∞
a
i
x
i
G(x)=a_0+a_1x+a_2x^2+\cdots=\sum_{i=0}^\infin a_ix^i
G(x)=a0+a1x+a2x2+⋯=i=0∑∞aixi
本原多项式:若一个多项式满足:1.
f
(
x
)
f(x)
f(x) 是既约的;2.
f
(
x
)
f(x)
f(x) 可以整除
2
m
+
1
2^m+1
2m+1 ,
m
=
2
n
−
1
m=2^n-1
m=2n−1 ;3.
f
(
x
)
f(x)
f(x) 除不尽
2
q
+
1
2^q+1
2q+1 ,
q
<
m
q<m
q<m ;则称该多项式为本原多项式。
当特征多项式为本原多项式时,反馈移存器的输出序列为 M 序列。
matlab 实现
% M 序列生成器
% @mg: M序列
% @f : 本原多项式
% a(k)=mod(sum_1^n(c(i)*a(k-i)),2)
function mg=m_generator(f)
f=dec2bin(f)-'0';
n=length(f)-1; %M序列阶数
N=2^n-1; %伪随机码的周期
register=[zeros(1,n-1),1];
newregister=zeros(1,n);
mg=zeros(1,N);
f_inv=f(end-1:-1:1);
for i=1:N
newregister(1)=mod(sum(f_inv.*register),2);
newregister(2:n)=register(1:n-1);
register=newregister;
mg(i)=register(n);
end
end
% 测试
clc,clear,close all
N=10; % 阶数
L=2^N-1;
F=primpoly(N,'all'); %显示该阶数下的全部本原多项式
mg=m_generator(F(1)); %使用一个本原多项式生成m序列
M序列测量滤波器频率响应
通常测量某一系统的幅频响应,是通过扫频来完成的,然而需要时间较长,还需要实现 DDS,比较复杂。而考虑到 m 序列具有白谱,因此测量经过滤波器后的 m 序列的值,并做 DFT,即可获知滤波器的幅频特性。
假设生成 10 阶 M 序列(1023 点),以 256 k H z 256kHz 256kHz 的比特速率发送该序列,经过滤波器后,以同样的采样率用 ADC 采样,则频率分辨率为 256 k H z / 1023 = ∼ 250 H z 256kHz/1023=\sim 250Hz 256kHz/1023=∼250Hz 。
%---------m序列--------
clc,clear,close all
N=10;
L=2^N-1;
F=primpoly(N,'all'); %显示该阶数下的全部本原多项式
mg=m_generator(F(1))'*2-1;
%% m序列幅频特性
fs=256e3;
freq=fs*(0:(L/2))/L;
mg_fft=fft(mg)/L*2;
plot(freq/1e3,abs(mg_fft(1:floor(L/2)+1))) %幅频显示为白谱
xlabel('f/kHz')
%% 滤波
mg2=smooth(mg);
mg2_fft=fft(mg2)/L*2;
figure
plot(freq/1e3,abs(mg2_fft(1:floor(L/2)+1))) %滤波器幅频特性曲线
xlabel('f/kHz')
测试中使用了平滑滤波器,测得的滤波器频率响应如下