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)
