#include <iostream.h>
#include <iomanip.h>
#include <fstream.h>
#include <math.h>
#include <stdlib.h>
#include <time.h>
#include <sys/types.h>
#include <project.h>

#define SWAP(a,b) tempr=(a);(a)=(b);(b)=tempr

class adpw
{
public:
  int reg_sig_len;
  float* reg_sig;
  int irreg_sig_len;
  float* irreg_sig;
  int* irreg_locs;
  int recon_sig_len;
  float* recon_sig;
  float* wt;
  int filt_len;
  float* filt;
  float* err;
  float max_err;
  float* Data;
  float* xd;
  
  int Done;
  int Fs;
  int Bw;
  int iter_num;
  int max_iter;
  float thresh;
  float relax;
  float total_time;
  
  adpw();
  adpw(float* rs,  int rsl,float* is,  int* il,  int isl,double  tsh,int mi,double bw,double  fs);    
  void reset(float* rs,  int rsl,float* is,  int* il,  int isl,double  tsh,int mi,double bw,double  fs);
  ~adpw();

  void four1(float data[], unsigned long int nn, int isign);
  void SetParam(int mi,float tol,float rlx);
  void SetWeights();
  void SetFilter(  int* index);
  float* GetReconstSig();
  float* GetErrorSig();
    int* GetBW(float* data,  int len,float frac);
  result  next_iter();
  void initialise();
  void start();
  friend ostream& operator << (ostream& os, adpw& p);
};

void adpw::four1(float data[], unsigned long nn, int isign)
{
	unsigned long n,mmax,m,j,istep,i;
	double wtemp,wr,wpr,wpi,wi,theta;
	float tempr,tempi;

	n=nn << 1;
	j=1;
	for (i=1;i<n;i+=2) {
		if (j > i) {
			SWAP(data[j],data[i]);
			SWAP(data[j+1],data[i+1]);
		}
		m=n >> 1;
		while (m >= 2 && j > m) {
			j -= m;
			m >>= 1;
		}
		j += m;
	}
	mmax=2;
	while (n > mmax) {
		istep=mmax << 1;
		theta=isign*(6.28318530717959/mmax);
		wtemp=sin(0.5*theta);
		wpr = -2.0*wtemp*wtemp;
		wpi=sin(theta);
		wr=1.0;
		wi=0.0;
		for (m=1;m<mmax;m+=2) {
		  for (i=m;i<=n;i+=istep) {
		    j=i+mmax;
		    tempr=wr*data[j]-wi*data[j+1];
		    tempi=wr*data[j+1]+wi*data[j];
		    data[j]=data[i]-tempr;
		    data[j+1]=data[i+1]-tempi;
		    data[i] += tempr;
		    data[i+1] += tempi;
		  }
		  wr=(wtemp=wr)*wpr-wi*wpi+wr;
		  wi=wi*wpr+wtemp*wpi+wi;
		}
		mmax=istep;
	}
}

void adpw::SetParam(int mi,float tol,float rlx)
{
    max_iter=mi;
    thresh=tol;
    relax=rlx;    
}

void adpw::SetWeights()
{
  for(  int i=1;i<irreg_sig_len-1;i++)
    wt[i]=(irreg_locs[i+1]-irreg_locs[i-1])/2.0;
  wt[irreg_sig_len-1]=(irreg_locs[0]-irreg_locs[irreg_sig_len-2]+recon_sig_len)/2.0;
  wt[0]=(irreg_locs[1]-irreg_locs[irreg_sig_len-1]+recon_sig_len)/2.0;
  //wt[irreg_sig_len-1]=-1.0;
  //wt[0]=-1.0;
}

float* adpw::GetReconstSig()
{
    return recon_sig;
}

float* adpw::GetErrorSig()
{
    return err;
}

void adpw::SetFilter(  int* index)
{
      int i;
    for(i=0;i<filt_len;i++)
	filt[i]=0.0;
    for(i=0;i<index[0];i++)
	{
	    filt[i]=1.0;
	}
  for(i=index[1];i<recon_sig_len;i++)
	{
	    filt[i]=1.0;
	}
}

