matlab中的turbo码,zz TURBO码MATLAB仿真（二）-CSDN博客

function L_all = sova(rec_s, g, L_a, ind_dec)

% This function implememts Soft Output Viterbi Algorithm in trace

back mode

% Input:

% rec_s: scaled received bits. rec_s(k) = 0.5 * L_c(k) * y(k)

% L_c = 4 * a * Es/No, reliability value of the channel

% y: received bits

% g: encoder generator matrix in binary form,

g(1,:) for feedback, g(2,:) for feedforward

% L_a: a priori information about the info. bits. Extrinsic info.

from the previous

% component decoder

% ind_dec: index of the component decoder.

% =1: component decoder 1; The trellis is terminated to all zero

state

% =2: component decoder 2; The trellis is not perfectly

terminated.

% Output:

% L_all: log ( P(x=1|y) ) / ( P(x=-1|y) )

% Copyright: Yufei Wu, Nov. 1998

% MPRG lab, Virginia Tech

% for academic use only

% Frame size, info. + tail bits

L_total = length(L_a);

[n,K] = size(g);

m = K - 1;

nstates = 2^m;

Infty = 1e10;

% SOVA window size. Make decision after 'delta' delay. Decide

bit k when received bits

% for bit (k+delta) are processed. Trace back from (k+delta) to

delta =

30;

% Set up the trellis defined by g.

[next_out, next_state, last_out, last_state] = trellis(g);

% Initialize path metrics to -Infty

for t=1:L_total+1

for state=1:nstates

path_metric(state,t) = -Infty;

end

% Trace forward to compute all the path metrics

path_metric(1,1) = 0;

for t=1:L_total

y = rec_s(2*t-1:2*t);

for state=1:nstates

sym0 = last_out(state,1:2);

sym1 = last_out(state,3:4);

state0 = last_state(state,1);

state1 = last_state(state,2);

Mk0 = y*sym0' - L_a(t)/2 + path_metric(state0,t);

Mk1 = y*sym1' + L_a(t)/2 + path_metric(state1,t);

if Mk0>Mk1

path_metric(state,t+1)=Mk0;

Mdiff(state,t+1) = Mk0 - Mk1;

prev_bit(state, t+1) = 0;

else

path_metric(state,t+1)=Mk1;

Mdiff(state,t+1) = Mk1 - Mk0;

prev_bit(state,t+1) = 1;

end

% For decoder 1, trace back from all zero state,

% for decoder two, trace back from the most likely state

if ind_dec == 1

mlstate(L_total+1) = 1;

else

mlstate(L_total+1) = find(

path_metric(:,L_total+1)==max(path_metric(:,L_total+1)) );

end

% Trace back to get the estimated bits, and the most likely

path

for t=L_total:-1:1

est(t) =

prev_bit(mlstate(t+1),t+1);

mlstate(t) =

last_state(mlstate(t+1), est(t)+1);

end

% Find the minimum delta that corresponds to a compitition

path with different info. bit

estimation. % Give the soft output

for t=1:L_total

llr = Infty;

for i=0:delta

if t+i

bit = 1-est(t+i);

temp_state = last_state(mlstate(t+i+1), bit+1);

for j=i-1:-1:0

bit = prev_bit(temp_state,t+j+1);

temp_state = last_state(temp_state, bit+1);

end

if bit~=est(t)

llr = min( llr,Mdiff(mlstate(t+i+1), t+i+1) );

end

L_all(t) = (2*est(t) - 1) *

llr;

end ----------------------------------------------------------

function [next_out, next_state, last_out, last_state] =

trellis(g)

% copyright Nov. 1998 Yufei Wu

% MPRG lab, Virginia Tech

% for academic use only

% set up the trellis given code generator g

% g given in binary matrix form. e.g. g = [ 1 1 1; 1 0 1 ];

% next_out(i,1:2): trellis next_out (systematic bit; parity

bit) when input = 0, state = i; next_out(i,j) = -1 or 1

% next_out(i,3:4): trellis next_out (systematic

bit; parity bit) when input = 1, state = i;

% next_state(i,1): next state when input = 0, state = i;

next_state(i,i) = 1,...2^m

% next_state(i,2): next state when input = 1, state = i;

% last_out(i,1:2): trellis last_out (systematic bit; parity bit)

