% Andrew Puryear
% 362 Project - vcode
% This matlab program simulates the performance of an convolutionaly encoded 
%   baseband antipodal image signaling over an additive white Gaussian noise channel.

clear;
load F:\MATLAB6p1\work\PSYCH221\Project\TestPictures\einstein.mat;
im = round(X/max(X(:))*256)-1;
Q_factor = 50;
qtable = 1./makeQTable(Q_factor);

% JPEG encode here
properSize = ceil(size(im)/8) * 8;
imfull = zeros(properSize);
imfull(1:size(im,1), 1:size(im,2)) = im;

% If the image values are between 0 and 1, scale to between 0 and
% 255
if (max(imfull(:))<=1)
  imfull = round(imfull*255);
end

% make the dct basis functions
n = 8;
c = [ 1/sqrt(2) 1 1 1 1 1 1 1 ];
j = 0:n-1;
dctMatrix = zeros(n,n);
for u = 0:n-1
  dctMatrix(u+1,:) = (2*c(u+1) / n)* cos( (2*j+1) * u * pi / (2*n));
end
idctMatrix = 4*dctMatrix';

% DCT transformation and quntization
dctQuant = zeros(size(imfull));

for i=1:8:size(imfull,1)
  for j=1:8:size(imfull,2)
    block = imfull(i:i+7, j:j+7);
    dctCoef = dctMatrix*block*dctMatrix';
    dctQuant(i:i+7, j:j+7) = round(dctCoef ./ qtable) .* qtable;
  end
end
im= dctQuant-min(dctQuant(:));
im= round(im/max(im(:))*255);
a=[];
for k=1:size(im,1)
    k
    for p=1:size(im,2)
        a=[a,dec2bin(im(k,p),8)];
    end
end
for k=1:length(a)
    sig(k)=bin2dec(a(k)); 
