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