% 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 errors
% actually occurred (this should be commented out 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