result adpw::next_iter()
{
    float iter_time=0.0;
    float mean_err = 0;
    if(Done== -1)
 	{
	   
	    float err_sum,xa_sum;
	      int i;
	    clock_t time1,time2;
	    iter_num++;
	    //if((iter_num%5)==0)
	    cout << iter_num << endl;
	    for(i=0;i<recon_sig_len;i++)
		{
		    xd[i]=0.0;
		    err[i]=0.0;
		}
	    time1=clock();
	    for(i=0;i<irreg_sig_len;i++)
		{
		    xd[irreg_locs[i]]=(irreg_sig[i]-recon_sig[irreg_locs[i]])*wt[i];
		}
	    for(i=0;i<recon_sig_len;i++)
		{
		    Data[2*i]=xd[i];
		    Data[2*i+1]=0.0;
		}
	    four1(Data-1,recon_sig_len,1);
	    for(i=0;i<recon_sig_len;i++)
		{
		    Data[2*i]=Data[2*i]*filt[i];
		    Data[2*i+1]=Data[2*i+1]*filt[i];
		}
	    //cout << "Computing Inverse FFT" << endl;
	    four1(Data-1,recon_sig_len,-1);
	    for(i=0;i<recon_sig_len;i++)
		{
		    recon_sig[i]=relax*Data[2*i]/recon_sig_len+recon_sig[i];
		}
	    time2=clock();
	    iter_time=(float)(time2-time1);
	    total_time=total_time+iter_time;
	    err_sum=0;xa_sum=0;
	    for(i=0;i<reg_sig_len;i++)
		{
		    err[i]=recon_sig[i]-reg_sig[i];
		    err_sum=err_sum+err[i]*err[i];
		    xa_sum=xa_sum+reg_sig[i]*reg_sig[i];
		}
	    err_sum=sqrt(err_sum);
	    mean_err=(err_sum/sqrt(xa_sum))*100.0;
	    cout << endl << " Error wrt Original signal  = " << mean_err << endl;
	    float pct=2;
	    if(max_err!=0.0) pct=fabs(100*(mean_err-max_err)/mean_err);
	    cout << " Percentage change  = " << pct << endl;
	    if((mean_err<thresh/100.0) || (iter_num>max_iter) || (pct<.0001))
		{
		  cout << "No. of Iter = " << iter_num << endl;
		  cout << endl << " Last error = " << mean_err << endl;
		  Done=1;
		}
	    if((max_err!=0.0) && (max_err<mean_err))
		{
		  cerr << "Diverging in ADPW quitting\n";
		  Done=1;
		}
	    max_err=mean_err;
	}    
    return result(recon_sig,recon_sig_len,iter_num,iter_time,total_time,mean_err,err,Done);
}

  int* adpw::GetBW(float* data,  int len,float frac)
{
  float energy;
  float temp;
    int i;
    int *bw1=new   int[2];
  energy=0.0;
  for(i=1;i<len/2;i++)
    {
      energy=energy+data[2*i]*data[2*i];
    }
  cout << energy << endl;
  temp=0.0;
  for(i=1;i<len/2;i++)
    {
      temp=temp+data[2*i]*data[2*i];
      if(frac<=(temp/energy)) break;
    }
  cout << temp << endl;
  bw1[0]=i;
  energy=0.0;temp=0.0;
  for(i=len-1;i>=len/2;i--)
    {
      energy=energy+data[2*i]*data[2*i];
    }
  cout << energy << endl;
  temp=0.0;
  for(i=len-1;i>=len/2;i--)
    {
      temp=temp+data[2*i]*data[2*i];
      if(frac<=(temp/energy)) break;
    }
  cout << temp << endl;
  bw1[1]=i;
  return bw1;
};

