%-- DEMOD.M: a complex-PAM receiver ---------------------------------------%
%--------------------------------------------------------------------------%

%-- Load Shared parameters generated by SEND.M ----------------------------%
%--------------------------------------------------------------------------%

load Share/alphabet;    M = length(alphabet);
load Share/parms;    baud = parms(1);      f0 = parms(2); Ltrain = parms(3);
                        L = parms(4); timelag = parms(5);   rows = parms(6);
		     cols = parms(7);
Rsample = 8192;         q = round(Rsample/baud); % 8192 = rcvr sampling rate

%-- Receiver parameters ---------------------------------------------------%
%--------------------------------------------------------------------------%
 
BT = 0.6;		% BT = {LPF bandwidth}/{baud rate}, 1/2 < BT < q / 2
Lsync = 30;		% number of training symbols used for synchronization
 
%-- Load received signal --------------------------------------------------%
%--------------------------------------------------------------------------%
 
r = auread('Share/recorded.au'); r = r / max(abs(r));
load Cheating/s; clf; subplot(221); psd(s,512,Rsample); hold on; 
                                    psd(r,512,Rsample); hold off;
ax = axis; axis([f0-3*8192/q, f0+3*8192/q, ax(3:4)]); 
title('Transmit and Receive PSD');  ylabel('dB'); drawnow;

%-- Downconvert -----------------------------------------------------------%
%--------------------------------------------------------------------------%

disp('Downconverting ...');
r_tilde = exp(-2i*pi*f0*((1:length(r))).' / Rsample) .* hilbert(r);

%-- Find start index (i0) of training tone --------------------------------%
%--------------------------------------------------------------------------%
disp('Synchronizing ...');

% regenerate env of training tone  ----

rand('state',0); train = alphabet(ceil(M*rand(Ltrain,1)));
trainchips = [train.'; zeros(q-1,Ltrain)];
tone_tilde = conv(sinc((-timelag*q : timelag*q) / q).', trainchips(:));

% correlate with complex env of r  ----

mf = conj(flipud(tone_tilde((timelag*q) + 1:(min(Lsync,Ltrain)*q))));
corr = conv(r_tilde(1:floor(length(r_tilde)/4)), mf);
[tmp, maxindex] = max(abs(corr));
i0 = maxindex - length(mf) + 1;

%-- Demodulate ------------------------------------------------------------%
%--------------------------------------------------------------------------%

tbyTs = (-timelag*q : timelag*q);
lpf = sinc(2 * BT * tbyTs / q);

disp('Demodulating ...');
for k = Ltrain+(1:L),
        z(k) = lpf * r_tilde(i0 + (k - 1)*q - tbyTs);
        dec(k-Ltrain) = alphabet((z(k) - alphabet)==min(z(k) - alphabet));
end;
subplot(223); plot(z(Ltrain+(1:L)),'.k'); axis square; axis equal;

%-- Convert decisions back to image ---------------------------------------%
%--------------------------------------------------------------------------%

Bits = zeros(length(dec),log2(M));
for i=1:length(alphabet), for j=1:log2(M),
        ix = dec==alphabet(i);
        if any(ix), Bits(ix,j) = bitget(i-1,j); end;
end; end;
Bits = Bits(:);
rand('state',0); scrambler = rand(length(Bits),1)>0.5;
Xr = reshape(reshape(xor(Bits(1:rows*cols*8),scrambler(1:rows*cols*8)),...
              rows*cols,8)*(2.^(0:7).'),rows,cols);
clear ix scrambler;
 
load Cheating/X;
subplot(222); imshow(uint8(X));  title('original');
subplot(224); imshow(uint8(Xr)); title('received');
 
%-- Calculate Bit-Error Rate ----------------------------------------------%
%--------------------------------------------------------------------------%
 
load Cheating/bits;
BER = sum(Bits ~= bits(:)) / length(Bits), clear bits

disp(['Bit rate = log2(', num2str(M), ') * (8192 / ', num2str(q), ') = ',...
num2str(log2(M)),' * ', num2str(8192/q),' = ', ...
num2str(log2(M)*(8192/q)), ' b/s.']);