%QPSK Transmitter and I-Q Correlation Receiver
%JC 12/23/05-revised 2/14/06
%Run from editor debug(F5)
%m-file for simulating a QPSK transmitter and receiver by modulating with a pseudo
%random bit stream. A serial to parallel conversion of the pseudo random
%bit stream is performed with mapping of two bits per symbol(phase). A cosine and
%sine carrier is configured and the I and Q symbols modulate these
%carriers via mixers. The I and Q carriers are combined and time and frequency domain
%plots are provided showing key waveforms at various positions in the QPSK
%transmitter and correlation receiver. A parallel to serial conversion is used on the
%output of the receiver. The simulation uses a serial "passband" approach.
%In  other words a carrier is used. Notes are provided at the end of the m-file.
%===================================================================
clear;
fcarr=5e3;         % Carrier frequency(Hz)
N=8;		        % Number of data bits(bit rate)
fs=20*1e3;		    % Sampling frequency
Fn=fs/2;            % Nyquist frequency
Ts=1/fs;	        % Sampling time = 1/fs
T=1/N;		        % Bit time
randn('state',0);   % Keeps PRBS from changing on reruns
td=[0:Ts:(N*T)-Ts]';% Time vector(data)(transpose)
%===================================================================
% The Transmitter.
%===================================================================
data=sign(randn(N,1))';%transpose
data1=ones(T/Ts,1)*data;
data2=data1(:);

%display input data bits in command window
data_2=data2';%transpose
data_2=data_2 >0;
Transmitted_data_bits=data_2(1:(fs)/N:end)

%Serial to parallel (alternating)
tiq = [0:Ts*2:(N*T)-Ts]';% Time vector for I and Q symbols(transpose)

bs1=data(1:2:length(data));%odd
symbols=ones(T/Ts,1)*bs1;
Isymbols=symbols(:);%I_waveform

bs2=data(2:2:length(data));%even
symbols1=ones(T/Ts,1)*bs2;
Qsymbols=symbols1(:);%Q_waveform

%generate carrier waves
%cosine and sine wave
%2 pi fc t is written as below
twopi_fc_t=(1:fs/2)*2*pi*fcarr/fs; 
a=1;
%phi=45*(pi/180)
phi=0;%phase error
cs_t = a * cos(twopi_fc_t + phi);
sn_t = a * sin(twopi_fc_t + phi);

cs_t=cs_t';%transpose
sn_t=sn_t';%transpose
si=cs_t.*Isymbols;%multiply I bitstream with cosine
sq=sn_t.*Qsymbols;%multiply Q bitstream with sine
sumiq=si+sq;%transmitter output
sumiq=.7*sumiq;%reduce gain to keep output at +/- one 

%=============================================================
%Noise
%=============================================================

noise=randn(size(sumiq));
SNR=10%set SNR in dB
constant=std(sumiq)/(std(noise)*10^(SNR/20));
sumiq=sumiq + noise*constant;
noise1=noise*constant;

%=============================================================
%Receiver(balanced modulators and low pass filters) 
%=============================================================
sig_rx1=sumiq.*cs_t;%cosine 
sig_rx1=.707.*sig_rx1;%keep output at 1Vp-p

%simple low pass filter
rc1=.01989316;%time constant
ht1=(1/rc1).*exp(-tiq/rc1);%impulse response
ycfo1=filter(sig_rx1,1,ht1)/fs;
Bit_rate=N
IFilterfreg_3dB=1/(2*pi*rc1)


sig_rx=sumiq.*sn_t;%sine
sig_rx=.707.*sig_rx;%keep output at 1Vp-p

%simple low pass filter
rc=.01989316;%time constant-
ht=(1/rc).*exp(-tiq/rc);%impulse response
ycfo=filter(sig_rx,1,ht)/fs;
Bit_rate=N
QFilterfreg_3dB=1/(2*pi*rc)
%=========================================================
% I CORRELATION RECEIVER COMPARATOR[ADC](after low pass filter)
%=========================================================
pt1=1.7e-8;%sets level where threshhold device comparator triggers
H=5;%(volts)
L=0;%(volts)
LEN=length(ycfo);
for ii=1:LEN;
    if ycfo(ii)>=pt1;%correlated output(ycfo) going above pt1 threshold setting
        pv1i(ii)=H;%I pulse voltage
    else;
        pv1i(ii)=L;
    end;
end ;
po1i=pv1i;%pulse out=pulse voltage