void adpw::initialise()
{
      int i,*index=new   int[2];
 
    cout  << "Regular Signal Length :  " << reg_sig_len << endl
	  << "Irregular Signal Length: " << irreg_sig_len << endl
	  << "Reconstructed Signal Length: " << recon_sig_len << endl
	  << "Filter Leength: " << filt_len << endl
	  << "Sampling Frequency: " << Fs << endl
	  << "BandWidth:  " << Bw  << endl
	  << "Relaxation Parameter:  " << relax << endl
	  << "Maximum number of Iterations: " << max_iter << endl
	  << "Threshhold Error: " << thresh << endl;
	
    for(i=0;i<recon_sig_len;i++)
	{
	    Data[2*i]=reg_sig[i];
	    Data[2*i+1]=0.0;
	}
    four1(Data-1,recon_sig_len,1);    
    index=GetBW(Data,recon_sig_len,0.99);
    cout << "index =" << index[0] << " " << index[1]  << endl;
    SetWeights();
      int temp=(  int) (filt_len*Bw/Fs);
    if(temp>index[0])
	{
	    index[0]=temp;
	    index[1]=filt_len-temp;
	}
    cout << "index =" << index[0] << " " << index[1]  << endl;
    SetFilter(index);
    for(i=0;i<recon_sig_len;i++)
	{
	    recon_sig[i]=0.0;
	}
    iter_num=0;
    max_err = 0;
 }

void adpw::start()
{
    //cout << "Starting int the Iterations..." << endl;
    result ans;
    while(Done==-1)
	{
	    ans=next_iter();
	}
    cout << total_time << endl;
}

ostream& operator << (ostream& os, adpw& p)
{
    for(  int i=0;i<p.recon_sig_len;i++)
	os << i << '\t' << p.recon_sig[i] << '\n';
    return os;
};

adpw::adpw(float* rs,  int rsl,float* is,  int* il,  int isl,double  tsh,int mi,double bw,double  fs) 
 {
   reg_sig_len=rsl;
   irreg_sig_len=isl;
   reg_sig=rs;
   irreg_sig=is;
   irreg_locs=il;
   recon_sig_len=rsl;
    
    Done=-1;
    Fs=(  int)(fs);
    Bw=(  int)(bw);
    iter_num=0;
    max_iter=mi;
    relax=1.0;
    thresh=(float)(tsh);    
    filt_len=recon_sig_len;
    total_time=0.0;

    recon_sig=new float[recon_sig_len];
    wt=new float[irreg_sig_len];
    filt=new float[filt_len];
    err=new float[recon_sig_len];
    initialise();
}

adpw::adpw()
 {
   Done=-1;
   iter_num=0;
   relax=1.0;
   recon_sig=NULL;
   wt=NULL;
   filt=NULL;
   err=NULL;
   max_err = 0;
 }


void adpw::reset(float* rs,  int rsl,float* is,  int* il,  int isl,double  tsh,int mi,double bw,double  fs)
{
    reg_sig_len=rsl;
    reg_sig=rs;
    irreg_sig_len=isl;
    irreg_sig=is;
    irreg_locs=il;
    recon_sig_len=rsl;
    
    Done=-1;
    Fs=(  int)(fs);
    Bw=(  int)(bw);
    iter_num=0;
    max_iter=mi;
    relax=1.0;
    thresh=(float)(tsh);    
    filt_len=recon_sig_len;
     total_time=0.0;

    if(xd!=NULL) delete xd;
    if(Data!=NULL) delete Data;
     if(recon_sig!=NULL) delete recon_sig;
     if(wt!=NULL) delete wt;
     if(filt!=NULL) delete filt;
     if(err!=NULL) delete err;
     recon_sig=new float[recon_sig_len];
     wt=new float[irreg_sig_len];
     filt=new float[filt_len];
     err=new float[recon_sig_len];
     xd=new float[recon_sig_len];
     Data=new float[2*recon_sig_len];
     initialise();
 };

adpw::~adpw()
{
    if(recon_sig!=NULL) delete recon_sig;
    if(wt!=NULL) delete wt;
    if(filt!=NULL) delete filt;
    if(err!=NULL) delete err;
    if(xd!=NULL) delete xd;
    if(Data!=NULL) delete Data;
}