clear;close all;
fcarr=5e3;         % Carrier frequency(Hz)
N=20;		        % 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)
%_________________________
figure(1)
subplot(1,3,1)
plot(td,data2)
axis([0 1 -2 2]);
grid on
xlabel('Time')
ylabel('Amplitude')
title('輸入數據')
%_________________________________
%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
%___________________________________
subplot(1,3,2)
plot(tiq,Isymbols)
axis([0 1 -2 2]);
grid on
xlabel('Time')
ylabel('Amplitude')
title('I 支路信號')
%_________________________________

bs2=data(2:2:length(data));%even
symbols1=ones(T/Ts,1)*bs2;
Qsymbols=symbols1(:);%Q_waveform
%__________________________________
subplot(1,3,3)
plot(tiq,Qsymbols)
axis([0 1 -2 2]);
grid on
xlabel('Time')
ylabel('Amplitude')
title('Q 支路信號')

%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 
%__________________________________________________________________________

figure(2)
subplot(1,3,1)
plot(tiq,si)
axis([.498 .502 -2 2]);
grid on
xlabel('Time')
ylabel('Amplitude')
title('I 支路調制后波形')

subplot(1,3,2)
plot(tiq,sq)
axis([.498 .502 -2 2]);
grid on
xlabel('Time')
ylabel('Amplitude')
title('Q 支路調制后波形')

subplot(1,3,3)
plot(tiq,sumiq)
axis([.498 .502 -2 2]);
grid on
xlabel('Time(s)')
ylabel('Amplitude')
title('QPSK 波形')

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

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

figure(3)
subplot(1,3,1)
plot(tiq,sumiq1)
axis([.498 .502 -2 2]);
grid on
xlabel(' Time(s)')
ylabel('Amplitude')
title('疊加噪聲后的QPSK 波形')




%=============================================================
%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)

subplot(1,3,2);
plot(tiq,ycfo1);
title('濾波后 I路信號');
grid on;

subplot(1,3,3);
plot(tiq,ycfo);
title('濾波后 Q路信號');
grid on;

%=========================================================
% I CORRELATION RECEIVER COMPARATOR[ADC](after low pass filter)
%=========================================================
pt1=1.7e-8;%sets level where threshhold device comparator triggers
H=5;%(volts)
L=-2;%(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=-2;%(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=[bit1];
bitout1=[bit2];
% 
% 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 +/- 1for i=1:8
%___________________________________________
figure(4)
subplot(1,3,1)
plot(tiq,bitout);
axis([0 1 -2 2]);
grid on;
title('經判決后的I路輸出')
xlabel('   Time')
ylabel('Voltage')

subplot(1,3,2)
plot(tiq,bitout1);
axis([0 1 -2 2]);
grid on;
title('經濾波判決后的Q路輸出')
xlabel('   Time')
ylabel('Voltage')

out=[];

for i=1:1000:10000
 out=[out bitout(i:i+999) bitout1(i:i+999)];
end
out1=[];
for i=1:1000:20000
    if sum(out(i:i+999))>0
        outa(i:i+999)=1;
        out1=[out1 1];
    else
        outa(i:i+999)=-1;
        out1=[out1 -1];
    end
end    
% out=[out bitout1(i) bitout(i)];
subplot(1,3,3);
plot(td,outa)
title('解調后的QPSK信號');
grid on;
axis([0 1 -2 2]);