when input = 0, state = i; last_out(i,j) = -1 or 1

% last_out(i,3:4): trellis last_out (systematic

bit; parity bit) when input = 1, state = i;

% last_state(i,1): previous state that comes to state i when info.

bit = 0;

% last_state(i,2): previous state that comes to state i when info.

bit = 1;

[n,K] = size(g);

m = K - 1;

max_state = 2^m;

% set up next_out and next_state matrices for systematic

code

for state=1:max_state

state_vector = bin_state(

state-1, m );

% when receive a 0

d_k = 0;

a_k = rem( g(1,:)*[0

state_vector]', 2 );

[out_0, state_0] =

encode_bit(g, a_k, state_vector);

out_0(1) = 0;

% when receive a 1

d_k = 1;

a_k = rem( g(1,:)*[1

state_vector]', 2 );

[out_1, state_1] =

encode_bit(g, a_k, state_vector);

out_1(1) = 1;

next_out(state,:) = 2*[out_0

out_1]-1;

next_state(state,:) =

[(int_state(state_0)+1) (int_state(state_1)+1)];

end

% find out which two previous states can come to present

state

last_state = zeros(max_state,2);

for bit=0:1

for state=1:max_state

last_state(next_state(state,bit+1), bit+1)=state;

last_out(next_state(state, bit+1), bit*2+1:bit*2+2) ...

= next_out(state, bit*2+1:bit*2+2);

end

-------------------------------------------------------------

% This script simulates the classical turbo encoding-decoding

system.

% It simulates parallel concatenated convolutional codes.

% Two component rate 1/2 RSC (Recursive Systematic Convolutional)

component encoders are assumed.

% First encoder is terminated with tails bits. (Info + tail) bits

are scrambled and passed to

% the second encoder, while second encoder is left open without

tail bits of itself.

% Random information bits are modulated into +1/-1, and transmitted

through a AWGN channel.

% Interleavers are randomly generated for each frame.

% Log-MAP algorithm without quantization or approximation is

used.

% By making use of ln(e^x+e^y) = max(x,y) +

ln(1+e^(-abs(x-y))),

% the Log-MAP can be simplified with a look-up table for the

correction function.

% If use approximation ln(e^x+e^y) = max(x,y), it becomes

MAX-Log-MAP.

% Copyright Nov 1998, Yufei Wu

% MPRG lab, Virginia Tech.

% for academic use only

clear all

% Write display messages to a text file

diary turbo_logmap.txt

% Choose decoding algorithm

dec_alg = input(' Please enter the decoding algorithm. (0:Log-MAP,

1:SOVA) default

0 ');

if isempty(dec_alg)

dec_alg = 0;

end

% Frame size

L_total = input(' Please enter the frame size (= info + tail,

default: 400) ');

if isempty(L_total)

L_total =

400; % infomation bits plus tail bits

end

% Code generator

g = input(' Please enter code generator: ( default: g = [1 1 1; 1 0

1 ]

) ');

if isempty(g)

g = [ 1 1 1;

1 0 1 ];

end

%g = [1 1 0 1; 1 1 1 1];

%g = [1 1 1 1 1; 1 0 0 0 1];

[n,K] = size(g);

m = K - 1;

nstates = 2^m;

%puncture = 0, puncturing into rate 1/2;

%puncture = 1, no puncturing

puncture = input(' Please choose punctured / unpunctured (0/1):

default

0 ');

if isempty(puncture)

puncture =

end

% Code rate

rate = 1/(2+puncture);

% Fading amplitude; a=1 in AWGN channel

a = 1;

% Number of iterations

niter = input(' Please enter number of iterations for each frame:

default

5 ');

if isempty(niter)

niter = 5;

end % Number of frame errors to count as a stop criterior

ferrlim = input(' Please enter number of frame errors to terminate:

default

15 ');

if isempty(ferrlim)

ferrlim = 15;

end

EbN0db = input(' Please enter Eb/N0 in dB : default

[2.0] ');

if isempty(EbN0db)

EbN0db = [2.0];

end

fprintf('\n\n----------------------------------------------------\n');

if dec_alg == 0

fprintf(' === Log-MAP decoder

=== \n');

else

fprintf(' === SOVA decoder ===

\n');

end

fprintf(' Frame size = %6d\n',L_total);

fprintf(' code generator: \n');

for i = 1:n

for j =

1:K

fprintf( '%6d', g(i,j));

end

fprintf('\n');

end if puncture==0

fprintf(' Punctured, code rate

= 1/2 \n');

else

fprintf(' Unpunctured, code

rate = 1/3 \n');

end

fprintf(' iteration number = %6d\n',

niter);

fprintf(' terminate frame errors = %6d\n', ferrlim);

fprintf(' Eb / N0 (dB) = ');

for i = 1:length(EbN0db)

fprintf('%10.2f',EbN0db(i));

end

fprintf('\n----------------------------------------------------\n\n');

fprintf('+ + + + Please be patient. Wait a while to get the result.

+ + + +\n');

for nEN = 1:length(EbN0db)

en =

10^(EbN0db(nEN)/10); % convert Eb/N0 from unit db to normal numbers

L_c = 4*a*en*rate;

% reliability value of the channel

sigma = 1/sqrt(2*rate*en);

% standard deviation of AWGN noise

% Clear bit error counter and frame error counter

errs(nEN,1:niter) =

zeros(1,niter);

nferr(nEN,1:niter) =

zeros(1,niter);

nframe =

0; % clear

counter of transmitted frames

while nferr(nEN,

niter)

nframe = nframe +

1; x = round(rand(1,

L_total-m)); % info. bits

[temp, alpha] =

sort(rand(1,L_total)); % random interleaver mapping

en_output = encoderm( x, g, alpha, puncture ) ; % encoder output

(+1/-1)

r = en_output+sigma*randn(1,L_total*(2+puncture)); % received

bits

yk = demultiplex(r,alpha,puncture); % demultiplex to get input for

decoder 1 and 2

% Scale the received

bits rec_s = 0.5*L_c*yk;

% Initialize extrinsic

information L_e(1:L_total) = zeros(1,L_total);

for iter = 1:niter

% Decoder one

L_a(alpha) = L_e; % a priori info.

if dec_alg == 0

L_all = logmapo(rec_s(1,:), g, L_a, 1); %

complete info.

else L_all = sova0(rec_s(1,:), g, L_a, 1); % complete

info.

end L_e = L_all - 2*rec_s(1,1:2:2*L_total) - L_a; %

extrinsic info.

% Decoder

two L_a = L_e(alpha); % a priori info.

if dec_alg == 0

L_all = logmapo(rec_s(2,:), g, L_a, 2); %

complete info. else

L_all = sova0(rec_s(2,:), g, L_a, 2); % complete

info.

end

L_e = L_all - 2*rec_s(2,1:2:2*L_total) - L_a; %

extrinsic info.

% Estimate the info.

bits xhat(alpha) = (sign(L_all)+1)/2;

% Number of bit errors in current iteration

err(iter) = length(find(xhat(1:L_total-m)~=x));

% Count frame errors for the current iteration

if err(iter)>0

nferr(nEN,iter) = nferr(nEN,iter)+1;

end end %iter

% Total number of bit errors for all iterations

errs(nEN,1:niter) = errs(nEN,1:niter) + err(1:niter);

if rem(nframe,3)==0 | nferr(nEN, niter)==ferrlim

% Bit error rate

ber(nEN,1:niter) = errs(nEN,1:niter)/nframe/(L_total-m);

% Frame error rate

fer(nEN,1:niter) = nferr(nEN,1:niter)/nframe;

% Display intermediate results in

process fprintf('************** Eb/N0 = %5.2f db **************\n',

EbN0db(nEN));

fprintf('Frame size = %d, rate 1/%d. \n', L_total,

2+puncture);

fprintf('%d frames transmitted, %d frames in error.\n', nframe,

nferr(nEN, niter));

fprintf('Bit Error Rate (from iteration 1 to iteration %d):\n',

niter);

for i=1:niter

fprintf('%8.4e ', ber(nEN,i));

end

fprintf('\n');

fprintf('Frame Error Rate (from iteration 1 to iteration %d):\n',

niter);

for i=1:niter

fprintf('%8.4e ', fer(nEN,i));

end

fprintf('\n');

fprintf('***********************************************\n\n');

% Save intermediate results

save turbo_sys_demo EbN0db ber fer

end

end %while

end %nEN

diary off

-----------------------------------------------------