%  This script simulates a digital communication system where noise and 
%  limited bandwidth are the primary sources of signal corruption
%  resulting in bit errors.  Binary sequences are randomly generated
%  and encoded into baseband line signals. White noise is added based
%  on a specified SNR and the channel band is set according to
%  a cut-off frequency for a low-pass Butterworth filter. A matched
%  filter is derived and implemented.  The resulting bit error rate is 
%  computed along with a plot showing where in the sequence bit error
%  actually occurred (this should be disable for excessively long bit
%  streams (>10000).  The bit error rate is printed in the title of
%  the figure.  This script uses the function MODULB (which must be
%  downloaded from EE422G web site as well) to generate various line
%  codes and extract the signal templates for the match
%  filter design and implementation.
%  
%   Written by Kevin D. Donohue (donohue@engr.uky.edu) October 2007
 
clear
fs = 25000; %  samples/sec, sampling rate on output waveform
dr = 1000;  %  bits/sec, bit rate, should be less than fs/10 for a good simulation
% Signal to noise ratio in dB for AWGN
snr = -3;
%  Channel bandwidth model
bw = 1.5*dr; %  Hz
[b,a] = butter(4, bw/(fs/2));  % use a Butterworth LPF as channel Model
%  Select line code (index for cell array below)
linecode = 10;  
%  Binary sequence:
seq = round(rand(1,5000));  %  Uncomment this for a random sequence
%seq = [1 0 1 0 0 0 1 1 0 1 0 1 0];  %  Uncomment this for a fixed sequence
 
%  Line code options in cell array:
lncds = {'unipolar_nrz'; ...
         'bipolar_nrz'; ...
         '4level_nrz'; ...
         'bipolar_rz'; ...
         'unipolar_rz'; ...
         'ami'; ...
         'manchester'; ...
         'miller'; ...
         'unipolar_nyquist'; ...
         'bipolar_nyquist'};
%
% Create line signal
%
[y,t] = modulb(seq,dr,fs,lncds{linecode});
%
%  Compute scale for noise at specified SNR
%
sigrms = sqrt(mean(y.^2));  %  Signal power in rms
% scale power in noise relative to signal
npow = sigrms*10^(-snr/20);  %  Noise power in rms
%
%  Distort signal based on limited bandwidth
%
yd = filtfilt(b,a,y);
ynd = yd+npow*randn(size(yd));  %  Add scaled noise to signal 
%  
%  Develop correlation receiver and implement as a matched filter for
%  decoding the bit signals.
%
template1 = modulb(1,dr,fs,lncds{linecode}); %  Template match for 1
template0 = modulb(0,dr,fs,lncds{linecode}); %  Template match for 0
% Reverse time sequence and scale down to normalize by signal length
template1 = fliplr(template1)/length(template1);  
% Reverse time sequence and scale down to normalize by signal length
template0 = fliplr(template0)/length(template0);
 
%  Apply matched filter
b0stat = filter(template0,1,ynd);
b1stat = filter(template1,1,ynd);
%  Create a while loop to sample at bit period intervals
koff = length(template1)-1;  %  Synchronize first sample bit interval
k = koff;  % Initialize first correlation sample index
% Create a counter for indexing through the sequence bits
sampinc = 0;
%  Loop through while correlation sample point is less than waveform lengths
%  and there are more bits to decode in sequence.
while k < length(b1stat) && sampinc < length(seq)
    %  Sample the correlated templates and line signal at bit intervals
    k = koff + round(sampinc*fs/dr); %  Increment signal samples to next bit period
    sampinc = sampinc + 1;  %  Increment bit interval counter to decode next bit
    detstat(1, sampinc) = b0stat(k);  %  Sample correlation decision statistic of template 0
    detstat(2,sampinc) = b1stat(k);   %  Sample correlation decision statistic of template 1 
end
 
% Make detection decision based on which template had the highest
% correlation statistic (best match)
for n=1:length(seq)
    if detstat(1,n) >= detstat(2,n)
        dseq(n) = 0;  %  Decide 0 if template 0 was more correlated
    else
        dseq(n) = 1;  %  Decide 1 if template 1 was more correlated
    end
end
%  Compute the bit error rate
ers = xor(dseq,seq);  %  identify difference between original and detected sequence 
%  Compute the errors per bit
erat = mean(ers);
 
% 
%  Plot the positions where errors occured in the detection sequence
figure(1)
stem(ers)
xlabel('Bit Index')
title(['1 indicates position where error occurred, bit error rate = ' num2str(erat) ])
 
% %  Uncomment these next set of lines to view noisy and distorted waveform
% %  received after passing through channel
%  
% %  Plot of line signal with noise and channel distortion
% figure(2)
% plot(t,ynd)
% xlabel('Seconds')
% ylabel('Volts')
% title([lncds{linecode},' line signal representing binary sequence with channel bandwidth of ' ...
%        num2str(bw), ' snr'],'Interpreter','none')
% %  Extend y-axis to see signal ends and put text on top
% ylow = min(ynd);  % Find signal low points 
% yhigh = max(ynd);  % Find signal high points
% set(gca,'Ylim', [ylow(1)-.15, yhigh(1)+.15])  %  Add a little to both ends of y-axis
% %  Write binary sequence above signal in the center of time interval
% for k=1:length(seq)
%     xpoint = length(template0)/(2*fs)+(k-1)/dr;
%     text(xpoint, max(y)+.05,int2str(seq(k)))
% end
% %  Write lines at boundaries of each bit signal
% hold on
% for k =1:length(seq)
%     xpoint = length(template0)/(2*fs)+ (k-1.5)/dr;
%     plot([xpoint, xpoint], [ylow(1)-.15 yhigh(1)+.15], 'g:')
% end
% hold off