turbo_sys_demo.m
上传用户:hnyfjx
上传日期:2013-06-30
资源大小:2149k
文件大小:6k
源码类别:

传真(Fax)编程

开发平台:

Matlab

  1. % This script simulates the classical turbo encoding-decoding system. 
  2. % It simulates parallel concatenated convolutional codes.
  3. % Two component rate 1/2 RSC (Recursive Systematic Convolutional) component encoders are assumed.
  4. % First encoder is terminated with tails bits. (Info + tail) bits are scrambled and passed to 
  5. % the second encoder, while second encoder is left open without tail bits of itself.
  6. %
  7. % Random information bits are modulated into +1/-1, and transmitted through a AWGN channel.
  8. % Interleavers are randomly generated for each frame.
  9. %
  10. % Log-MAP algorithm without quantization or approximation is used.
  11. % By making use of ln(e^x+e^y) = max(x,y) + ln(1+e^(-abs(x-y))),
  12. % the Log-MAP can be simplified with a look-up table for the correction function.
  13. % If use approximation ln(e^x+e^y) = max(x,y), it becomes MAX-Log-MAP.
  14. %
  15. % Copyright Nov 1998, Yufei Wu
  16. % MPRG lab, Virginia Tech.
  17. % for academic use only
  18. clear all
  19. % Write display messages to a text file
  20. diary turbo_logmap.txt
  21. % Choose decoding algorithm 
  22. dec_alg = input(' Please enter the decoding algorithm. (0:Log-MAP, 1:SOVA)  default 0    ');
  23. if isempty(dec_alg)
  24.    dec_alg = 0;
  25. end
  26. % Frame size
  27. L_total = input(' Please enter the frame size (= info + tail, default: 400)   ');
  28. if isempty(L_total)
  29.    L_total = 400;  % infomation bits plus tail bits
  30. end
  31. % Code generator
  32. g = input(' Please enter code generator: ( default: g = [1 1 1; 1 0 1 ] )      ');
  33. if isempty(g)
  34.    g = [ 1 1 1;
  35.          1 0 1 ];
  36. end
  37. %g = [1 1 0 1; 1 1 1 1];
  38. %g = [1 1 1 1 1; 1 0 0 0 1];
  39. [n,K] = size(g); 
  40. m = K - 1;
  41. nstates = 2^m;
  42. %puncture = 0, puncturing into rate 1/2; 
  43. %puncture = 1, no puncturing
  44. puncture = input(' Please choose punctured / unpunctured (0/1): default 0     ');
  45. if isempty(puncture) 
  46.     puncture = 0;
  47. end
  48. % Code rate
  49. rate = 1/(2+puncture);   
  50. % Fading amplitude; a=1 in AWGN channel
  51. a = 1; 
  52. % Number of iterations
  53. niter = input(' Please enter number of iterations for each frame: default 5       ');
  54. if isempty(niter) 
  55.    niter = 5;
  56. end   
  57. % Number of frame errors to count as a stop criterior
  58. ferrlim = input(' Please enter number of frame errors to terminate: default 15        ');
  59. if isempty(ferrlim)
  60.    ferrlim = 15;
  61. end   
  62. EbN0db = input(' Please enter Eb/N0 in dB : default [2.0]    ');
  63. if isempty(EbN0db)
  64.    EbN0db = [2.0];
  65. end
  66. fprintf('nn----------------------------------------------------n'); 
  67. if dec_alg == 0
  68.    fprintf(' === Log-MAP decoder === n');
  69. else
  70.    fprintf(' === SOVA decoder === n');
  71. end
  72. fprintf(' Frame size = %6dn',L_total);
  73. fprintf(' code generator: n');
  74. for i = 1:n
  75.     for j = 1:K
  76.         fprintf( '%6d', g(i,j));
  77.     end
  78.     fprintf('n');
  79. end        
  80. if puncture==0
  81.    fprintf(' Punctured, code rate = 1/2 n');
  82. else
  83.    fprintf(' Unpunctured, code rate = 1/3 n');
  84. end
  85. fprintf(' iteration number =  %6dn', niter);
  86. fprintf(' terminate frame errors = %6dn', ferrlim);
  87. fprintf(' Eb / N0 (dB) = ');
  88. for i = 1:length(EbN0db)
  89.     fprintf('%10.2f',EbN0db(i));
  90. end
  91. fprintf('n----------------------------------------------------nn');
  92.     
  93. fprintf('+ + + + Please be patient. Wait a while to get the result. + + + +n');
  94. for nEN = 1:length(EbN0db)
  95.    en = 10^(EbN0db(nEN)/10);      % convert Eb/N0 from unit db to normal numbers
  96.    L_c = 4*a*en*rate;  % reliability value of the channel
  97.    sigma = 1/sqrt(2*rate*en);  % standard deviation of AWGN noise
  98. % Clear bit error counter and frame error counter
  99.    errs(nEN,1:niter) = zeros(1,niter);
  100.    nferr(nEN,1:niter) = zeros(1,niter);
  101.    nframe = 0;    % clear counter of transmitted frames
  102.    while nferr(nEN, niter)<ferrlim
  103.       nframe = nframe + 1;    
  104.       x = round(rand(1, L_total-m));    % info. bits
  105.       [temp, alpha] = sort(rand(1,L_total));        % random interleaver mapping
  106.       en_output = encoderm( x, g, alpha, puncture ) ; % encoder output (+1/-1)
  107.           
  108.       r = en_output+sigma*randn(1,L_total*(2+puncture)); % received bits
  109.       yk = demultiplex(r,alpha,puncture); % demultiplex to get input for decoder 1 and 2
  110.       
  111. % Scale the received bits      
  112.       rec_s = 0.5*L_c*yk;
  113. % Initialize extrinsic information      
  114.       L_e(1:L_total) = zeros(1,L_total);
  115.       
  116.       for iter = 1:niter
  117. % Decoder one
  118.          L_a(alpha) = L_e;  % a priori info. 
  119.          if dec_alg == 0
  120.             L_all = logmapo(rec_s(1,:), g, L_a, 1);  % complete info.
  121.          else   
  122.             L_all = sova0(rec_s(1,:), g, L_a, 1);  % complete info.
  123.          end   
  124.          L_e = L_all - 2*rec_s(1,1:2:2*L_total) - L_a;  % extrinsic info.
  125. % Decoder two         
  126.          L_a = L_e(alpha);  % a priori info.
  127.          if dec_alg == 0
  128.             L_all = logmapo(rec_s(2,:), g, L_a, 2);  % complete info.  
  129.          else
  130.             L_all = sova0(rec_s(2,:), g, L_a, 2);  % complete info. 
  131.          end
  132.          L_e = L_all - 2*rec_s(2,1:2:2*L_total) - L_a;  % extrinsic info.
  133.          
  134. % Estimate the info. bits        
  135.          xhat(alpha) = (sign(L_all)+1)/2;
  136. % Number of bit errors in current iteration
  137.          err(iter) = length(find(xhat(1:L_total-m)~=x));
  138. % Count frame errors for the current iteration
  139.          if err(iter)>0
  140.             nferr(nEN,iter) = nferr(nEN,iter)+1;
  141.          end   
  142.       end %iter
  143.       
  144. % Total number of bit errors for all iterations
  145.       errs(nEN,1:niter) = errs(nEN,1:niter) + err(1:niter);
  146.       if rem(nframe,3)==0 | nferr(nEN, niter)==ferrlim
  147. % Bit error rate
  148.          ber(nEN,1:niter) = errs(nEN,1:niter)/nframe/(L_total-m);
  149. % Frame error rate
  150.          fer(nEN,1:niter) = nferr(nEN,1:niter)/nframe;
  151. % Display intermediate results in process  
  152.          fprintf('************** Eb/N0 = %5.2f db **************n', EbN0db(nEN));
  153.          fprintf('Frame size = %d, rate 1/%d. n', L_total, 2+puncture);
  154.          fprintf('%d frames transmitted, %d frames in error.n', nframe, nferr(nEN, niter));
  155.          fprintf('Bit Error Rate (from iteration 1 to iteration %d):n', niter);
  156.          for i=1:niter
  157.             fprintf('%8.4e    ', ber(nEN,i));
  158.          end
  159.          fprintf('n');
  160.          fprintf('Frame Error Rate (from iteration 1 to iteration %d):n', niter);
  161.          for i=1:niter
  162.             fprintf('%8.4e    ', fer(nEN,i));
  163.          end
  164.          fprintf('n');
  165.          fprintf('***********************************************nn');
  166. % Save intermediate results 
  167.          save turbo_sys_demo EbN0db ber fer
  168.       end   %iter
  169.       
  170.    end %while
  171. end  %nEN
  172. diary off