M序列的生成

matlab就提供了产生M序列的专门函数

idinput函数

产生系统辨识经常使用的典型信号。

格式：u = idinput(N,type,band,levels)

[u,freqs] = idinput(N,'sine',band,levels,sinedata)

N

产生的序列的长度，若是N=[N nu]，则nu为输入的通道数，若是N=[P nu M]，则nu指定通道数，P为周期，M*P为信号长度。默认状况下，nu=1，M=1，即一个通道，一个周期。

Type

指定产生信号的类型，可选类型以下

‘rgs’

高斯随机信号

‘rbs’ （默认）

二值随机信号

‘prbs’

二值伪随机信号（M序列）

‘sine’

正弦信号和

Band

指定信号的频率成分。对于’rgs’、’rbs’、’sine’，band = [wlow, whigh]指定通带的范围，若是是白噪声信号，则band=[0, 1]，这也是默认值。指定非默认值时，至关于有色噪声。

对于’prbs’，band=[0, B]，B表示信号在一个间隔1/B（时钟周期）内为恒值，默认为[0, 1]。

Levels

指定输入的水平。Levels=[minu, maxu]，在type=’rbs’、’prbs’、’sine’时，表示信号u的值老是在minu和maxu之间。对于type=’rgs’，minu指定信号的均值减标准差，maxu指定信号的均值加标准差，对于0均值、标准差为1的高斯白噪声信号，则levels=[-1, 1]，这也是默认值。

说明

对于PRBS信号，若是M>1，则序列的长度和PRBS周期会作调整，使PRBS的周期为对应必定阶数的最大值（即2^n-1，n为阶数）；若是M=1，PRBS的周期是大于N的相应阶数的值。在多输入的情形时，信号被最大平移，即P/nu为此信号能被估计的模型阶次的上界。

上面的意思可以下理解：对于M=1时，

ms = idinput(12, 'prbs', [0 1], [0 1]);

figure

stairs(ms)

title('M序列')

ylim([-0.5 1.5])

ms = idinput(12, 'prbs', [0 1], [0 1]);

figure

stairs(ms)

title('M序列')

ylim([-0.5 1.5])

结果以下

同时，matlab给出以下警告

Warning: The PRBS signal delivered is the 12 first values of a full sequence of length 15.

即函数的输出为周期为15（大于12的第一个2^n-1的值）PRBS信号的前12个值组成的序列。如

ms = idinput(15, 'prbs', [0 1], [0 1]);

figure

stairs(ms)

title('M序列')

ylim([-0.5 1.5])

能够看到指定12时的序列为指定15时的序列的前面部分。

对于M>1时，

ms = idinput([12,1,2], 'prbs', [0 1], [0 1]);

figure

stairs(ms)

title('M序列')

ylim([-0.5 1.5])

结果以下

Matlab给出的响应警告为

Warning: The period of the PRBS signal was changed to 7. Accordingly, the length of the

generated signal will be 14.

对于正弦信号和的产生，貌似用的很少，语法还挺复杂，等用的时候再看吧。

方法

产生’rgs’信号的带通讯号使用的是一个8阶巴特沃斯滤波器，使用idfilt作的非因果滤波，这个是可信赖的方法。

对于’rbs’信号，使用的是相同的滤波器，可是是在二值化以前，这意味着频率成分并不保证是精确的。

产生高斯随机信号

clc

clear all

close all

% 高斯随机信号

u = idinput(1000, 'rgs');

figure

stairs(u)

title('高斯随机信号')

figure

hist(u, -4:4)

title('高斯随机信号的分布')

产生二值随机信号

clc

clear all

close all

% 二值随机信号

u = idinput(100, 'rbs');

figure

stairs(u)

title('二值随机信号')

ylim([-1.5 1.5])

产生二值伪随机信号（M序列）

合理的选择输入激励信号，能有效的激励起系统的动态信号。白噪声的平稳谱的性质决定了它是一个很好的输入信号，但它在工程中不易实现，而M序列具备近似白噪声的性质，可保证良好的辨识精度。

clc

clear all

close all

% 二值伪随机信号（M序列）

n = 8; % 阶次

p = 2^n -1; % 循环周期

ms = idinput(p, 'prbs');

figure

stairs(ms)

title('M序列')

ylim([-1.5 1.5])

结果

验证M序列的性质以下

• -1和1的个数差1

sum(ms==1) % 1的个数

sum(ms==-1) % -1的个数

ans =127

ans =128

• 存在直流份量

mean(ms) % 直流份量

ans =-0.0039

• 相关函数

a = zeros(length(ms)*10, 1); % 采样

for i = 1:10

a(i:10:end) = ms;

end

c = xcorr(a, 'coeff'); % 自相关函数

figure

plot(c)

title('相关函数')

自相关函数接近于δ函数。

• 谱密度

figure

pwelch(a) % 谱密度

说明M序列不含基频的整数倍的频率成分。

产生逆M序列

谱分析代表，M序列含有直流成分，将形成对辨识系统的“净扰动”，这一般不是所但愿的。而逆M序列将克服这一缺点，是一种比M序列更为理想的伪随机码序列。

clc

clear all

close all

% 二值伪随机信号（M序列）

n = 8; % 阶次

p = 2^n -1; % 循环周期

ms = idinput(p, 'prbs', [], [0 1]);

figure

stairs(ms)

title('M序列')

ylim([-0.5 1.5])

% 产生逆M序列

s = 0;

ims = zeros(2*p, 1);

mstemp = [ms; ms];

for i = 1:2*p

ims(i) = xor(mstemp(i), s);

s = not(s);

end

ims(ims==0) = -1;

figure

stairs(ims)

title('逆M序列')

ylim([-1.5 1.5])

• -1和1的个数差1

sum(ims==1) % 1的个数

sum(ims==-1) % -1的个数

ans =

255

ans =

255

• 无直流份量

mean(ims) % 直流份量

ans =

0

• 相关函数

a = zeros(length(ims)*10, 1); % 采样

for i = 1:10

a(i:10:end) = ims;

end

c = xcorr(a, 'coeff');

figure

plot(c)

• 谱密度

figure

pwelch(a)

 

转自：http://blog.sina.com.cn/s/blog_6163bdeb0102e2qj.html