* *
M序列是工程中常用的输入信号,它的性质类似于白噪声,而白噪声是理论上最好的输入信号,可见M序列的价值。下面介绍M序列的matlab产生方法。
idinput函数
产生系统辨识常用的典型信号。
格式
N
u = idinput(N,type,band,levels)
[u,freqs] = idinput(N,'sine',band,levels,sinedata)
产生的序列的长度,如果N=[N nu],则nu为输入的通道数,如果N=[P nu M],则
nu指定通道数,P为周期,M*P为信号长度。默认情况下,nu=1,M=1,即一个通道,一个周期。
Type
指定产生信号的类型,可选类型如下
‘rgs’ ‘rbs’ (默认) ‘prbs’ ‘sine’ 高斯随机信号 二值随机信号 二值伪随机信号(M序列) 正弦信号和 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])
结果如下
* *
同时,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])
* *
对于M>1时,
可以看到指定12时的序列为指定15时的序列的前面部分。
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)
* *