%=========================================================
% Q CORRELATION RECEIVER COMPARATOR[ADC](after low pass filter)
%=========================================================
pt2=1.7e-8;%sets level where threshhold device comparator triggers
H=5;%(volts)
L=0;%(volts)
LEN=length(ycfo1);
for ii=1:LEN;
    if ycfo1(ii)>=pt2;%correlated output(ycfo1) going above pt2 threshold setting
        pv2q(ii)=H;% Q pulse voltage
    else;
        pv2q(ii)=L;
    end;
end ;
po1q=pv2q;%pulse out=pulse voltage


bit1=sign(po1q);%0 and 1
bit2=sign(po1i);%0 and 1
bit3=bit1 >0;%0 and 1
bit4=bit2 >0;%0 and 1
bitout=[bit3];
bitout1=[bit4];
%==========================================================================
%Parallel to serial bitstream(uses interleaving and concatenation{joining})
%==========================================================================
    
bitout2=[bitout];
x=1128;%x=fs/N;%This is a cluge way to program but x is required to make the parallel
%to serial converter work if one changes the basic parameters such as N,fs,etc.
%x=N*(bit3 # 1's or 0's in first bit time)-fs:x=(8*2641)-20000=1128
bitout2=bitout2(1:(fs+x)/N:end);
bitout2=[bitout2];
bitout3=[bitout1];
bitout3=bitout3(1:(fs+x)/N:end);
bitout3=[bitout3];
bitfinalout=[bitout2;bitout3];
bitfinalout=bitfinalout(1:end);

%display received output data bits in command window
Received_data_bits=bitfinalout

%Received data output
data1a=ones(T/Ts,1)*bitfinalout;
bitfinal1=data1a(:);
bitfinal1=bitfinal1-mean(bitfinal1);
bitfinal1=2*bitfinal1;%get to +/- 1

%=====================================================================
%Plots
%======================================================================
figure(1)
subplot(3,2,1)
plot(td,data2)
axis([0 1 -2 2]);
grid on
xlabel('                                                          Time')
ylabel('Amplitude')
title('Input Data')

subplot(3,2,3)
plot(tiq,Isymbols)
axis([0 1 -2 2]);
grid on
xlabel('                                                          Time')
ylabel('Amplitude')
title('I Channel(one bit/symbol(phase)) Data')

subplot(3,2,5)
plot(tiq,Qsymbols)
axis([0 1 -2 2]);
grid on
xlabel('                                                          Time')
ylabel('Amplitude')
title('Q Channel(one bit/symbol(phase)) Data')

subplot(3,2,2)
plot(tiq,si)
axis([.498 .502 -2 2]);
grid on
xlabel('                                                          Time')
ylabel('Amplitude')
title('I Channel Modulated Waveform')

subplot(3,2,4)
plot(tiq,sq)
axis([.498 .502 -2 2]);
grid on
xlabel('                                                          Time')
ylabel('Amplitude')
title('Q Channel Modulated Waveform')

subplot(3,2,6)
plot(tiq,sumiq)
axis([.498 .502 -2 2]);
grid on
xlabel('                                                          Time')
ylabel('Amplitude')
title('QPSK Output Waveform')


%========================================================================
%Take FFT of modulated carrier
%========================================================================
y=sumiq;
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;
%=========================================================================
%Plot frequency domain
%=========================================================================
figure(2)
subplot(3,1,1); plot(f1,MY);xlabel('');ylabel('AMPLITUDE');
axis([4500 5500 -.5 1]);%zoom in/out
title('Frequency domain plots');
grid on

subplot(3,1,2); plot(f1,20*log10(abs(MY).^2));xlabel('FREQUENCY(Hz)');ylabel('DB');
axis([4000 6000 -80 10]);%zoom in/out
grid on
title('Modulated QPSK carrier')


figure(3)
subplot(3,2,1);
plot(td,bitfinal1)
title('Received output data');
grid on;
axis([0 1 -2 2]);

subplot(3,2,3);
plot(tiq,ycfo1);
title('Filtered I Channel Data');
grid on;

subplot(3,2,5);
plot(tiq,ycfo);
title('Filtered Q Channel Data');
grid on;

subplot(3,2,2);
plot(tiq,sig_rx1);
grid on;
axis([0 1 -1 1]);%zoom in/out
title('Unfiltered I Channel Output');


subplot(3,2,4);
plot(tiq,sig_rx);
grid on;
axis([0 1 -1 1]);%zoom in/out
title('Unfiltered Q Channel Output');

phasevl=atan2(ycfo1,ycfo);%phases of transmitter output
                          %phase deviation---(-2.355 radians)=-135 degrees(-1-1)
                          %phase deviation---(+0.785 radians)=+45 degrees(1 1)                          %1)
                          %phase deviation---(-0.785 radians)=-45 degrees(-1 1)
                          %phase deviation---(+2.355 radians)=+135 degrees(1-1)
