%Run from editor Debug(F)
%This m file(ASK.m) analyzes a coherent amplitute shift keyed(ASK) and a binary
%phase shift keyed(BPSK) communication system. The receiver uses a correlator
%(mixer-integrator [lpf]) configuration with BER measurements comparing measured
%and theoretical results. The band pass and low pass filters used in the receiver are
%constructed using z transforms. M files on BPF and LPF design using z
%transforms can be found in the author index (CJ) mathworks file exchange. A fundemental 
%question to ask is: "Why does BPSK show a 3dB improvement(as you will see using the 
%program)in BER over ASK?". A simple answer is that the signal for ASK is being
%transmitted only half the time. A reference is provided at the end of the program that
%was used in writing the program. Always remember a journey of a thousand
%miles requires a first one small step.
%========================================
% Set universal parameters
% =======================================
clear
fs=8e5;%sampling frequency
fm=20e3;%square wave modulating frequency(NRZ)= 40KHz bit rate
n=2*(6*fs/fm);
final=(1/fs)*(n-1);
fc=2e5; % carrier frequency
t=0:1/fs:(final);
Fn=fs/2;%nyquist frequency
%=========================================
%Generate square wave by using cosine wave
%==========================================
% cosine wave
% 2 pi fc t is written as below
twopi_fc_t=2*pi*fm*t; 
A=1;
phi=0;
x = A * cos(twopi_fc_t + phi);

% square wave
am=1;
x(x>0)=am;
%x(x<0)=0;%use for ASK-comment out BPSK
x(x<0)=-1;%use for BPSK-remember ASK variables become BPSK

subplot(321);
plot(t,x);
axis([1e-4 3e-4 -2 2]);
title('Square Wave Modulating Input To Transmitter');
grid on

car=sin(2*pi*fc*t);%Sinewave carrier waveform
ask=x.*car;%modulate carrier(ASK or BPSK)

subplot(322);
plot(t,ask);
axis([0 100e-6 -2 2]);
title('Modulated Sinewave Carrier Waveform');
grid on;
%=====================================================
%Noise generator SNR=Eb/No=20log(Signalrms/Noiserms)
%======================================================
%vn=0;
vn=.1;%set noise level 0.1~=6db=SNR=Eb/No
noise=vn*(randn(size(t)));%noise generator
subplot(323);
plot(t,noise);
grid on;
title('Noise Level');
axis([0 .6e-3 -1 1]);


askn=(ask+noise);%modulated carrier plus noise
subplot(324);
plot(t,askn);
axis([0 100e-6 -2 2]);
title('Modulated Carrier Waveform Plus Noise');
grid on;

%======================================================================
%Receiver bandpass filter(two poles two zeros) 
%======================================================================
fBW=40e3;
f=[0:3e3:4e5];
w=2*pi*f/fs;
z=exp(w*j);
BW=2*pi*fBW/fs;
a=.8547;%BW=2(1-a)/sqrt(a)
p=(j^2*a^2);
gain=.135;
Hz=gain*(z+1).*(z-1)./(z.^2-(p));

subplot(325);
plot(f,abs(Hz));
title('Receiver Bandpass Filter Response');
grid on;


Hz(Hz==0)=10^(8);%avoid log(0)
subplot(326);
plot(f,20*log10(abs(Hz)));
grid on;
title('Receiver -3dB Filter Response');
axis([1e5 3e5 -3 1]);

%filter coefficients
a=[1 0 .7305];%[1 0 p]
b=[.135 0 -.135];%gain*[1 0 -1]
faskn=filter(b,a,askn);
figure;
subplot(321);
plot(t,faskn);
axis([0 100e-6 -2 2]);
title('Receiver BPF Output');
grid on;


cm=faskn.*car;%multiply modulated carrier with unmodulated carrier
subplot(322);
plot(t,cm);
axis([0 100e-6 -2 2]);
grid on;
title('Receiver Multiplier Output');

%===================================================================
%Low pass filter(one pole one zero)
%==================================================================
p=.72;
gain1=.14;%gain=(1-p)/2
Hz1=gain1*(z+1)./(z-(p));
%0.65=53KHz
%0.7=45KHz
%0.72=40KHZ=best for lowest BER-Theory says set equal to the bit rate
%0.75=37KHz
%0.8=28KHz

%subplot(323);
%plot(f,abs(Hz1));
%title('Receiver LPF Response');
%grid on;

subplot(323);
Hz1(Hz1==0)=10^(-8);%avoid log(0)
plot(f,20*log10(abs(Hz1)));
grid on;
title('LPF -3dB response');
axis([0 5e4 -3 1]);

%filter coefficients
a1=[1 -.72];%(z-(p))
b1=[.14 .14];%gain*[1 1]

