jpeg_turbo_sys_demo.m

结合Turbo和OFDM

% JPEGdemo.m
% Prototype JPEG compression algorithm demostration
%
% copyright (c) 1997-2002 by Yu Hen Hu
% 
% This algorithm only demonstrate the basic 
% JPEG functionalities.
% It is not necessarily a faithful 
% implementation of JPEG.
% Its output will not be binary bit streams 
% either, but rather
% an integer stream of 0 and 1s
% Only gray scale picture is considered
% 
% Last modification: 11/6/2002
clear all
clc
% Load data
 chos=0;  % default choice

    %load p64int.txt;f=p64int; clear p64int;  
    load lena.mat
    f=x;%(1+128:128+128,1+128:128+128);
    imshow(mat2gray(f))
    clear x

%echo on

% level shift by 128
 f=f-128; 
%pause
drawnow
[mf,nf]=size(f); mb=mf/8; nb=nf/8;  
% size of f, # of blocks of f

% Step 1. 2D separable DCT on each 8x8 
% blocks 
    Ff=blkproc(f,[8 8],'dct');  
    % apply DCT to each column of each block of f
    Ff=blkproc(Ff',[8 8],'dct');
    % apply DCT to each row of each block of Ff
    Ff=round(Ff');
    % transpose back to proper orientation

%pause

% Perceptual scaler quantization
%
Q =[16 11 10 16  24  40  51  61
    12 12 14 19  26  58  60  55
    14 13 16 24  40  57  69  56
    14 17 22 29  51  87  80  62
    18 22 37 56  68 109 103  77
    24 35 55 64  81 104 113  92
    49 64 78 87 103 121 120 101
    72 92 95 98 112 100 103 99];
% this is the quantization matrix shown in figure 8.37 in the textbook
%pause
% Now perform rounding

    Fq=round(blkproc(Ff,[8 8],'divq',Q));

%pause
%echo off
% DPCM of DC component, scaned row-wise 
if mb*nb > 1,
   fdc=reshape(Fq(1:8:mf,1:8:nf)',mb*nb,1);   
   fdpcm=dpcm(fdc,1);
else
   fdpcm=Fq(1,1);
end
dccof=[];
for i=1:mb*nb,
   dccof=[dccof jdcenc(fdpcm(i))];
end
%pause
%echo on

% Zig-Zag scanning of AC coefficients
z=[1   2   6   7  15  16  28  29
   3   5   8  14  17  27  30  43
   4   9  13  18  26  31  42  44
  10  12  19  25  32  41  45  54
  11  20  24  33  40  46  53  55
  21  23  34  39  47  52  56  61
  22  35  38  48  51  57  60  62
  36  37  49  50  58  59  63  64];
%pause
%echo off
acseq=[];
for i=1:mb
  for j=1:nb
    tmp(z)=Fq(8*(i-1)+1:8*i,8*(j-1)+1:8*j); 
    % tmp is 1 by 64
    eobi=max(find(tmp~=0)); %end of block index
                    % eob is labelled with 999    
    acseq=[acseq tmp(2:eobi) 999];
  end
end
accof=jacenc(acseq);
fillzero=zeros(1,173);
outbit=[dccof accof fillzero];
jepgout=reshape(outbit,590,100);
jepgout=jepgout';


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 = 593;	 % 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 0 1; 1 1 1 1];
end
%; g = [ 1 1 1; 1 0 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 = 0;
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


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(' 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);

   hang = 0;    % clear counter of transmitted frames
  while hang<100
      hang = hang + 1;    
      x = jepgout(hang,:);    % 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);

    
% 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;

% Save intermediate results 
       %  save turbo_sys_demo EbN0db ber fer
    fprintf('第 %d 行\n', hang);
   turboout(hang,:)=xhat(1:L_total-m);
   end		%while
end 		%nEN

turboout=turboout';

turbobit=reshape(turboout,1,59000);
error=xor(turbobit,outbit);
sum(error,2)


% Inverse JPEG i.e reconstruction of image lena
%clear,clc
% accof and dccof are from jpegdemo.m , run it first
dcarr=jdcdec(turbobit(1:length(dccof)));
afrom=length(dccof)+1;
acarr=jacdec(turbobit(afrom:58827));
% Assumed that image size is 256 X 256, recostruction begins
load lena.mat % To find MSE, we need to have original image
subplot 121
imshow(mat2gray(x)),title(' 原始图象 ')
drawnow
Q =[16 11 10 16  24  40  51  61
    12 12 14 19  26  58  60  55
    14 13 16 24  40  57  69  56
    14 17 22 29  51  87  80  62
    18 22 37 56  68 109 103  77
    24 35 55 64  81 104 113  92
    49 64 78 87 103 121 120 101
    72 92 95 98 112 100 103 99];

z=[1   2   6   7  15  16  28  29
   3   5   8  14  17  27  30  43
   4   9  13  18  26  31  42  44
  10  12  19  25  32  41  45  54
  11  20  24  33  40  46  53  55
  21  23  34  39  47  52  56  61
  22  35  38  48  51  57  60  62
  36  37  49  50  58  59  63  64];
z=z(:);
mb=256/8; nb=256/8;  % Number of blocks

Eob=find(acarr==999);
kk=1;ind1=1;n=1;
for ii=1:mb
    for jj=1:nb
        ac=acarr(ind1:Eob(n)-1);
        ind1=Eob(n)+1;
        n=n+1;
        ri(8*(ii-1)+1:8*ii,8*(jj-1)+1:8*jj)=dezz([dcarr(kk) ac zeros(1,63-length(ac))]);
        kk=kk+1;
    end
end

iFq=round(blkproc(ri,[8 8],'idivq',Q));
iFf=blkproc(iFq,[8 8],'idct2');  
iFf=round(iFf+128);
subplot 122
imshow(mat2gray(iFf)),title(' 还原图象 ')
% Calculate MSE , SNR
MSE=mean(mean((x-iFf).^2))                 % Doubt about formulae 
SNR=10*log10(255^2/MSE)       % Doubt about formulae

MSE =

   34.8418


SNR =

   32.7098