% 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 onlyclear all% Write display messages to a text filediary 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 sizeL_total = input(' Please enter the frame size (= info + tail, default: 128)   ');if isempty(L_total)   L_total = 128;	 % infomation bits plus tail bitsend% Code generatorg = input(' Please enter code generator: ( default: g = [1 1 1; 1 0 1  ] )      ');if isempty(g)   g = [ 1 0 0 1 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];%g = [ 1 1 1 0 1 ;   1 0 0 1 1];% g = [1 1 0 0 1; 1 0 1 1 1 ];g = [1 1 1 ;    1 0 1];[n,K] = size(g); m = K - 1;nstates = 2^m;%puncture = 0, puncturing into rate 1/2; %puncture = 1, no puncturingpuncture = input(' Please choose punctured / unpunctured (0/1): default 0     ');if isempty(puncture)     puncture = 0;end% Code raterate = 1/(2+puncture);   % Fading amplitude; a=1 in AWGN channela = 1; % Number of iterationsniter = 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 criteriorferrlim = 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];endfprintf('\n\n----------------------------------------------------\n'); if dec_alg == 0   fprintf(' === Log-MAP decoder === \n');else   fprintf(' === SOVA decoder === \n');endfprintf(' 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');endfprintf(' 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));endfprintf('\n----------------------------------------------------\n\n');    fprintf('+ + + + Please be patient. Wait a while to get the result. + + + +\n');%**************************** Preparation part *****************************sr   = 256000.0;                                                    % symbol rateml   = 1;                                                           % number of modulation levelsbr   = sr * ml;                                                     % bit rate%nd   = 100;                                                         % number of symbol%************************** Filter initialization **************************irfn   = 21;                                                        % number of filter tapsIPOINT =  8;                                                        % number of oversamplealfs   =  0.5;                                                      % roll off factor[xh]   = hrollfcoef(irfn,IPOINT,sr,alfs,1);                         % T FILTER FUNCTION[xh2]  = hrollfcoef(irfn,IPOINT,sr,alfs,0);                         % R FILTER FUNCTION%********************** Spreading code initialization **********************user  = 2;                                                          % number of usersseq   = 2;                                                          % 1:M-sequence  2:Gold  3:Orthogonal Goldstage = 3;                                                          % number of stagesptap1 = [1 3];                                                      % position of taps for 1stptap2 = [2 3];                                                      % position of taps for 2ndregi1 = [1 1 1];                                                    % initial value of register for 1stregi2 = [1 1 1];                                                    % initial value of register for 2nd%******************** Generation of the spreading code *********************%switch seq%case 1                                                              % M-sequence%    code = mseq(stage,ptap1,regi1,user);%case 2                                                              % Gold sequence%    m1   = mseq(stage,ptap1,regi1);%    m2   = mseq(stage,ptap2,regi2);%    code = goldseq(m1,m2,user);%case 3                                                              % Orthogonal Gold sequence%    m1   = mseq(stage,ptap1,regi1);%    m2   = mseq(stage,ptap2,regi2);%    code = [goldseq(m1,m2,user),zeros(user,1)];%end%code = code * 2 - 1;%**************************************************************************code=ones(user,10);%code(:,13:2:40)=-1;code(2,1:9:10)=-1;%code=code/40.0;clen = length(code);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*a*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)<ferrlim      nframe = nframe + 1;          x = round(rand(user, L_total-m));    % info. bits                 for i=1:user      [temp, alpha(i,:)] = sort(1:2*L_total);        % random interleaver mapping      en_output(i,:) = encoderm( x(i,:), g, alpha, puncture ) ; % encoder output (+1/-1)      interleave_output(i,:)=en_output(i,alpha(i,:));      end      nd=L_total*(2+puncture);           % code=mseq(stage,ptap1,regi1,user);      ich1=spread(interleave_output,code);      ich2=compoversamp2(ich1,IPOINT);%over sample      ich3=compconv2(ich2,xh); %filter            if user==1          ich4=ich3;      else          ich4=sum(ich3);      end      spow=sum(rot90(ich3.^2 ))/nd;      attn=sqrt(0.5*spow*sr/(br)*10^(-EbN0db(nEN)/10));      ich6=comb2(ich4,attn); %add gaussian AWGN      ich7=compconv2(ich6,xh2);      sampl=irfn*IPOINT+1;      ich8=ich7(:,sampl:IPOINT:IPOINT*nd*clen+sampl-1);      r=despread(ich8,code);      % Initialize extrinsic information            L_e = zeros(user,L_total*2);           L_a = zeros(user,L_total*2);      L_check=zeros(user,L_total*2);      rec_s = zeros(user,L_total*(2+puncture));      L_k = zeros(user,L_total*(2+puncture));      L_all = zeros(user,L_total*2);            rec_s = r/10;   for iter = 1:niter     L_check=multi_dete(rec_s,L_e,sigma);    %xhat=sign(L_check+L_e);% Number of bit errors in current iteration    %     err(2*iter-1) = length(find(xhat(1,1:1:2*(L_total-m))~=interleave_output(1,1:1:2*(L_total-m))));  for i=1:user   %Decoder one               L_a(i,alpha(i,:)) = L_check(i,:);%-L_e(i,:);  % a priori info. %         if dec_alg == 0            L_all(i,:) = logmapo(L_a(i,:), g, L_a(i,:), 1);  % complete info.            %         else   %            L_all(i,:) = sova0(rec_s(2*i-1,:), g, L_a(i,beta(i,:)), 1);  % complete info.               L_k(i,:) = L_all(i,:)- L_a(i,:);                        L_e(i,:) = L_k(i,alpha(i,:));            end        % extrinsic info.              %   L_all = L_all+L_a;% Decoder two         %         L_a(i,:) = L_e(i,alpha);  % a priori info.%         if dec_alg == 0%            [L_all(i,:),L_check(i,:)] = logmapo(rec_s(2*i,:), g, L_a(i,:), 2);  % complete info.  %         else%            L_all(i,:) = sova0(rec_s(2*i,:), g, L_a(i,:), 2);  % complete info. %         end%         L_e(i,beta(i,:)) = L_all(i,:) - 2*rec_s(2*i,1:2:2*L_total) - L_a(i,:);  % extrinsic info.               %end i     % Estimate the info. bits                % xhat = (sign(L_all)+1)/2;        xhat=sign(L_all);% Number of bit errors in current iteration         err(iter) = length(find(xhat(1,1:2:2*(L_total-m))~=en_output(1,1:2:2*(L_total-m))));% 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,1)==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		%whileend 		%nENdiary off