インパルス応答の測定方法の一つとして使われている 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.
It's not true that those who can't do, teach (some of the best hackers I know
are professors), but it is true that there are a lot of things that those who
teach can't do. Research imposes constraining caste restrictions. In any
academic field there are topics that are ok to work on and others that aren't.
Unfortunately the distinction between acceptable and forbidden topics is
usually based on how intellectual the work sounds when described in research
papers, rather than how important it is for getting good results. The extreme
case is probably literature; people studying literature rarely say anything
that would be of the slightest use to those producing it.
Though the situation is better in the sciences, the overlap between the kind
of work you're allowed to do and the kind of work that yields good languages
is distressingly small. (Olin Shivers has grumbled eloquently about this.) For
example, types seem to be an inexhaustible source of research papers, despite
the fact that static typing seems to preclude true macros-- without which, in
my opinion, no language is worth using.
-- Paul Graham
-- "The Hundred-Year Language" ( http://www.paulgraham.com/hundred.html )
<rindolf> What should I do now?
<rindolf> I'll work on Text-Qantor.
<rindolf> It's so great not to have a job.
<Zuu> yeah, if someone else pays for the food it sure is :D
<Zuu> also, i dont really understand much of what you just told me
:P
* Zuu puts a stick into the Text-Qantor
<rindolf> Zuu: Qantor == Qantor ain't no TeX/Troff oh really.
<rindolf> It's a typesetting system I'm working on.
* Zuu hates the name
<Zuu> it makes me kinda mad actually :/
<rindolf> Zuu: :-)
<rindolf> Zuu: maybe it will grow on you.
<rindolf> Zuu: some people I know named a browser suckass.
<Zuu> :(
<rindolf> I refused to work on it.
<Zuu> see that's a name!
<rindolf> Zuu: heh.
<Zuu> i didnt mean that btw :)
<Zuu> suckass is kinda... unkind
<rindolf> OK, now I should write an
http://www.shlomifish.org/humour/bits/facts/XSLT/
transformation.
<rindolf> I'll start from something I already have.
<Zuu> But the "X ain't no <something related>" is just a lame naming
convention imho
<Zuu> yeah, work on some XSLT facts :D
<rindolf> Zuu: just call it Qantor then.
<rindolf> Without the mnemonics.
<Zuu> but anyone interrested will learn that it's an abbreviation
<Zuu> just by the fact that it's recursive makes me want to kill
myself a little bit more :P
<rindolf> Zuu: do me a break and kill yourself.
<Zuu> :>
<rindolf> Less Zuus - more grass for evil reindeers like me to feed on.
-- What is Qantor?
-- ##programming, Freenode