so=filter(b1,a1,cm);
so=so*10;%add gain
so=so-mean(so);%removes DC component
subplot(324);
plot(t,so);
axis([1e-4 3e-4 -2.5 2.5]);
title('Receiver Output Signal-LPF');
grid on;
%======================================================
%Comparator
%======================================================
High=2.5;
Low=-2.5;
vt=0;%sets comparator threshold
error=0;
%calculate the number of elements in matrix so, the output of the filter
len1=length(so); %number of elements in so
for ii=1:len1
  if so(ii) >= vt
  Vs(ii)=High;
  else
  Vs(ii)=Low;
  end
end
Vo=Vs+2.5;
subplot(325);
plot (t,Vo), title('Receiver Regenerated signal'),
axis([1e-4 3e-4 -8 8])
grid on;
xlabel('Time (s)'), ylabel('Magnitude(V)'),
%========================================================
%Generate delayed squarewave
%=========================================================
%cosine wave delayed to match receiver regenerated signal square waves.
%filters have delay
% 2 pi fc t is written as below
twopi_fc_t1=2*pi*fm*(t-.59e-4);%adds delay 
%A=1;
phi1=0;
x1 = A * cos(twopi_fc_t1 + phi1);

% square wave
x1(x1>0)=am;
x1(x1<0)=0;

sr=x1;
sf=5*sr;
subplot(326);
plot (t,sf), title('Input signal-delayed'),
axis([1e-4 3e-4 -8 8]);
grid on;
xlabel('Time (s)'), ylabel('Magnitude(V)'),

a=5*xor(sf,Vo); %show where errors occur against time.
figure(3);
plot (t,a), title('Error pulses'),
axis([0 8e-4 -8 8]);
grid on;
xlabel('Time (s)'), ylabel('Magnitude(V)'),
  %for ii=1:len1;
  for ii=20:1:len1%during startup,errors are generated. The first 20 samples of len1
                  %(480 samples) are bypassed and no noise gives zero errors. 480-20=460 
  if sf(ii) ~= Vo(ii) % if output and input are not the same then increment error
  error=error+1;
  end
end
error %display number of errors

figure(4);
err=[0 0 0 0 0 0 error/460 0 0];%Measured at different values of SNR by setting vn
snr=10.^ ( [0:0.1:12]./10);
Pb= 0.5*erfc(.707*sqrt(snr)); % BER (Theoretical-ASK)
Pb1= 0.5*erfc(sqrt(snr)); % BER (Theoretical-BPSK)
semilogy([0:8],err,'*',[0:0.1:12],Pb,'-',[0:0.1:12],Pb1,'--'); % plot
grid on; xlabel('SNR=Eb/No(dB)'); ylabel('BER');
title('Simulation of BER for ASK and BPSK');
legend('BER','ASK Theoretical','BPSK Theoretical');
BER=error/460%display BER
%Could have set up a for loop using vn to calculate all points(err) but would
%complicate matters.Easier to run by hand and observe BER plot. Presently
%set for SNR=6dB.
%===================================================
%Frequency domain plots
%====================================================
%Lets take a look at askn(unfiltered carrier + noise) and
%faskn(band pass filtered carrier + noise) output on a spectrum analyzer.
%This fft code is somewhat complex but comes from the mathworks. Works
%great as it shows one what they would actually see using a spectrum analyzer.

%y=askn;%unfiltered carrier + noise 
y=faskn;%band pass filtered carrier + noise
NFFY=2.^(ceil(log(length(y))/log(2)));
FFTY=fft(y,NFFY);%pad with zeros
NumUniquePts=ceil((NFFY+1)/2); 
FFTY=FFTY(1:NumUniquePts);
MY=abs(FFTY);
MY=MY*2;
MY(1)=MY(1)/2;
MY(length(MY))=MY(length(MY))/2;
MY=MY/length(y);
f1=(0:NumUniquePts-1)*2*Fn/NFFY;

figure(5)
subplot(2,2,1); plot(f1,MY);xlabel('FREQUENCY');ylabel('AMPLITUDE');
axis([0 3e5 -.5 1]);%zoom in/out
title('Modulated Tx carrier plus noise');
grid on;
subplot(2,2,2); plot(f1,20*log10(abs(MY).^2));xlabel('FREQUENCY');ylabel('DB');
axis([0 3e5 -60 5]);
grid on;
title('Frequency domain plots')

%Reference
%http://soe.unn.ac.uk/ocr/people/ghassemlooy/
%Prof. Z Ghassemlooy's website
%Good site for theory and m files. I structured this m file from his work
%and have added a few things. I constructed filters using z transforms
%since I don't have the filter toolbox functions.