end
disp('Hard Decision Decoding');
for db=[2:2:10]   
    windL=15;
    out=[];
    ts = 0; 			% The timing error in 10^-1 seconds 
    
    % Create binary image stream.
    inter=[de2bi(conv(sig,[1,0,1]));de2bi(conv(sig,[1,1,1]))];  
    
    % convolutionaly encode and interleave data sequence to create sequence for transmission
    encode(1:2:2*no_sigs)=inter(1:no_sigs,1);
    encode(2:2:2*no_sigs)=inter((no_sigs+3):(2*no_sigs+2),1);
    
    % modulate encoded data sequence with root raised cosine pulse shaping
    alpha = .5;
    signal = zeros (1,no_sigs*10); 
    signal (1:10:no_sigs*2*10) = round(encode-.5);
    t = [-5: 0.1 :5]; 
    u = (4.*alpha.*cos((1+alpha) .*pi.*t)+pi.*(1-alpha) .*sinc((1-alpha) .*t)) ./(pi.*(1-(4. *alpha.*t) .^2)+eps); 
    if 1/(4*alpha)*10 == round(1/(4*alpha) *10), % Take care of 'undefined' values 
        u(1/(4*alpha)*10+51)=-(2*cos(pi/(4*alpha)+pi/4)-sin(pi/(4*alpha)+pi/4)*pi)*alpha/pi;
        u(-1/(4*alpha)*10+51)=-(2*cos(pi/(4*alpha)+pi/4)-sin(pi/(4*alpha)+pi/4)*pi)*alpha/pi;
    end; 
    u=u/norm(u); % Normalize to unit energy
    pulse=filter(u,1,signal); % Create analog transmission signal 
    
    % model transmission over channel
    var = 1/sqrt(2*10^(db/10)); 
    AWGN = randn(1,length(pulse))*var;
    get = pulse + AWGN;                 % Add noise to signal 
    
    % model receiver
    rec=filter(u,1,get); % Matched filter receiver 
    no_sigs=no_sigs-ts;
    htrans=sign(rec(101+ts:10:no_sigs*2*10-9+ts))/2+0.5; %create phase estimation error and hard decision decode.
    N=length(htrans)/2;
    
    % Viterbi hard decision decoding
    mem=zeros(4,windL);
    cur=zeros(4,1);
    % create and initalize 'current state'
    cur(4)=xor(htrans(1),1)+xor(htrans(2),1)+xor(htrans(3),1)+xor(htrans(4),0);
    cur(3)=xor(htrans(1),0)+xor(htrans(2),0)+xor(htrans(3),1)+xor(htrans(4),1);
    cur(2)=xor(htrans(1),1)+xor(htrans(2),1)+xor(htrans(3),0)+xor(htrans(4),1);
    cur(1)=xor(htrans(1),0)+xor(htrans(2),0)+xor(htrans(3),0)+xor(htrans(4),0);
    count=1;
    for cnt=5:2:2*N
        % Shift columns of mem
        mem(:,2:windL)=mem(:,1:windL-1);
        %figure out which transition is the best
        if (cur(3)+xor(htrans(cnt),1)+xor(htrans(cnt+1),0))<(cur(4)+xor(htrans(cnt),0)+xor(htrans(cnt+1),1))
            mem(4,1)=3;
        else
            mem(4,1)=4;
        end
        if (cur(1)+xor(htrans(cnt),1)+xor(htrans(cnt+1),1))<(cur(2)+xor(htrans(cnt),0)+xor(htrans(cnt+1),0))
            mem(3,1)=1;
        else
            mem(3,1)=2;
        end
        if (cur(3)+xor(htrans(cnt),0)+xor(htrans(cnt+1),1))<(cur(4)+xor(htrans(cnt),1)+xor(htrans(cnt+1),0))
            mem(2,1)=3;
        else
            mem(2,1)=4;
            end
        if (cur(1)+xor(htrans(cnt),0)+xor(htrans(cnt+1),0))<(cur(2)+xor(htrans(cnt),1)+xor(htrans(cnt+1),1))
            mem(1,1)=1;
        else
            mem(1,1)=2;
        end
        % update hamming distance for each state
        curn(4)=min([cur(3)+xor(htrans(cnt),1)+xor(htrans(cnt+1),0),cur(4)+xor(htrans(cnt),0)+xor(htrans(cnt+1),1)]);
        curn(3)=min([cur(1)+xor(htrans(cnt),1)+xor(htrans(cnt+1),1),cur(2)+xor(htrans(cnt),0)+xor(htrans(cnt+1),0)]);
        curn(2)=min([cur(3)+xor(htrans(cnt),0)+xor(htrans(cnt+1),1),cur(4)+xor(htrans(cnt),1)+xor(htrans(cnt+1),0)]);
        curn(1)=min([cur(1)+xor(htrans(cnt),0)+xor(htrans(cnt+1),0),cur(2)+xor(htrans(cnt),1)+xor(htrans(cnt+1),1)]);
        cur=curn;
        flipud(mem);
        % pop out data point from end of window
        if cnt>(windL*2+1)
            if cur(4)==max(cur)
                sstate=4;
            end
            if cur(3)==max(cur)
                sstate=3;
            end
            if cur(2)==max(cur)
                sstate=2;
            end
            if cur(1)==max(cur)
                sstate=1;
            end
            for cnter=[1:windL]
                sstate1=sstate;
                sstate=mem(sstate,cnter);
            end
            if sstate1>sstate
                out(count)=1;
            end
            if sstate1<sstate
                out(count)=0;    
            end
            if sstate1==4&sstate==4
                out(count)=1;
            end
            if sstate1==1&sstate==1
                out(count)=0;            
            end
            count=count+1;
        end
        % display picture
        s=0;
        for k=1:size(im,1)
            for p=1:size(im,2)
                s=s+1;
                image(k,p)=out(s*8-7)*2^7+out(s*8-6)*2^6+out(s*8-5)*2^5+out(s*8-4)*2^4+out(s*8-3)*2^3+out(s*8-2)*2^2+out(s*8-1)*2^1+out(s*8);
            end
        end
        image=idct2(image);
        image=image-min(image(:));
        image=round(image/max(image(:))*255)
        imshow(image/255,255);            
    end
    runs=runs+1;
end
%******************************************************************************************************%