インパルス応答の測定方法の一つとして使われている 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.
A classic is something that everyone wants to have read
and nobody wants to read.
-- Mark Twain, "The Disappearance of Literature"
<pmurias> rindolf: what is Park/Spark?
<rindolf> pmurias:
http://www.shlomifish.org/open-source/projects/Park-Lisp/
<rindolf> pmurias: it's still incomplete.
<rindolf> And I haven't updated it.
<pmurias> rindolf: if you like lisp/perl6 projects you might consider
helping with a common lisp elf backend
<rindolf> pmurias: Common Lisp.
<rindolf> pmurias: thing is I think both CL and Scheme suck.
<rindolf> I like Lisp as a concept.
<rindolf> Arc is nice, but has too many implementation problems.
<rindolf> And missing features.
<vixey> Arc is not nice
<rindolf> I want to give a presentation to the Perl Mongers about
"Foreign Languages: Lisp"
<rindolf> vixey: I like it.
<rindolf> Though I hate that "(not)" has become "(no)"
<rindolf> it's so non-English.
<vixey> it's just TCL with horrible syntax
<rindolf> vixey: but it's missing a lot of exciting features.
<rindolf> Which PG deemed as unnecessary.
<rindolf> Doesn't look like the 100-years language to me.
<rindolf> Which is why - Spark!
<pmurias> why not just write an s-expression p6 dialect?
<rindolf> pmurias: could be.
<rindolf> pmurias: it's another approach.
<rindolf> But some things make sense in Lisp and not in p6.
<rindolf> For example, Perl does not like to use + for string or list
concat.
<rindolf> While Python does and it seems to be OK in Arc too.
<rindolf> And in CL you have (concatenate) (yuck!).
<pbuetow> (((hehe)))
<pmurias> + for strings sucks
<Auzon> seconded.
<vixey> rindolf: If you don't like CONCATENATE you can just rename it
<rindolf> vixey: yeah.
<rindolf> vixey: but I'd rather not rename concatenate because then
people won't understand my code.
<rindolf> vixey: as TimToady said people hate abstractions.
<vixey> yes they will rindolf
<rindolf> They want things to work out of the box.
<vixey> A program is many many totally newly defined procedures
<vixey> just renaming one thing is nothing in the context of a big
program
<rindolf> vixey: "let's spend 3 days creating a new language, and 1 day
implementing the solution with it."
<TimToady> if it would take 10 days without the new language, it's worth
it
<rindolf> TimToady: yeah.
<rindolf> TimToady: but this is the CL mentality.
<vixey> no it's not
<rindolf> Sometimes you can take 1 day to write an API.
<vixey> CL is too diverse you cannot generalize like that
<rindolf> vixey: I meant a common idiom there.
<rindolf> I think I'll /quit and do something productive.
<rindolf> Like work on Spark.
<vixey> another quote:
<vixey> how to write any computer program in two easy stages:
<vixey> Design and implement the programming language which would be
best for solving the problem.
<vixey> Write the program in the language you���ve just implemented.
<rindolf> vixey: or just use Perl which is the best for everything.
<vixey> heh
<TimToady> the second step is obvious--the best language for the job is
one that does the job on a null input
* pmichaud notes that vixey's algorithm is somewhat recursive
<TimToady> "All rules of thumb are false, including this one."
<pmurias> rindolf: when you feel like writing Common Lisp backends,
contact me or mncharity ;)
-- Lisp Dialects (Scheme, Common Lisp, Arc, Spark) Mentality and Usability
-- #perl6, Freenode