subplot(3,2,6);
plot(tiq,phasevl);
grid on;
title('Output phase deviation(radians) levels')
xlabel('   Time')
ylabel('Radians')

figure(4);
subplot(4,2,1);
plot(tiq,noise1);%plot noise times a constant
grid on;
title('Noise output')
xlabel('   Time')
ylabel('Voltage')

subplot(4,2,3)
plot(tiq,po1i);
axis([0 1 -1 6]);
grid on;
title('I channel ADC output')
xlabel('   Time')
ylabel('Voltage')

subplot(4,2,5)
plot(tiq,po1q);
axis([0 1 -1 6]);
grid on;
title('Q channel ADC output')
xlabel('   Time')
ylabel('Voltage')

figure(5);%These values were obtained from program discussed in given reference with
%500000 symbols on each run.
err=[.2921 .2450 .1970 .1511 .1099 .0739 .0457 .0248 .0121 .0047 .0017 3.6e-4 7.8e-5];
snr=10.^ ( [0:0.1:12]./10);
Pbsym=erfc(.707*sqrt(snr)); % symbolBER (Theoretical-QPSK)
Pbbit=0.5*erfc(sqrt(snr)); % bitBER (Theoretical-BPSK)
semilogy([0:12],err,'*',[0:0.1:12],Pbsym,'-',[0:0.1:12],Pbbit,'--'); % plot
grid on; xlabel('SNR(dB)'); ylabel('BER');
title('Simulation of symbol and bit BER for QPSK and BPSK');
legend('QPSK Simulated','QPSKsymbol Theoretical','BPSKbit Theoretical');



%NOTES
%Serial to parallel conversion of a serial bit stream and mapping of
%two bits to a symbol(phase) can sometimes be confusing. I will try and explain
%with an example.
%Suppose you have a serial bit stream of ten  0 0 1 1 0 1 1 0 1 1 even # of bits      
                                %odd bits     0   1   0   1   1
                                %even bits      0   1   1   0   1
                                
%The odd bits are the I Channel Data at one half the original serial bit stream
%bit rate. Notice that the amplitudes are +/- one as shown in figure 1.
%The possible combinations(changing zero to -1) are -1 -1, 1 1, -1 1, 1 -1(four phases or four symbols).
%The amplitudes, in theory, should be held at +/- 0.707 to keep the summed output of the 
%QPSK transmitter at a  constant amplitude of +/- one.
%The even bits are the Q Channel Data at one half the original serial bit
%stream bit rate. The odd bits are the I Channel Data at one half the
%original serial bit streamm bit rate.

%Things to do:
%Implement BER code to prove that the BER of the output of either the I or
%Q channel is equal to the BER of BPSK. Also prove that the symbol
%BER(combined output) is ~ 3DB poorer than the I or Q channel output.
%Implement Gray coding and prove that the symbol BER approaches the BER of
%the I or Q channels at high S/N ratios.
%Implement different types of low pass filters for best BER. Appropriate TX and RX
%bandpass filters could also be added for best BER.

%==========================================================================
%Conclusions:
%==========================================================================
%I've come to the conclusion that accurate BER's can't be performed for
%QPSK at passband using this m-file.If someone figures out how to do it, let me know.
%Simulations must be setup using a baseband approach m-file construction. In other words,
%complex symbol inputs must be used. I've observed that even the Communication's Toolbox does
%not do a passband simulation for QPSK. This m-file, as configured, gives a good overview
%of how a real world QPSK systom would be configured. I've added several slicers
%(I call them one bit ADC'S) for drill and provided theoritical BER curves
%for bitBPSK and symbolQPSK systems. The following reference provides a baseband m-file 
%that uses a complex symbol input baseband approach and gives an accurate(theoretical
%matches simulated) symbol BER for QPSK. I have run the program with 500000
%symbols and put the results on figure 5 with *,s.
%Ref:www.cs.tut.fi/kurssit/ 83050/matlab_tehtavat/83050_M2.pdf 
%The math for probability of error for QPSK using distance(signal space) criteria from
%constellations shows the maximum distance between symbols for BPSK and
%QPSK to be the same and gives the same probability of error for each
%which means you can get something for nothing which is not true. QPSK is
%~3 dB poorer than BPSK because BPSK has only one possible bit for error where QPSK has 2
%bit per symbol that can be in error. If you use Gray coding QPSK can
%approach the BER of BPSK at higher SNR's. In a nutshell, always use a QPSK
%design with Gray coding(assuming you can handle the added design
%complexity) as you can send twice as many bits in the same amount of RF
%spectrum. QPSK is much more efficient than BPSK.