Hi Søren,
To which package do you plan on
adding your functions?
Communication package
To keep the number of inactive developers down we have a policy that new
developers should post code on the list before they are given access to
SVN. Could you post you code once you're ready? Perhaps somebody can
even help with the speed issues :-)
Sure. I attached the code. I didn't do it because it isn't publishable
code now and some comments would be nice (I will not commit the attached
code as it is!). I'm not so confident regarding the speed issues. I
vectorized the main bottleneck, but the problem is, that you have two
different task, if you uses a Viterbi algorithm. This task can only be
calculated interleaved or you blow up the memory.
whats the problem?:
- first, you must calculate a path metric (it indicates the
probabilities of the code words),
- second you must decode the most probable message,
then you must start again with the next bit (encoder input word).
Sure, i'm not the guru of fast octave programming, but i think it's
nearly impossible to make it really fast in octave only. There is a big
big calculation efford and it's hard to vectorize the code, because you
have data dependencies and you need many iterations (It's a little bit
like writing a computer chess is octave).
If you write C functions, it can be much faster, because you lost the
interpreting overhead. I'm sure this is the way to do. But the fist step
should be working scripts, that will be used as reference for C
implementation later.
regards,
Martin
## -*- texinfo -*-
## @deftypefn {Function File} {...@var{code} =} vitdec (@var{message},@var{trellis},@var{punctPattern})
##
## Return the @var{code} of the @var{message} coded with given @var{trellis} description and punctured with the optional @var{punctPattern}.
##
## @seealso{poly2trellis, vitdec}
## @end deftypefn
function code = convenc(msg, t, punctPattern);
if (nargin < 2 || nargin > 3)
print_usage();
end
if (nargin == 3)
error("Sorry, puncturing isn't implemented at the moment!");
end
len = length(msg);
state = 0;
numCodeBits = floor(log2(t.numOutputSymbols));
code = zeros(numCodeBits, len);
for n=1:length(msg)
output = t.outputs(state + 1, msg(n) + 1);
for m=1:numCodeBits
code(m, n) = bitand(bitshift(output, m - numCodeBits), 1);
end
state = t.nextStates(state + 1, msg(n) + 1);
end
endfunction
function trel = poly2trellis(influence, poly, feedback);
if nargin < 2 || nargin > 3
error("need two or three arguments!");
end
if !isscalar(influence)
error("Implementation isn't complete! Only one encoder input supported now!");
end
if influence <= 1
error("influence must be grater than one!");
end
if rows(poly) != 1
error("Implementation isn't complete! poly must look like this [7 5] now!");
end
if nargin == 2
feedback = zeros(1,1);
else
if rows(feedback) != 1
error("Implementation isn't complete! feedback must look like this [357] now!");
end
end
numStates = 2^(influence - 1);
trel = struct('numInputSymbols', 2, ...
'numOutputSymbols', 2 ^ columns(poly), ...
'numStates', numStates,...
'nextStates', zeros(numStates, 2),...
'outputs', zeros(numStates, 2));
decpoly = oct2dec(poly);
outputValues = 2.^(columns(poly) - 1:-1:0);
decfeedback = oct2dec(feedback);
for n = 0:numStates-1
feedbackValue = mod(sum(bitget(bitand(decfeedback, n),1:influence-1)),2);
trel.nextStates(n+1,1) = floor((feedbackValue * numStates + n) / 2);
trel.nextStates(n+1,2) = floor(((1 - feedbackValue) * numStates + n) / 2);
temp0 = bitand(decpoly, n);
next0 = zeros(1,columns(temp0));
%keyboard
while sum(temp0) > 0
next0 = next0 + bitand(temp0, 1);
temp0 = bitshift(temp0, -1);
end
temp1 = bitand(decpoly, n + numStates);
next1 = zeros(1,columns(temp1));
while sum(temp1) > 0
next1 = next1 + bitand(temp1, 1);
temp1 = bitshift(temp1, -1);
end
trel.outputs(n+1,1) = mod(next0, 2) * outputValues';
trel.outputs(n+1,2) = mod(next1, 2) * outputValues';
end
function bits = vitdec2(trellis, symbols, tblen, mode, N);
t = trellis;
s = symbols;
numInputBits = round(log2(t.numOutputSymbols));
numOutputSymbols = t.numInputSymbols;
metricTime = 0;
tbTime = 0;
completeTime = time;
if (rows(s) != numInputBits)
error("Your trellis description need %d input bits, i got %d!", numInputBits, rows(s));
end
switch (mode)
case 'hard'
if sum((s(:) == 0) + (s(:) == 1)) != rows(s(:))
error("Symbols must be 0 or 1");
end
case 'soft'
if (min(s(:)) < 0) || (max(s(:)) > (N - 1))
error("Symbols must be in range from 0 to N-1");
end
otherwise
error("Mode must consists of \"soft\" or \"hard\"");
endswitch
expectedBits = zeros(rows(t.outputs), columns(t.outputs), numInputBits);
for m = 1:rows(t.outputs)
for n = 1:numOutputSymbols
symbol = t.outputs(m, n);
for o = 1:numInputBits
expectedBits(m, o, n) = (N - 1) * bitand(bitshift(symbol, o - numInputBits), 1);
% expectedBits(m, n, numInputBits - o + 1) = bitand(bitshift(symbol, (1 - o)), 1);
end
end
end
lastStates = zeros(rows(t.nextStates), 2) - 1;
for m = 1:rows(t.nextStates)
if (lastStates(t.nextStates(m,1) + 1, 1) < 0)
lastStates(t.nextStates(m,1) + 1, 1) = m; % should be m - 1 but we need m later, because octave starts indexing with 1
else
lastStates(t.nextStates(m,1) + 1, 2) = m; % the same
end
if (lastStates(t.nextStates(m,2) + 1, 1) < 0)
lastStates(t.nextStates(m,2) + 1, 1) = m; % the same
else
lastStates(t.nextStates(m,2) + 1, 2) = m; % the same
end
end
metric = zeros(t.numStates, tblen + 1);
bits = [];
dataLen = columns(s);
numInputBits = rows(s);
tbWindowStart = 0;
idx = 1;
while tbWindowStart < dataLen
tempTime = time;
metric = [metric(:,2:tblen + 1) zeros(t.numStates, 1)];
while idx <= tblen
inBits = getInputBitVector(s, tbWindowStart + idx);
for x=1:t.numStates
in(:,x) = inBits;
end
%in = inBits
for outputSymbol = 1:numOutputSymbols;
%for
state = 1:t.numStates;
%ex0 = reshape(expectedBits(state, outputSymbol, :)(:),t.numStates,numInputBits)'
ex0 = expectedBits(state, :, outputSymbol)';
%keyboard;
ns0 = t.nextStates(state, outputSymbol) + 1;
sm0 = metric(state, idx) + sum((N - 1) - abs(ex0 - in))';
%for tempidx = 1:t.numStates
% metric(ns0(tempidx), idx + 1) = max(metric(ns0(tempidx), idx + 1), sm0(tempidx));
%end
metric(ns0(1:2:end), idx + 1) = max(metric(ns0(1:2:end), idx + 1), sm0(1:2:end));
metric(ns0(2:2:end), idx + 1) = max(metric(ns0(2:2:end), idx + 1), sm0(2:2:end));
%end
end
idx = idx + 1;
end
metricTime = metricTime + (time - tempTime);
%metric
tempTime = time;
currentState = find(metric(:,tblen + 1) == max(metric(:,tblen + 1)))(1);
states = zeros(1,tblen + 1);
states(tblen + 1) = 0;
for idx=tblen + 1:-1:2
ls0 = lastStates(currentState, 1);
ls1 = lastStates(currentState, 2);
ls0v = metric(ls0, idx - 1);
ls1v = metric(ls1, idx - 1);
if (ls0v >= ls1v)
currentState = ls0;
else
currentState = ls1;
end
end
if (tbWindowStart > 0)
bits = [bits bitshift(currentState - 1, 1 - log2(t.numStates))];
end
tbTime = tbTime + (time - tempTime);
tbWindowStart = tbWindowStart + 1;
idx = tblen;
end
completeTime = time - completeTime;
%printf("complete time: %f s metric: %f s tb: %f s\n", completeTime, metricTime, tbTime);
endfunction
function bitVector = getExpectedBitVector(expectedBits, state, bitValue);
bitVector = expectedBits(state, bitValue, :)(:);
endfunction
function bitVector = getInputBitVector(symbols, idx);
if (columns(symbols) >= idx)
bitVector = symbols(:,idx);
else
bitVector = [0 ; 0; 0];
end
endfunction------------------------------------------------------------------------------
Stay on top of everything new and different, both inside and
around Java (TM) technology - register by April 22, and save
$200 on the JavaOne (SM) conference, June 2-5, 2009, San Francisco.
300 plus technical and hands-on sessions. Register today.
Use priority code J9JMT32. http://p.sf.net/sfu/p
_______________________________________________
Octave-dev mailing list
[email protected]
https://lists.sourceforge.net/lists/listinfo/octave-dev
