www.pudn.com > communicationmatlab.rar > SVITERBI.M
function [sys, x0, str, ts] = sviterbi(t, x, u, flag, tran_func, leng, tran_prob, plot_flag);
%SVITERBI SIMULINK file for convolution decoding using viterbi algorithm.
% This file is designed to be used in a SIMULINK S-Function block.
% The function requires the system inputs
% code_param = [N, K, M, num_state, num_std_sta]; % prepared by simviter
% tran_func = [A, B; C, D]; % prepared by simviter
% leng % memory length
% tran_prob % transfer probability
% % when it is not a three row matrix, take the code as hard decision one.
% plot_flag % should plot or not.
% % when it is not a positive integer, don't plot.
% % when it is a positive integer, keep the plot have the time length
% % as required.
% This function keeps a number of discrete-time states:
% figure number. -Inf: to be opened; 0 not need to plot; positive: handle
% figure posit. Current figure position.
% trace_num. current trace number.
% trace_flag. flag to indicate that computation is not under leng.
% stater. variable in keeping the very beinging of state.
% trace. a LENG-by-NUM_STD_STA matrix storage the traces.
% solut. a LENG-by-NUM_STD_STA matrix storage the possible msg.
% expense. a 2-by-NUM_STD_STA matrix storage the expense.
%
% Wes Wang
% Copyright (c) 1995-96 by The MathWorks, Inc.
if (flag == 2) % refresh discrete-time states
% t2 = clock;
% the major processing routine.
if u(length(u)) < .2
% in the case of no signal, no processing.
return;
end;
% otherwise, there is a signal, have to process.
u = u(:)';
% initial parameters.
[A, B, C, D, N, K, M] = gen2abcd(tran_func);
num_state = M;
n_std_sta = 2^M;
PowPowM = n_std_sta^2;
K2 = 2^K;
inp_pre = de2bi([0:K2-1]', K);
cur_sta_pre = de2bi([0:n_std_sta-1], M);
%%%%%begin plot related %%%%%
if x(1) < 0
% initialize the figure.
[sl_name, block] = get_param;
[n_b, m_b]= size(block);
if n_b < 1
error('Cannot delete block while simulating')
elseif n_b > 1
error('Something wrong in get_param, You don''t have the current SIMULINK')
end;
% findout if the graphics window exist
ind = find(sl_name == setstr(10));
dash = '_';
sl_name(ind)=dash(ones(size(ind)));
% in the case of the figure exists
Figures = get(0, 'Child');
new_figure = 1;
i = 1;
while ((new_figure) & (i <= length(Figures)))
if strcmp(get(Figures(i), 'Type'), 'figure')
if strcmp(get(Figures(i), 'Name'), sl_name)
h_fig = Figures(i);
handles = get(h_fig,'UserData');
new_figure = 0;
h_axis = handles(1);
set(0,'CurrentF',h_fig);
set(h_axis,'Xlim', [0, plot_flag-1], 'Ylim', [0, n_std_sta-1]);
delete(get(h_axis,'child'));
handles = handles(1);
end;
end;
i = i + 1;
end;
if new_figure
h_fig = figure;
set(h_fig, 'Name', sl_name);
handles = axes('Xlim',[0, plot_flag-1], 'Ylim',[0, n_std_sta-1]);
end;
handles = get(h_fig, 'Child');
set(handles, 'Visible', 'Off');
set(handles(1),'NextPlot','add');
set(h_fig,'NextPlot','add', 'Clipping', 'off');
handles(2) = plot(0,0,'y-','Erasemode','normal');
handles(3) = plot(0,0,'r-','Erasemode','normal','LineWidth',2);
handles(4) = plot(-10,-10, 'go', 'Erasemode','normal');
for i = 1 : n_std_sta
handles(i+4) = text(0,0,'0');
set(handles(i+4), 'Color',[0,1 1]);
end;
for i = 0 : n_std_sta - 1
handles(5+n_std_sta+i) = text(-plot_flag/8, i, ['State ', num2str(bi2de(fliplr(de2bi(i, num_state))))]);
end;
handles(5+2*n_std_sta) = text(plot_flag/2.1, - n_std_sta/15, 'Time');
set(h_fig, 'UserData',handles);
set(h_axis,'NextPlot','new');
set(h_fig,'NextPlot','new');
sys(1) = h_fig;
x(1) = h_fig;
end;% if x(1) < 0
plot_flag_test = get(0,'child');
if isempty(plot_flag_test)
plot_flag_test = 0;
elseif isempty(find(plot_flag_test==x(1)))
plot_flag_test = 0;
else
plot_flag_test = 1;
handles = get(x(1), 'UserData');
if length(handles) ~= (5+2*n_std_sta)
plot_flag_test = 0;
end;
end;
%%%%% end plot related %%%%%
[n_tran_prob, m_tran_prob] = size(tran_prob);
if n_tran_prob == 3
expen_flag = 1; % expense flag == 0; BSC code. == 1, with TRAN_PROB.
expen_work = zeros(1, N); % work space.
tran_prob(2:3, :) = log10(tran_prob(2:3, :));
else
expen_flag = 0;
end;
% veraible to record the weight trace back to the first state.
% the column number is the state number - 1.
% the first line is the previous. The second line is the current.
starter = x(5);
solu = zeros(leng, PowPowM);
expense = solu;
code = zeros(leng, N);
loc_tmp = 6;
expense(:) = x(loc_tmp : leng*PowPowM + loc_tmp - 1);
solu(:) = x(loc_tmp + leng * PowPowM : 2 * leng * PowPowM + loc_tmp - 1);
code(:) = x(loc_tmp + 2 * leng * PowPowM : loc_tmp + 2 * leng * PowPowM -1 + leng * N);
%%%%% begin plot related %%%%%
fig_position = x(2) + 1;
if (x(1) > 0) & (rem(fig_position-leng, plot_flag - leng) == 0) & (fig_position >= plot_flag) & plot_flag_test
set(handles(1), 'XLim', [fig_position-leng, fig_position + plot_flag - leng]);
plot_curv_x = get(handles(2), 'XData');
plot_curv_y = get(handles(2), 'YData');
index = find(plot_curv_x == -(fig_position - leng));
% index = find(plot_curv_x <= (fig_position - leng -1));
leng_x = length(plot_curv_x);
if isempty(index)
% no move
else
% index = index(length(index));
plot_curv_x = plot_curv_x(index(1)+2: leng_x);
plot_curv_y = plot_curv_y(index(1)+2: leng_x);
end;
set(handles(2), 'XData', plot_curv_x, 'YData', plot_curv_y);
for i = 0:n_std_sta - 1
tmp = get(handles(5+n_std_sta+i), 'Position');
tmp(1) = fig_position-leng-plot_flag/8;
set(handles(5+n_std_sta+i), 'Position', tmp);
end;
tmp = get(handles(5+2*n_std_sta), 'Position');
tmp(1) = fig_position-leng+plot_flag/2.1;
set(handles(5+2*n_std_sta), 'Position', tmp);
plot_curv_x = get(handles(4), 'XData');
plot_curv_y = get(handles(4), 'YData');
index = find(plot_curv_x <= fig_position-leng);
index = max([index(length(index)) 0]) + 1;
leng_x = length(plot_curv_x);
plot_curv_x = plot_curv_x(index : leng_x);
plot_curv_y = plot_curv_y(index : leng_x);
set(handles(4), 'XData', plot_curv_x, 'YData', plot_curv_y);
% drawnow;
end;
trace_num = x(3) + 1;
trace_flag = x(4);
code(trace_num, :) = u(1:length(u)-1);
if (trace_flag == 0) & (trace_num == leng)
trace_flag = 1;
end;
%%%%% end plot related %%%%%
trace_pre = rem(trace_num - 2 + leng, leng) + 1;
% fill up trace and solut
if (trace_flag == 0) & (trace_num == 1)
pre_state = starter + 1;
else
pre_state = [];
for j2 = 1 : n_std_sta
if max(~isnan(expense(trace_pre, [1-n_std_sta : 0] + j2*n_std_sta)))
pre_state = [pre_state, j2];
end;
end;
end;
%%%%% begin plot related %%%%%
if (x(1) > 0) & plot_flag_test
plot_curv_x = [get(handles(2), 'XData'), NaN, -fig_position, NaN];
plot_curv_y = [get(handles(2), 'YData'), NaN, -1000, NaN];
dot_x = get(handles(4), 'XData');
dot_y = get(handles(4), 'YData');
end;
%%%%% end plot related %%%%%
expense(trace_num, :) = expense(trace_num,:) * NaN;
if expen_flag
for j = 1 : length(pre_state)
jj = pre_state(j) - 1; % index number - 1 is the state.
cur_sta = cur_sta_pre(pre_state(j),:)';
indx_j = (pre_state(j) - 1) * n_std_sta;
for num_N = 1 : N
expen_work(num_N) = max(find(tran_prob(1,:) <= code(trace_num, num_N)));
end;
for num_k = 1 : K2
inp = inp_pre(num_k, :)';
if isempty(C)
% tran_indx = pre_state(j) + (num_k -1) * K2;
tran_indx = pre_state(j) + (num_k -1) * n_std_sta;
nex_sta = A(tran_indx, :)';
out = B(tran_indx, :)';
else
out = rem(C * cur_sta + D * inp,2);
nex_sta = rem(A * cur_sta + B * inp, 2);
end;
nex_sta_de = bi2de(nex_sta') + 1;
% find the expense by the transfer probability
expen_0 = find(out' <= 0.5);
expen_1 = find(out' > 0.5);
loca_exp = sum([tran_prob(2,expen_work(expen_0)) 0])...
+sum([tran_prob(3,expen_work(expen_1)) 0]);
tmp = (nex_sta_de-1)*n_std_sta + pre_state(j);
if isnan(expense(trace_num, tmp))
expense(trace_num, tmp) = loca_exp;
solu(trace_num, nex_sta_de + indx_j) = num_k;
elseif expense(trace_num, tmp) < loca_exp
expense(trace_num, tmp) = loca_exp;
solu(trace_num, nex_sta_de + indx_j) = num_k;
end;
end;
end;
else %expen_flag = 0
% t0 = clock;
for j = 1 : length(pre_state)
jj = pre_state(j) - 1; % index number - 1 is the state.
cur_sta = cur_sta_pre(pre_state(j),:)';
indx_j = (pre_state(j) - 1) * n_std_sta;
for num_k = 1 : K2
inp = inp_pre(num_k, :)';
if isempty(C)
% tran_indx = pre_state(j) + (num_k -1) * K2;
tran_indx = pre_state(j) + (num_k -1) * n_std_sta;
nex_sta = A(tran_indx, :)';
out = B(tran_indx, :)';
else
out = rem(C * cur_sta + D * inp,2);
nex_sta = rem(A * cur_sta + B * inp, 2);
end;
nex_sta_de = bi2de(nex_sta') + 1;
loca_exp = sum(rem(code(trace_num,:) + out', 2));
tmp = (nex_sta_de-1)*n_std_sta + pre_state(j);
if isnan(expense(trace_num, tmp))
expense(trace_num, tmp) = loca_exp;
solu(trace_num, nex_sta_de + indx_j) = num_k;
elseif expense(trace_num, tmp) > loca_exp
expense(trace_num, tmp) = loca_exp;
solu(trace_num, nex_sta_de + indx_j) = num_k;
end;
end;
end;
% disp(etime(clock, t0))
end;
aft_state = [];
for j2 = 1 : n_std_sta
if max(~isnan(expense(trace_num, [1-n_std_sta : 0] + j2*n_std_sta)))
aft_state = [aft_state, j2];
end
end;
%%%%% begin plot related %%%%%
if (x(1) > 0) & plot_flag_test
% plot only; no computation function here.
flipped = zeros(1, n_std_sta);
for i = 1 : n_std_sta
flipped(i) = bi2de(fliplr(de2bi(i-1, M)));
end;
for j2 = 1 : length(pre_state)
jj = pre_state(j2) - 1;
for j = 1 : n_std_sta
tmp = expense(trace_num, (j-1)*n_std_sta+1: j*n_std_sta);
if expen_flag
expen_rec(j) = max(tmp(find(~isnan(tmp))));
else
expen_rec(j) = min(tmp(find(~isnan(tmp))));
end;
for j2 = 1 : length(tmp)
if ~isnan(tmp(j2))
plot_curv_x = [plot_curv_x, fig_position-1+.1, fig_position, NaN];
plot_curv_y = [plot_curv_y, ...
flipped(j2), ...
flipped(j), NaN];
end;
end;
end;
end;
for j2 = 1 : length(aft_state)
jj = aft_state(j2) - 1;
set(handles(jj+5), 'Position', [fig_position, bi2de(fliplr(de2bi(jj,num_state))), 0],...
'String', num2str(expen_rec(jj+1)));
tmp = get(handles(jj+5), 'Position');
dot_x = [dot_x, tmp(1)];
dot_y = [dot_y, tmp(2)];
end;
set(handles(2), 'XData', plot_curv_x, 'YData', plot_curv_y);
set(handles(4), 'XData', dot_x, 'YData', dot_y);
if ~trace_flag
drawnow;
end;
end;
%%%%% end plot related %%%%%
% decision making.
% t0 = clock;
if trace_flag
sol = expense(trace_num,:);
trace_eli = rem(trace_num, leng) + 1;
% strike out the unnecessary.
for j_k = 1 : leng - 2
j_pre = rem(trace_num - j_k - 1 + leng, leng) + 1;
sol = vitshort(expense(j_pre, :), sol, n_std_sta, expen_flag);
end;
tmp = (ones(n_std_sta,1) * expense(trace_eli, [starter+1:n_std_sta:PowPowM]))';
sol = sol + tmp(:)';
if expen_flag
loc_exp = max(sol(find(~isnan(sol))));
else
loc_exp = min(sol(find(~isnan(sol))));
end
dec = find(sol == loc_exp);
dec = dec(1);
dec = rem((dec - 1), n_std_sta);
num_k = solu(trace_eli, starter*n_std_sta+dec+1);
inp = inp_pre(num_k, :)';
if isempty(C)
% tran_indx = pre_state(j) + (num_k -1) * K2;
tran_indx = starter + 1 + (num_k -1) * n_std_sta;
out = B(tran_indx, :)';
else
cur_sta = cur_sta_pre(starter+1, :)';
out = rem(C * cur_sta + D * inp,2);
end;
if expen_flag
% find the expense by the transfer probability
expen_0 = find(out' <= 0.5);
expen_1 = find(out' > 0.5);
loc_exp = sum([tran_prob(2,expen_work(expen_0)) 0])...
+sum([tran_prob(3,expen_work(expen_1)) 0]);
else
loc_exp = sum(rem(code(trace_eli, :) + out', 2));
end
%%%%%begin plot related %%%%%
if plot_flag & plot_flag_test
plot_curv_x = [get(handles(3), 'XData'), NaN];
plot_curv_y = [get(handles(3), 'YData'), NaN];
if length(plot_curv_x) > 3 * plot_flag
leng_x = length(plot_curv_x);
limi = leng_x - 3 * plot_flag + 1;
plot_curv_x = plot_curv_x(limi:leng_x);
plot_curv_y = plot_curv_y(limi:leng_x);
end;
set(handles(3), 'XData', [plot_curv_x, fig_position-leng+.1, fig_position-leng+1],...
'YData', [plot_curv_y, bi2de(fliplr(de2bi(starter, num_state))),...
bi2de(fliplr(de2bi(dec, num_state)))]);
drawnow
end;
%%%%% end plot related %%%%%
starter = dec;
output = [inp(:); loca_exp];
else %(if trace_flag)
output = zeros(K+1, 1);
starter = 0;
end; %(if trace_flag)
% disp(['decision etime ', num2str(etime(clock, t0))])
trace_num = rem(trace_num, leng);
sys = [x(1); fig_position; trace_num; trace_flag; starter; expense(:); solu(:); code(:); output(:)];
% disp(['total etime ', num2str(etime(clock, t2))])
elseif flag == 3 % output
if u(length(u)) > 0.2
K = tran_func(2, size(tran_func, 2));
len_x = length(x);
sys = x(len_x - K: len_x);
end;
elseif flag == 0
[A, B, C, D, N, K, M] = gen2abcd(tran_func);
num_state = M;
n_std_sta = 2^M;
PowPowM = n_std_sta^2;
[n_tran_prob, m_tran_prob] = size(tran_prob);
if n_tran_prob == 3
expen_flag = 1; % expense flag == 0; BSC code. == 1, with TRAN_PROB.
else
expen_flag = 0;
end;
% veraible to record the weight trace back to the first state.
% the column number is the state number - 1.
% the first line is the previous. The second line is the current.
expense = ones(leng, PowPowM) * NaN;
expense(leng, [1:n_std_sta]) = zeros(1, n_std_sta);
solu = zeros(leng, PowPowM);
% the column number is the state + 1.
% the contents is the state it transfered from.
% the starter is always single number.
starter = 0;
% the solution in each of the transfer as above.
% ith row is the transition input (msg) from (i-1)th row of trace to ith row
% of trace. When i = 1, the transfer is from starter to trace.
if plot_flag > 0
x0 = -Inf;
else
x0 = 0;
end;
code = zeros(leng, N);
y = zeros(K+1, 1);
% add output storage.
x0 = [x0; 0; 0; 0; starter; expense(:); solu(:); code(:); y];
% | | | | | | | \- output
% \ \ \ \ \ \ \ decode
% \ \ \ \ \start \-weight, expense
% \ \ \ \-trace_flag
% \ \ \-trace_num
% \ \-fig_position
% \-figure handle
sys = [0; length(x0); K+1; N+1; 0; 0; 1];
ts = [-1, 0];
elseif flag == 9
if x(1) > 0
plot_flag_test = get(0,'child');
if isempty(plot_flag_test)
plot_flag_test = 0;
elseif isempty(find(plot_flag_test==x(1)))
plot_flag_test = 0;
else
plot_flag_test = 1;
handles = get(x(1), 'UserData');
if length(handles) < (5+2*n_std_sta)
plot_flag_test = 0;
end;
end;
if plot_flag_test
plot_curv_x = get(handles(4), 'XData');
plot_curv_y = get(handles(4), 'YData');
set(handles(4), 'Erasemode','Normal',...
'XData', plot_curv_x,...
'YData', plot_curv_y);
drawnow
end;
end;
sys = [];
else
sys = [];
end;