%******************************
% 	QPSK Demodulation Test		*
%******************************
%
%	File Name: 	qpsk_test.m
%	Author:		Qin Zhang
%	Last Modified on: 6/9/1999
%
%*******************************

%************************
%	Block Phase Locking	*
%************************

clear;
load sc1005.dat;
%load a:\lecroy_1.dir\sc1007.dat;
wave=sc1005';
length(wave)
fs=1e9;

% Set the estimated frequency draft;
% equal to 1 to 100. 
% 1 corresponding to low draft, 100, high LO draft.
frequency_draft=1;		

%***********************************************
% generate local carrier for down converting
%***********************************************   

   f_lo=(480-50)*1e6;   % carrier
   t=[0:length(wave)-1]/fs;
	carrier_a=sin(2*pi*f_lo*t); 	% carrier sine;

	if_data=wave.*carrier_a;	% down converting   
      
   %************************************
   %		Design A Low-Pass Filter 		*
   %************************************
   
   f_lowpass=[   0,    65/500,   80/500 ,    1.     ];
        m=   [   1.,   1.0 ,       1e-5,       1e-5   ];

n1=80;
b=firls(n1, f_lowpass, m);

if_data=filter(b,1,if_data);   

	% save if_50_data for demod

if_50_data=if_data;				   

	%***************************
   %		Times 2					*
   %***************************
   
   if_data=if_data.*if_data;
      
   %***************************
   %		Times 4					*
   %***************************
   
   if_data=if_data.*if_data;
      
   
   %*********************************
   %		Design A Low-Pass Filter	*
   %*********************************
      b=fir1(400,[213/500 220/500]);

a=1;

if_data=filter(b,a,if_data);

	%***************************
   %		sign  					*
   %***************************
   
   if_data=sign(if_data);   
   
if_data=filter(b,a,if_data);


%***********************
% count cycle
%***********************

carrier_ref=if_data(1000:5000-100); 
length_ref=length(carrier_ref);
cycle_count=0;
for i=1:length_ref-1;
   if carrier_ref(i)*carrier_ref(i+1) <=0
      cycle_count=cycle_count+1;
   end
end

carrier_cycle=(length_ref-1)/fs/cycle_count*2;
carrier_frequency=1/carrier_cycle;

true_carrier_frequency=carrier_frequency/4;

   %*******************
	%	Generate LO	   *
   %*******************   
   length_wave=length(if_50_data);
   t=[0:length_wave-1]/fs;
   lo_cos    =cos(2*pi*true_carrier_frequency*t);
   lo_sin=cos(2*pi*true_carrier_frequency*t+pi*22.5/180);
   
   %***************************
   %		Find LO Phase			*
   %***************************
   
   
	%***************************
   %		Times 2					*
   %***************************
   
   if_data=if_50_data.^2;
   lo_cos=lo_cos.^2;
   lo_sin=lo_sin.^2;
   %***************************
   %		Times 4					*
   %***************************
   
   if_data=if_data.^2;
   lo_cos=lo_cos.^2;
   lo_sin=lo_sin.^2;

      
   
   %*********************************
   %		Design A Band-Pass Filter	*
   %*********************************
   b=fir1(400,[213/500 220/500]);
   a=1;
if_data=filter(b,a,if_data);
lo_cos=filter(b,a,lo_cos);
lo_sin=filter(b,a,lo_sin);

a_cos=if_data.*lo_cos;
a_sin=if_data.*lo_sin;

for section=1:frequency_draft;
   step=5000/frequency_draft;
   
	a_cos_a=sum(a_cos(1+(section-1)*step: section*step));
	a_sin_a=sum(a_sin(1+(section-1)*step: section*step));
	offset_theta=angle(a_cos_a+1i*a_sin_a)*180/pi;
   
   %****************
   %	Down Convert LO
   %***************
   
   ttt=(offset_theta/4+45)*pi/180;
   lo_cos_ref((1+(section-1)*step: section*step))  = ...
      cos(2*pi*true_carrier_frequency*t((1+(section-1)*step: section*step))+ttt);
   lo_sin_ref((1+(section-1)*step: section*step))	= ...
      sin(2*pi*true_carrier_frequency*t((1+(section-1)*step: section*step))+ttt);
   
end

data_length=length(lo_cos_ref);

   i_data=if_50_data(1:data_length).*lo_cos_ref;
   q_data=if_50_data(1:data_length).*lo_sin_ref;
   
   
   %*********************************
   %		Design A Low-Pass Filter	*
   %*********************************
      b=fir1(100,25/500);

i_data=filter(b,1,i_data);
q_data=filter(b,1,q_data);

%******************
%	Timing Loop		*
%******************

f_symbol=19.51e6;
n_ratio=fs/f_symbol;
iq_offset=400;
eye_max=0;

for timing_offset=0:50
   
   
	n(1)=0;
	n_sample=0;

   for i=2:4000
      
      n(i)=round( (i+timing_offset)/n_ratio);
      if n(i) > n(i-1)+.5		% take sample;
         n_sample=n_sample+1;
         i_sample(n_sample)=i_data(iq_offset+i);	
         q_sample(n_sample)=q_data(iq_offset+i);
         %[i n_sample]
         %pause
      end
      
   end
   
   eye(timing_offset+1)=sum(i_sample.^2+q_sample.^2);
	   
   if eye(timing_offset+1) > eye_max
      eye_max=eye(timing_offset+1);
      max_offset=timing_offset;
      i_sample_out=i_sample;
      q_sample_out=q_sample;
   end
end

iq_sample=i_sample_out+1i*q_sample_out;

plot(iq_sample, '*');
grid
title(' Demodulation Results');
pause;

%***************************
%	correct phase error		*
%***************************




ref_45=45*pi/180;
ref_135=135*pi/180;
phase_correction=0;
rate=1;

for kk=1:n_sample-2
   
   angle_rad=angle(iq_sample(kk));
   sign_i=sign(real(iq_sample(kk)));
   sign_q=sign(imag(iq_sample(kk)));
   
   switch sign_i*2+sign_q
      
   case  3, eye_phase_error(kk)=angle_rad-ref_45;
   case -1, eye_phase_error(kk)=angle_rad-ref_135;
   case -3, eye_phase_error(kk)=ref_135+angle_rad;
   case  1, eye_phase_error(kk)=ref_45 +angle_rad;
   otherwise '!!!! Error !!!!';
   end
   
   %phase_correction=eye_phase_error(kk)*rate+phase_correction*(1-rate);
	%phase_shift=cos(phase_correction)-1i*sin(phase_correction);
	%iq_sample(kk+1)= iq_sample(kk+1)*phase_shift*2;
   
end


%plot(eye_phase_error*180/pi, '*');
%grid
%title(' Phase Error. ');