インパルス応答の測定方法の一つとして使われている M系列を題材にScilabやOctaveを使いこなしてみる。後の項目では、前の項目の関数に依存している部分もあるので、最初から順に実行する。
MLSの生成方法は他に色々なページで述べられているので、ここではScilabとOctaveでの一例を挙げる。 この例では、ScilabとOctaveの両方ともで同じ記述となった。
function list = mls(n)
switch n
case 2
taps=2; tap1=1; tap2=2;
case 3
taps=2; tap1=1; tap2=3;
case 4
taps=2; tap1=1; tap2=4;
case 5
taps=2; tap1=2; tap2=5;
case 6
taps=2; tap1=1; tap2=6;
case 7
taps=2; tap1=1; tap2=7;
otherwise
disp('input bits not supported (must be between 2~7)');
return
end
buff = ones(1,n);
list = zeros(1,2^n-1);
for i = (2^n)-1:-1:1
xorbit = bitxor(buff(tap1),buff(tap2));
buff = [xorbit buff(1:n-1)];
list(i) = xorbit;
end
endfunction
mls(3)
mls(4)
--> mls(3) ans = 1. 1. 1. 0. 0. 1. 0. --> mls(4) ans = 1. 1. 1. 1. 0. 0. 0. 1. 0. 0. 1. 1. 0. 1. 0.
function list = mls(n)
switch n
case 2
taps=2; tap1=1; tap2=2;
case 3
taps=2; tap1=1; tap2=3;
case 4
taps=2; tap1=1; tap2=4;
case 5
taps=2; tap1=2; tap2=5;
case 6
taps=2; tap1=1; tap2=6;
case 7
taps=2; tap1=1; tap2=7;
otherwise
disp('input bits not supported (must be between 2~7)');
return
end
buff = ones(1,n);
list = zeros(1,2^n-1);
for i = (2^n)-1:-1:1
xorbit = bitxor(buff(tap1),buff(tap2));
buff = [xorbit buff(1:n-1)];
list(i) = xorbit;
end
endfunction
mls(3)
mls(4)
>> mls(3) ans = 1 1 1 0 0 1 0 >> mls(4) ans = 1 1 1 1 0 0 0 1 0 0 1 1 0 1 0
MLSを順に並べて生成した行列をMLS matrixと呼ぶことがあり、アダマール行列(Hadamard matrix)と似た性質がある。こちらについては別項目を参照。
// only for n=2,3,4,5 function mmat = mlsmatrix(n) seq = mls(n); mmat = zeros((2^n)-1,(2^n)-1); for i = 1:1:(2^n)-1 mmat(i,:) = seq; seqd = bin2dec(strcat(string(seq))); seqd = bitand( bitor( seqd*2, int(seqd/2^(2^n-2)) ), (2^(2^n-1))-1); seq = strtod(strsplit(dec2bin( seqd ,(2^n)-1)))' end endfunction mlsmatrix (3)
--> mlsmatrix (3) ans = 1. 1. 1. 0. 0. 1. 0. 1. 1. 0. 0. 1. 0. 1. 1. 0. 0. 1. 0. 1. 1. 0. 0. 1. 0. 1. 1. 1. 0. 1. 0. 1. 1. 1. 0. 1. 0. 1. 1. 1. 0. 0. 0. 1. 1. 1. 0. 0. 1.
function mmat = mlsmatrix(n) seq = mls(n); mmat = zeros((2^n)-1); for i = 1:1:(2^n)-1 mmat(i,:) = seq; seqd = bin2dec(char(seq+'0')); seqd = bitand(bitor(bitshift(seqd,1),bitshift(seqd,(2^n-2)*-1)),(2^(2^n-1))-1); seq = dec2bin(seqd,(2^n)-1)-'0'; end endfunction mlsmatrix (3)
>> mlsmatrix (3) ans = 1 1 1 0 0 1 0 1 1 0 0 1 0 1 1 0 0 1 0 1 1 0 0 1 0 1 1 1 0 1 0 1 1 1 0 1 0 1 1 1 0 0 0 1 1 1 0 0 1
積が単位行列に近い形となる。これを利用して、畳み込みによりインパルス応答を演算できることがわかる。但し、DC成分が乗ってしまうのが注意点である。
mlsmatrix(3)*mlsmatrix(3) mlsmatrix(3)*mlsmatrix(3) / 2^(3-2) -1 norm(mlsmatrix(4)*mlsmatrix(4) / 2^(4-2) -1 - eye(2^4-1,2^4-1)) norm(mlsmatrix(5)*mlsmatrix(5) / 2^(5-2) -1 - eye(2^5-1,2^5-1))
ans = 4. 2. 2. 2. 2. 2. 2. 2. 4. 2. 2. 2. 2. 2. 2. 2. 4. 2. 2. 2. 2. 2. 2. 2. 4. 2. 2. 2. 2. 2. 2. 2. 4. 2. 2. 2. 2. 2. 2. 2. 4. 2. 2. 2. 2. 2. 2. 2. 4. ans = 1. 0. 0. 0. 0. 0. 0. 0. 1. 0. 0. 0. 0. 0. 0. 0. 1. 0. 0. 0. 0. 0. 0. 0. 1. 0. 0. 0. 0. 0. 0. 0. 1. 0. 0. 0. 0. 0. 0. 0. 1. 0. 0. 0. 0. 0. 0. 0. 1. ans = 0. ans = 0.
[d,e]=convol(mls(3),flipdim(mls(3),2)); [d,e]=convol(mls(3),flipdim(mls(3),2),e); disp(d([2^3-1,1:2^3-2])) [d,e]=convol(mls(4),flipdim(mls(4),2)); [d,e]=convol(mls(4),flipdim(mls(4),2),e); disp(d([2^4-1,1:2^4-2]))
4. 2. 2. 2. 2. 2. 2. 8. 4. 4. 4. 4. 4. 4. 4. 4. 4. 4. 4. 4. 4. 4.
MLSとその順序を逆にしたものを使って巡回畳み込みを行うと、インパルス応答が演算できる。時間領域で畳み込みを行うとかなり時間がかかるので、FFTを利用した高速な畳み込み方を使う事が多い。 その他に高速アダマール変換を利用する方法もあり、これは後述する。
testsig0 = [10,2,1,-1,0,5,2,1,-3,0,0,-1,8,1,-2]
mls4 = mls(4)*2-1; // 0/1 to -1/+1
[f,e]=convol(mls4,testsig0);
[f,e]=convol(mls4,testsig0,e); // testsig(*)mls
// deconvolution:
// convolve with backward sequence to get original signal
[g,e]=convol(f,flipdim(mls4,2));
[g,e]=convol(f,flipdim(mls4,2),e); // circular convolution
h = (g([2^4-1,1:2^4-2]) + sum(g)) / 2^4;
h = round(h) // omit arithmetic error
if ( testsig0 == h ) then
disp('MLS deconvolution was successful.');
else
disp('MLS deconvolution was unsuccessful.');
end
testsig0 = 10. 2. 1. -1. 0. 5. 2. 1. -3. 0. 0. -1. 8. 1. -2. h = 10. 2. 1. -1. 0. 5. 2. 1. -3. 0. 0. -1. 8. 1. -2. MLS deconvolution was successful.
アダマール変換を利用した方法を使う場合、信号順を置換してアダマール変換に合う順序として高速アダマール変換等で処理した後、また信号順を置換する必要がある。 詳細については下記論文を参照すると良い。置換の数列を２つ用意するが、片方はMLS matrixの上ビット分、もう片方は上ビット分が単位行列になるように MLS matrixから列を探してビット分並べると、その２つの行列の積がMLS matrixになることがわかる。
function [matr,matc] = mlsperm(n) mlsmat = mlsmatrix(n); matc = mlsmat(1:n,:); matr = zeros((2^n)-1,n); matu = eye(n,n); for i = 1:1:n for j = 1:1:(2^n)-1 if ( mlsmat(1:n,j) == matu(:,i) ) then matr(:,i) = mlsmat(:,j) end end end endfunction function rlist = mlspermr(matr,n) rlist = bin2dec(strsplit(strcat(string(matr')),[n*(1:(2^n)-2)]))' endfunction function clist = mlspermc(matc,n) clist = bin2dec(strsplit(strcat(string(matc)),[n*(1:(2^n)-2)]))' endfunction [mr,mc] = mlsperm(3) cc = modulo(mr*mc,2) permr = mlspermr(mr,3) permc = mlspermc(mc,3)
mc = 1. 1. 1. 0. 0. 1. 0. 1. 1. 0. 0. 1. 0. 1. 1. 0. 0. 1. 0. 1. 1. mr = 1. 0. 0. 0. 1. 0. 0. 0. 1. 1. 1. 0. 0. 1. 1. 1. 1. 1. 1. 0. 1. cc = 1. 1. 1. 0. 0. 1. 0. 1. 1. 0. 0. 1. 0. 1. 1. 0. 0. 1. 0. 1. 1. 0. 0. 1. 0. 1. 1. 1. 0. 1. 0. 1. 1. 1. 0. 1. 0. 1. 1. 1. 0. 0. 0. 1. 1. 1. 0. 0. 1. permr = 4. 2. 1. 6. 3. 7. 5. permc = 7. 6. 4. 1. 2. 5. 3.
先述のMLSの逆順を畳み込むことと、MLS matrixとの行列積を求めることは同義で、行列積を高速に処理したい時は例えばATLAS等のBLASのDGEMVを使ったりしても良いが、 先の論文で提案されているように、MLS matrixは上で求めた置換用の数列を使って変形することでアダマール行列にすることができるので、入力ベクトルの順序を置換すれば、 バタフライ演算を使用してO(n log(n))の計算量で処理できるようになる。 変形にはいくつか手順が必要で、アダマール行列の1と-1は0と1に置換し、MLS行列の1行目と1列目には0を追加する。ScilabやMATLABやOctaveでは、行列の置換の記述方法が非常に楽で、 左辺に置換数列を指定することもできるので、このような時は使いやすいと思う。MATLABはMatrix Laboratoryの略というのがよく分かる。
function hmatrix=hadamard(n) h=[1,1;1,-1] hmatrix=h; for n = 1:1:(log2(n)+1)-2 hmatrix = hmatrix.*.h; end endfunction mls3matex = [0,zeros(1,7);zeros(7,1)mlsmatrix(3)] (hadamard(8)-1)/-2 // Hadamard matrix -> MLS matrix permutation hadama8p1 = ((hadamard(8)-1)/-2)([1,permr+1],[1,permc+1]) hadama8p2 = ((hadamard(8)-1)/-2)([1,permc+1],[1,permr+1]) mls3matex - hadama8p1 mls3matex - hadama8p2
mls3matex = 0. 0. 0. 0. 0. 0. 0. 0. 0. 1. 1. 1. 0. 0. 1. 0. 0. 1. 1. 0. 0. 1. 0. 1. 0. 1. 0. 0. 1. 0. 1. 1. 0. 0. 0. 1. 0. 1. 1. 1. 0. 0. 1. 0. 1. 1. 1. 0. 0. 1. 0. 1. 1. 1. 0. 0. 0. 0. 1. 1. 1. 0. 0. 1. --> (hadamard(8)-1)/-2 ans = 0. 0. 0. 0. 0. 0. 0. 0. 0. 1. 0. 1. 0. 1. 0. 1. 0. 0. 1. 1. 0. 0. 1. 1. 0. 1. 1. 0. 0. 1. 1. 0. 0. 0. 0. 0. 1. 1. 1. 1. 0. 1. 0. 1. 1. 0. 1. 0. 0. 0. 1. 1. 1. 1. 0. 0. 0. 1. 1. 0. 1. 0. 0. 1. --> mls3matex - hadama8p1 ans = 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. --> mls3matex - hadama8p2 ans = 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.
// MLS matrix -> Hadamard matrix permutation lmr([1,permr+1]) = 1:8 lmc([1,permc+1]) = 1:8 mls3matex(lmr,lmc) - (hadamard(8)-1)/-2 mls3matex(lmc,lmr) - (hadamard(8)-1)/-2
> lmr([1,permr+1]) = 1:8 lmr = 1. 4. 3. 6. 2. 8. 5. 7. --> lmc([1,permc+1]) = 1:8 lmc = 1. 5. 6. 8. 4. 7. 3. 2. --> mls3matex(lmr,lmc) - (hadamard(8)-1)/-2 ans = 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. --> mls3matex(lmc,lmr) - (hadamard(8)-1)/-2 ans = 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.
// [2^3] impulse input n = 3; [mr,mc] = mlsperm(n); permr = mlspermr(mr,n); permc = mlspermc(mc,n); mlsmat = mlsmatrix(n); hadamardMatrix = hadamard(2^n); sigmlsperm = zeros(1,2^n-1) sigmlsperm(permc) = mlsmat(1,:) hinput = [0,sigmlsperm] houtput = hinput*hadamardMatrix/(2^n) hout = houtput(2:2^3) signal = hout(permr)*(-2)
--> sigmlsperm = zeros(1,2^n-1) sigmlsperm = 0. 0. 0. 0. 0. 0. 0. --> sigmlsperm(permc) = mlsmat(1,:) sigmlsperm = 0. 0. 0. 1. 1. 1. 1. --> hinput = [0,sigmlsperm] hinput = 0. 0. 0. 0. 1. 1. 1. 1. --> houtput = hinput*hadamardMatrix/(2^n) houtput = 0.5 0. 0. 0. -0.5 0. 0. 0. --> hout = houtput(2:2^3) hout = 0. 0. 0. -0.5 0. 0. 0. --> signal = hout(permr)*(-2) signal = 1. 0. 0. 0. 0. 0. 0.
// Test signal 1 testsig1 = [1,-1,0.5,-0.5,0.1,-0.1,0] input = testsig1*mlsmat // convolution input sigmlsperm = zeros(1,2^n-1); sigmlsperm(permc) = input hinput = [0,sigmlsperm] houtput = hinput*hadamardMatrix/(2^n) hout = houtput(2:2^3) signal1 = hout(permr)*(-2)
--> testsig1 = [1,-1,0.5,-0.5,0.1,-0.1,0] testsig1 = 1. -1. 0.5 -0.5 0.1 -0.1 0. --> input = testsig1*mlsmat // convolution input input = 0.4 0.1 0.4 0.5 -1.5 1.1 -1. --> sigmlsperm(permc) = input sigmlsperm = 0.5 -1.5 -1. 0.4 1.1 0.1 0.4 --> hinput = [0,sigmlsperm] hinput = 0. 0.5 -1.5 -1. 0.4 1.1 0.1 0.4 --> houtput = hinput*hadamardMatrix/(2^n) houtput = 0. -0.25 0.5 -0.05 -0.5 0. 0.25 0.05 --> hout = houtput(2:2^3) hout = -0.25 0.5 -0.05 -0.5 0. 0.25 0.05 --> signal1 = hout(permr)*(-2) signal1 = 1. -1. 0.5 -0.5 0.1 -0.1 0.
// Test signal 2 testsig2 = [4,-3,3,-9,-3,2,1] input = testsig2*mlsmat // convolution input sigmlsperm = zeros(1,2^n-1); sigmlsperm(permc) = input hinput = [0,sigmlsperm] houtput = hinput*hadamardMatrix/(2^n) hout = houtput(2:2^3) signal2 = hout(permr)*(-2)
--> testsig2 = [4,-3,3,-9,-3,2,1] testsig2 = 4. -3. 3. -9. -3. 2. 1. --> input = testsig2*mlsmat // convolution input input = 6. -1. -2. 3. -13. -5. -8. --> sigmlsperm(permc) = input sigmlsperm = 3. -13. -8. -2. -5. -1. 6. --> hinput = [0,sigmlsperm] hinput = 0. 3. -13. -8. -2. -5. -1. 6. --> houtput = hinput*hadamardMatrix/(2^n) houtput = -2.5 -1.5 1.5 1.5 -2. -0.5 4.5 -1. --> hout = houtput(2:2^3) hout = -1.5 1.5 1.5 -2. -0.5 4.5 -1. --> signal2 = hout(permr)*(-2) signal2 = 4. -3. 3. -9. -3. 2. 1.
// [2^4]
n = 4
[mr,mc] = mlsperm(n)
permr = mlspermr(mr,n)
permc = mlspermc(mc,n)
hadamardMatrix = hadamard(2^n);
rand('seed',0)
testsig3 = rand(1,2^n-1)
input = testsig3 * mlsmatrix(n);
sigmlsperm = zeros(1,2^n-1);
sigmlsperm(permc) = input;
hinput = [0,sigmlsperm];
houtput = hinput*hadamardMatrix/(2^n);
hout = houtput(2:2^n);
signal3 = hout(permr)*(-2)
if ( testsig3 == signal3 ) then
disp('MLS deconvolution was successful.');
else
disp('MLS deconvolution was unsuccessful.');
end
testsig3 = column 1 to 9 0.2113249 0.7560439 0.0002211 0.3303271 0.6653811 0.6283918 0.8497452 0.685731 0.8782165 column 10 to 15 0.068374 0.5608486 0.6623569 0.7263507 0.1985144 0.5442573 signal3 = column 1 to 9 0.2113249 0.7560439 0.0002211 0.3303271 0.6653811 0.6283918 0.8497452 0.685731 0.8782165 column 10 to 15 0.068374 0.5608486 0.6623569 0.7263507 0.1985144 0.5442573 MLS deconvolution was successful.
But, for my own part, it was Greek to me.
-- William Shakespeare, "Julius Caesar"
The difference between $399 and $220, i.e., $179, is called consumer surplus.
It's the extra value that those rich consumers got from their purchase that
they would have been perfectly happy to do without.
It's sort of like if you were all set to buy that new merino wool sweater, and
you thought it was going to cost $70, which is well worth it, and when you got
to Banana Republic it was on sale for only $50! Now you have an extra $20 in
found money that you would have been perfectly happy to give to the Banana
Republicans!
Yipes!
That bothers good capitalists. Gosh darn it, if you're willing to do without
it, well, give it to me! I can put it to good use, buying a SUV or condo or
Mooney or yacht one of those other things capitalists buy!
In economist jargon, capitalists want to capture the consumer surplus.
-- Joel Spolsky
-- "Camels and Rubber Duckies" ( http://www.joelonsoftware.com/articles/CamelsandRubberDuckies.html )