va999.m

是個能夠做附力葉分解動作相當有用的一套程式碼希望能有所幫助

clc;
clear all;
close all;
op_sequence=[0.0241 -0.3733 1.8132 0.6410 -0.1521 0.7423 -0.4910 0.9656 1.2019 0.3474 -0.621 -0.83];
h=[0.3 0.8 -0.1];

% %%%%%%%%%%%%%%%%% I/P & register state %%%%%%%%%%%%%%%%%%%%


ip=[-1 1]';
seq1=[-1 -1];
seq2=[1 1];
seq_up=[];
seq_down=[];
for i=1:(length(h)-1)
    seq_up=[seq1' ip];
    seq_down=[seq2' ip];
    seq1=[seq1 seq1];
    seq2=[seq2 seq2];
    ip_register_state=[seq_up;seq_down];
    ip=ip_register_state;
end
ip_register_state
%%%%%%%%%%%%%%%%%%%%%% O/P state %%%%%%%%%%%%%%%%%%%%
metric_seq=[];
op_seq=ip_register_state*h'
k=0;
m=[];
min_path_metric=[];
path_metric=[];
metric_nan_seq=nan(2^length(h),1)
for i=1:length(op_sequence)
    metric_op_seq=(op_seq-op_sequence(i)).^2
    metric_seq=[metric_seq metric_op_seq]
    r=1;
    q=(2^length(h))/(2^i)
    for j=1:(2^length(h))
        if q<1
            q=1;
        end
        metric_nan_seq(r)=metric_op_seq(r);
        r=r+q;
        if r>(2^length(h))
            break
        end
    end
    n=length(op_sequence)-i;
    z=2^n+1;
    if n<length(h)
        for t=1:(2^(length(h)))
            metric_nan_seq(z)=nan;
            z=z+1;
            if z>2^(length(h))
                break;
            end
        end
    end
    path_metric=[path_metric metric_nan_seq];
    for j=1:2:(2^length(h))
        k=k+1;
        metric(k)=min(metric_nan_seq(j),metric_nan_seq(j+1));
    end
    k=0;
    min_path_metric=[min_path_metric metric'];
end
metric_seq
path_metric
min_path_metric
%%%%%%%%%%%%%%%%min_path_metric & path state %%%%%%%%%%%%%%%%%%%
metric_sum=nan((2^(length(h)-1)),2*length(op_sequence))
g=1;
for j=1:length(op_sequence)
    f=1;
    d=1;
    for i=1:(2^(length(h)-1))
        metric_sum(i,g+1)=min(path_metric(f,j),path_metric(f+1,j));
        if metric_sum(i,g+1)==path_metric(f,j)
            metric_sum(i,g)=1;
        else
            metric_sum(i,g)=-1;
        end
        if (isnan(path_metric(f,j))==1)&&(isnan(path_metric(f+1,j))==1)
            metric_sum(i,g+1)=nan;
            metric_sum(i,g)=nan;    
        end
        if j>1
            m=metric_sum(d,g-1)+path_metric(f,j);
            n=metric_sum(d+1,g-1)+path_metric(f+1,j);
            
            if m>n
                metric_sum(i,g+1)=n;
                metric_sum(i,g)=-1;
            else
                metric_sum(i,g+1)=m;
                metric_sum(i,g)=1;
            end 
            if isnan(path_metric(f,j))==1
                metric_sum(i,g+1)=n;
                metric_sum(i,g)=-1;
            end
            if isnan(path_metric(f+1,j))==1
                metric_sum(i,g+1)=m;
                metric_sum(i,g)=1;
            end
            if (isnan(path_metric(f,j))==1)&&(isnan(path_metric(f+1,j))==1)
                metric_sum(i,g+1)=nan;
                metric_sum(i,g)=nan;
            end
         end
        d=d+2;
        if d>(2^(length(h)-1)-1)
            d=1;
        end
        f=f+2;
    end
    
    g=g+2;
    if g>(2*length(op_sequence))
        break
    end
end  
metric_sum
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
min_metric_sum=nan((2^(length(h)-1)),2*length(op_sequence));
v=2*length(op_sequence);
min_metric_sum(:,v)=metric_sum(:,v);
min_metric_sum(:,v-1)=metric_sum(:,v-1);
v=2*length(op_sequence)-2
d=2*length(op_sequence)+1
for i=1:length(op_sequence)
    if v<2
        break
    end
    if d<0
        break
    end
    b=1;
    f=-3;
    t=2;
    w=-2;
    d=d-2;
    if min(min_metric_sum(:,d))>0
        min_metric_sum(:,d)
        for r=1:2^(length(h)-1)
            if isnan(min_metric_sum(r,d))==0
                if r<=2^(length(h)-2)
                    min_metric_sum(b,v)=metric_sum(b,v)
                    min_metric_sum(b,(v-1))=metric_sum(b,(v-1))
                   
                    if b>2^(length(h)-1)
                        break
                    end
                elseif r>2^(length(h)-2)                    
                    min_metric_sum(f,v)=metric_sum(f,v)
                    min_metric_sum(f,(v-1))=metric_sum(f,(v-1))                   
                    if f>2^(length(h)-1)
                        break
                    end
                end

            end
             b=b+2
             f=f+2
        end
    end
    if min(min_metric_sum(:,d))<0
        t=2;
        for r=1:2^(length(h)-1)
            if isnan(min_metric_sum(r,d))==0                 
                if r<=2^(length(h)-2)
                    min_metric_sum(t,v)=metric_sum(t,v)
                    min_metric_sum(t,(v-1))=metric_sum(t,(v-1))
                    
                    if t>2^(length(h)-1)
                        break
                    end
                elseif r>2^(length(h)-2)                    
                    min_metric_sum(w,v)=metric_sum(w,v)
                    min_metric_sum(w,(v-1))=metric_sum(w,(v-1))
                    if w>2^(length(h)-1)
                        break
                    end
                end
            end
            t=t+2
            w=w+2
        end
    end
    v=v-2
end
min_metric_sum
path_matrix=[];
for i=1:length(op_sequence)
    path_matrix=[path_matrix  min_metric_sum(:,(2*i))];
end
path_matrix

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%plot(0:1:L+m, survivor_stage_index, '-kO');
%plot(0:1:L+m, kron(ones(1, L+m+1), [1:num_stage]'), '+');
optimal_path_seq=[1];
optimal_metric_seq=[0];
for i=1:length(op_sequence)
    for j=1:2^(length(h)-1)
        if isnan(path_matrix(j,i))==0
            z=j
            q=path_matrix(j,i)
        end
    end
    optimal_path_seq=[optimal_path_seq z];
    optimal_metric_seq=[optimal_metric_seq q];
end
path_matrix
optimal_path_seq
optimal_metric_seq
clf;
figure(1); hold on;
for i=1:length(optimal_path_seq)
    if optimal_path_seq(i)==1
        r=4
    end
    if optimal_path_seq(i)==2
        r=3
    end
    if optimal_path_seq(i)==3
        r=2
    end
    if optimal_path_seq(i)==4
        r=1
    end
    optimal_path_seq(i)=r;
end
plot(0:1:length(op_sequence), kron(ones(1, length(op_sequence)+1), [1:2^(length(h)-1)]'), '+');
plot(0:1:length(op_sequence), optimal_path_seq, '-kO');
for kk = 0 : 1 : length(op_sequence)
    text(kk-0.1, optimal_path_seq(kk+1)+0.2, num2str(round(optimal_metric_seq(kk+1)*100)/100) );
end
xlabel('Time index');
ylabel('Stage index');
title(['Viterbi Homework  ']);        % show tiltle in vector form
flops;