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

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

class adpw
{
 private :
    long int reg_sig_len;
    float* reg_sig;
    long int irreg_sig_len;
    float* irreg_sig;
    long int* irreg_locs;
    long int recon_sig_len;
    float* recon_sig;
    float* wt;
    long int filt_len;
    float* filt;
    float* err;
    float max_err;
    
    int Done;
    long int Fs;
    long int Bw;
    int iter_num;
    int max_iter;
    float thresh;
    float relax;
    float total_time;

 public:
    adpw(float* rs,long int rsl,float* is,long int* il,long int isl,double  tsh,int mi,double bw,double  fs);    
    void reset(float* rs,long int rsl,float* is,long int* il,long int isl,double  tsh,int mi,double bw,double  fs);
    ~adpw();

    void four1(float data[], unsigned long nn, int isign);
    void SetParam(int mi,float tol,float rlx);
    void SetWeights();
    void SetFilter(long int* index);
    float* GetReconstSig();
    float* GetErrorSig();
    long int* GetBW(float* data,long 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;
	float 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(long 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(long int* index)
{
    long int i;
    for(i=0;i<filt_len;i++)
	filt[i]=0.0;
    for(i=0;i<index[1];i++)
	{
	    filt[i]=1.0;
	}
  for(i=index[2];i<recon_sig_len;i++)
	{
	    filt[i]=1.0;
	}
}

result adpw::next_iter()
{
    float iter_time=0.0;
    if(Done==0)
 	{
	    float* xd;
	    xd=new float[recon_sig_len];
	    float* data;
	    data=new float[2*recon_sig_len];
	    float max_err;
	    static float prev_err=0;
	    float err_sum,xa_sum;
	    long 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]=reg_sig[i]-recon_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);
	    max_err=err_sum/sqrt(xa_sum);
	    cout << endl << " Error wrt Original signal  = " << max_err << endl;
	    float pct=2;
	    if(prev_err!=0.0) pct=fabs(100*(max_err-prev_err)/prev_err);
	    cout << " Percentage change  = " << pct << endl;
	    if((max_err<thresh) || (iter_num>max_iter) || (pct<.0001))
		{
		    cout << "No. of Iter = " << iter_num << endl;
		    cout << endl << " Last error = " << max_err << endl;
		    Done=1;
		}
	    if((prev_err!=0.0) && (prev_err<max_err))
		{
		    cerr << "Diverging...Therefore Terminated\n";
		    Done=1;
		}
	    prev_err=max_err;
	    delete(xd);
	    delete(data);
	}    
    return result(recon_sig,recon_sig_len,iter_num,iter_time,total_time,max_err,err,Done);
}

long int* adpw::GetBW(float* data,long int len,float frac)
{
  float energy;
  float temp;
  long int i;
  long int *bw1=new long 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[1]=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[2]=i;
  return bw1;
};

void adpw::initialise()
{
    long int i,*index=new long 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;
	
    float* data;
    data=new float[2*recon_sig_len];
    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[1] << " " << index[2]  << endl;
    SetWeights();
    long int temp=(long int) (filt_len*Bw/Fs);
    if(temp>index[1])
	{
	    index[1]=temp;
	    index[2]=filt_len-temp;
	}
    cout << "index =" << index[1] << " " << index[2]  << endl;
    SetFilter(index);
    for(i=0;i<recon_sig_len;i++)
	{
	    recon_sig[i]=0.0;
	}
    iter_num=0;
    delete(data);
 }

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

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

adpw::adpw(float* rs,long int rsl,float* is,long int* il,long 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=0;
    Fs=(long int)(fs);
    Bw=(long 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();
}

void adpw::reset(float* rs,long int rsl,float* is,long int* il,long 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=0;
    Fs=(long int)(fs);
    Bw=(long int)(bw);
    iter_num=0;
    max_iter=mi;
    relax=1.0;
    thresh=(float)(tsh);    
    filt_len=recon_sig_len;
     total_time=0.0;

     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];
    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;
}