/*This program implements Costas Loop in C++.The carrier frq is 6.9 MHz & Data
Rate is 300 KHZ . The Sampling frequency is taken as 30 MHz ( comes from the
relation Fs>=2(2Fc+DR+Fdrift) .The natural frequency of oscilation is 500Hz*/
#include<iostream.h>
#include<conio.h>
#include<math.h>
#include<stdlib.h>
#include<dos.h>
#include<stdio.h>
const float PI=3.1416;
const float E=1.0; //signal power in Watt.
//Variance of Gaussian Random variable
float q[16]={2.336797,2.082673,1.856184,1.654326,1.474420,1.314078,1.171173,1.043809,0.930296,0.829127,0.738960,0.658599,0.586977,0.523144,0.466252,0.415548};
//coefficients of Butterworth BPF having lower Fc=6.59MHz ,Higher Fc=7.21 MHz
float Prefilt_num[]= {7.055046e-14,0,-2.11651e-13,0,2.11651e-13,0,-7.055046e-14};
float Prefilt_den[]={2.923078e-10,-2.107514e-10,8.51759e-10,-3.898526e-10,7.812128e-10,-1.77213e-10,2.2540959e-10};
//coefficients of I & Q Arm LPF having Fc=0.31MHz
float Arm_num[]={3.20942e-5,9.628e-5,9.628e-5,3.20942e-5};
float Arm_den[]={1,-2.8701,2.7486,-0.878204};
//coefficients of loop filter
float Loop_num[]={7.000041666e-3,-6.999958333e-3};
float Loop_den[]={1,-1};
float Fc=6900000; //carrier frequency =6.9 MHZ
float DR=300000; //data rate =0.3 MHZ
float Fs=30000000; //Sampling frequency =30 MHz
float LF_out1=0;
main()
{
clrscr();
FILE*fp3=fopen("snr_cst.txt","w+");
FILE*fp4=fopen("ber_cst.txt","w+");
long no_of_bits=100000000; //100 MBits taken,100 samples per bit
float nob=no_of_bits;
float VCO_gain=100000; //VCO gain=100KHz/Volt
int curr_bit=0; //current bit from random binary src
randomize();
long j=0,k=0;
float curr_samp; // current sample
float Fdrift=0,Ph_drift=0,LF_out=0;//Fdrift=carrier freq drift,Ph_drift=
//carrier phase drift,LF_out=loop filter output;
float VCO_I=0,VCO_Q=0; //Inphase and Quadrature output of VCO
float I_mult=0,Q_mult=0; //Inphase and Quadrature arm multipliers
//inputs and outputs of Butterworth bpf
float bpf_in[7]={0};
float bpf_out[7]={0};
float result=0; //result sample for current bit
int decision=0;
float bpf_output=0;
//inputs and outputs of I & Q Arm lpf
float Iarmlpf_in[4]={0};
float Iarmlpf_out[4]={0};
float Qarmlpf_in[4]={0};
float Qarmlpf_out[4]={0};
float Iarm_output=0,Qarm_output=0;
//inputs and outputs of Loop filter
float Loop_in[2]={0};
float Loop_out[2]={0};
//third multiplier
float IQmult=0;
cout <<"\n enter the frequency drift in Hz ";
cin>>Fdrift;
cout <<"\n enter the Phase drift in Radian ";
cin>>Ph_drift;
int a_curr=0;
float no_of_err=0;
float N;//noise psd
float a1,b,theta,R,C;
float value[100]={0},value1=0,value2=0;
for(int i1=0;i1<16;i1++)
{
N=2*q[i1]*q[i1]; //Noise PSD
fprintf(fp3,"%f%c",10*log(E/N)/log(10)+10*log10(100/23),' ');
//as carrier frequency is 6.9 MHz and Sampling frequency is 30 MHz ,
//so the ratio is 30/6.9
no_of_err=0; //initialize no of erroneous bits to 0
int a_prev=0; //previous bit
for( j=0;j<no_of_bits*100;j++)
{
if(j%100==0)
{ a_prev=a_curr;
curr_bit=random(2);
//cout<<" "<<curr_bit<<" ";
a_curr=curr_bit;
if(curr_bit==0) //NRZ signalling
curr_bit=-1;
}
//the current sample of bpsk modulator
curr_samp=curr_bit*sin( 2*PI*(Fc+Fdrift)*(float)j/Fs +Ph_drift);
//this current sample is passed through butterworth BPF
//generate noise sample for current sample
a1=0.01*random(100);
b=0.01*random(100);
theta=2*PI*b;
R=sqrt(2*q[i1]*q[i1]*log(1/(1-a1)));
C=R*cos(theta); //noise sample
bpf_in[0]=curr_samp+C;
bpf_out[0]=0;
for(k=0;k<7;k++)
bpf_out[0]+=bpf_in[k]*Prefilt_num[k];
for(k=1;k<7;k++)
bpf_out[0]-=bpf_out[k]*Prefilt_den[k];
bpf_out[0]/=Prefilt_den[0];
bpf_output=bpf_out[0];
//shifting the coefficients
for(k=0;k<6;k++)
{ bpf_in[6-k]=bpf_in[5-k];
bpf_out[6-k]=bpf_out[5-k];
}
/*----
---------------------------------------------------------*/
LF_out1+=LF_out;
//inphase and quad signals from VCO
VCO_I=2*sin( 2*PI*(Fc)*(float)j/Fs+2*PI*VCO_gain*LF_out1/Fs);
VCO_Q=2*cos( 2*PI*(Fc)*(float)j/Fs+2*PI*VCO_gain*LF_out1/Fs);
/*-------------------------------------------------------------*/
//I & Q arm Multipliers
I_mult=bpf_output*VCO_I;
Q_mult=bpf_output*VCO_Q;
/*-------------------------------------------------------------*/
//passing the multiplier outputs through I & Q arm filters
Iarmlpf_in[0]=I_mult;
Iarmlpf_out[0]=0;
for(k=0;k<4;k++)
Iarmlpf_out[0]+=Iarmlpf_in[k]*Arm_num[k];
for(k=1;k<4;k++)
Iarmlpf_out[0]-=Iarmlpf_out[k]*Arm_den[k];
Iarmlpf_out[0]/=Arm_den[0];
Iarm_output=Iarmlpf_out[0];
for(k=0;k<3;k++)
{ Iarmlpf_in[3-k]=Iarmlpf_in[2-k];
Iarmlpf_out[3-k]=Iarmlpf_out[2-k];
}
value[(j%100)]=Iarm_output;
//taking the decision on current received bit by passing it throungh I&D
//simpson rule for integration is applied here
if((j+1)%100==0)
{
for(int l=1;l<100;l++)
{ if(!(l%2))
value1+=value[l];
else
value2+=value[l];
}
result=(value[0]+4*value2+2*value1)/3;
if(result<0)
decision=0;
else
decision=1;
value1=0;
value2=0;
if(decision!=a_prev)
no_of_err++;
}
/*----------------------------------------------------------------*/
Qarmlpf_in[0]=Q_mult;
Qarmlpf_out[0]=0;
for(k=0;k<4;k++)
Qarmlpf_out[0]+=Qarmlpf_in[k]*Arm_num[k];
for(k=1;k<4;k++)
Qarmlpf_out[0]-=Qarmlpf_out[k]*Arm_den[k];
Qarmlpf_out[0]/=Arm_den[0];
Qarm_output=Qarmlpf_out[0];
//shifting the coefficients
for(k=0;k<3;k++)
{ Qarmlpf_in[3-k]=Qarmlpf_in[2-k];
Qarmlpf_out[3-k]=Qarmlpf_out[2-k];
}
/*----------------------------------------------------------------*/
//multiplying the outputs of I & Q arm Filters
IQmult=Iarm_output*Qarm_output;
//passing the multiplier through Loop filter
Loop_in[0]=IQmult;
Loop_out[0]=0;
for(k=0;k<2;k++)
Loop_out[0]+=Loop_in[k]*Loop_num[k];
Loop_out[0]-=Loop_out[1]*Loop_den[1];
Loop_out[0]/=Loop_den[0];
LF_out=Loop_out[0];
//shifting the coefficients
Loop_out[1]=Loop_out[0];
Loop_in[1]=Loop_in[0];
}
//verify the result
fprintf(fp4,"%f%c",(no_of_err/nob),' '); cout<<" "<<(no_of_err/nob)<<endl;
}
cout<<"\n\n\n done ";
fclose(fp3);
fclose(fp4);
getch();
return 0;
}
