va_demod_zhaduitest.asv

来自「短波信道抗多音干扰的性能分析及其仿真」· ASV 代码 · 共 137 行

% function pb=va_demod(Eb_to_Nj_in_dB,rho_in_dB,BPH,number_of_states,Q)
% this function simulates the viterbi demodulation of
% the DFH system

% rho_in_dB=5;
% BPH=2;
% number_of_states=16;
Eb_to_Nj_in_dB=17;
rho_in_dB=13.35;
BPH=2;
number_of_states=16;
Q=8;
N=1000;
fanout=2^BPH;
rho=10^(rho_in_dB/10);
Eb_to_Nj=10^(Eb_to_Nj_in_dB/10);
L=floor(log(number_of_states)/log(fanout));
times=100;
num_of_err=zeros(

source=[randint(1,N,fanout),zeros(1,L)];
dsource=zeros(1,N*BPH);
if(BPH~=1)
    for i=1:N
        dsource((i-1)*BPH+1:i*BPH)=deci2change(source(i),BPH,2);
    end
else
    dsource=source(1:N);
end

nextstate=zeros(number_of_states,fanout); 
output=zeros(number_of_states,fanout);
input=zeros(number_of_states,number_of_states);
number_of_out=number_of_states*fanout;
for i=0:number_of_states-1
    for j=0:fanout-1
        [next_state,out_put]=G_func1(i,j,L,fanout);
        nextstate(i+1,j+1)=next_state;
        output(i+1,j+1)=out_put;
        input(i+1,next_state+1)=j;
    end
end 

depth_of_trellis=length(source);
Eb=1;
Es=Eb*BPH;
thyta=2*pi*rand;
Ej0=(Es*number_of_out)/(BPH*Q*Eb_to_Nj);
sgma=sqrt(Eb/(BPH*2*rho));
demod_input=zeros(number_of_out,depth_of_trellis);
f=zeros(1,depth_of_trellis);
D=0;
for i=1:depth_of_trellis % i表示网格图的时间走势
    f(i)=output(D+1,source(i)+1); % f(i)是i时刻的分支转移输出,即i时刻的跳频频率号
%   J=randint(1,Q,number_of_out); % J 矩阵中存放Q个干扰音所在的频率号
    J=gen_multijammer(Q,number_of_out);% J 矩阵中存放Q个干扰音所在的频率号,干扰音所在频率号范围也是[0,number_of_out-1]而不是[1,number_of_out]
    for j=0:number_of_out-1
        if (j==f(i))
            rc(j+1)=sqrt(Es)+sgma*randn;
            rs(j+1)=sgma*randn;
        else
            rc(j+1)=sgma*randn;
            rs(j+1)=sgma*randn;
        end
    end
    for k=1:Q
        for j=0:number_of_out-1
            if (j==J(k))
                rc(j+1)=rc(j+1)+sqrt(Ej0)*cos(thyta);
                rs(j+1)=rs(j+1)+sqrt(Ej0)*sin(thyta);
            end
        end
    end
    for j=0:number_of_out-1
        demod_input(j+1,i)=sqrt(rc(j+1)^2+rs(j+1)^2);
    end
    D=nextstate(D+1,source(i)+1);
end

state_metric=zeros(number_of_states,2);
survivor_state=zeros(number_of_states,depth_of_trellis+1);
for i=1:depth_of_trellis-L
    flag=zeros(1,number_of_states);
    if i<=L+1
        step=2^((L+1-i)*BPH);
    else
        step=1;
    end
    for j=0:step:number_of_states-1
        for m=0:fanout-1
            branch_metric=demod_input(output(j+1,m+1)+1,i);
            if((state_metric(nextstate(j+1,m+1)+1,2)<state_metric(j+1,1)...
                    +branch_metric)|flag(nextstate(j+1,m+1)+1)==0)
                state_metric(nextstate(j+1,m+1)+1,2)=state_metric(j+1,1)+branch_metric;
                survivor_state(nextstate(j+1,m+1)+1,i+1)=j;
                flag(nextstate(j+1,m+1)+1)=1;
            end
        end
    end
    state_metric=state_metric(:,2:-1:1);
end
for i=depth_of_trellis-L+1:depth_of_trellis
    flag=zeros(1,number_of_states);
    last_stop=number_of_states/(2^((i-depth_of_trellis+L-1)*BPH));
    for j=0:last_stop-1
        branch_metric=demod_input(output(j+1,m+1)+1,i);
        if((state_metric(nextstate(j+1,1)+1,2)<state_metric(j+1,1)...
            +branch_metric)|flag(nextstate(j+1,1)+1)==0)
            state_metric(nextstate(j+1,1)+1,2)=state_metric(j+1,1)+branch_metric;
            survivor_state(nextstate(j+1,1)+1,i+1)=j;
            flag(nextstate(j+1,1)+1)=1;
        end
    end
    state_metric=state_metric(:,2:-1:1);
end
state_sequence=zeros(1,depth_of_trellis+1);
for i=1:depth_of_trellis
    state_sequence(1,depth_of_trellis-i+1)=survivor_state((state_sequence(1,depth_of_trellis+2-i)...
        +1),depth_of_trellis-i+2);
end
decoder_output=zeros(1,BPH*(depth_of_trellis-L));
for i=1:depth_of_trellis-L
    dec_output_deci=input(state_sequence(1,i)+1,state_sequence(1,i+1)+1);
    if(BPH~=1)
        dec_output_bin=deci2change(dec_output_deci,BPH,2);
        decoder_output((i-1)*BPH+1:i*BPH)=dec_output_bin;
    else
        decoder_output(1,i)=dec_output_deci;
    end
end

for i=1:N*BPH
    if(dsource(i)~=decoder_output(i))
       num_of_err(rep)=num_of_err(rep)+1;
   end
end
pb=num_of_err/(N*BPH);