//File: setieasy_v133b.cpp
//Source by:
//  Robert L. Samuell IV
//  Robert A. Lodder Ph.D
//  Dan Kennedy
//  Kah-Siew Ho
//                                             
//Software Name: SETIEasy
//Version: 1.33b
//Last Modified: Friday, March 30, 2001

#include<fstream.h>
#include<math.h>
#include<stdio.h>
#include<stdlib.h>
#include<iomanip.h>
#include<time.h>
#include<errno.h>
#include<string>
#define SWAP(a,b) tempr=(a); (a)=(b); (b)=tempr

class Statpack1{
      public:

        //CONSTRUCTOR
        Statpack1();

        //FUNCTIONS FOR FILES
        char* string_to_Cstring(string s1);
        string Cstring_to_string(char* s1);
        char* get_input_filename(int num);
        char* get_outData_filename();
        char* get_grafData_filename();
	char* get_factor();
	void set_input_filename(char* fname1, char* fname2);
        void open_input_file();
        void reopen_input_file();
        unsigned long get_nChunkSize(FILE* &infile);
        void open_output_file(char *outData, char *graphData);
        void close_input_file();
        void close_output_file() { fout.close(); gout.close(); }
	void print_bit();
        void put_output_file_header();
        void put_output_file_footer(time_t *TIME);
        void put_dechirp_factor_header(char* dechirp_factor, time_t *TIME);
        void print_fftblock_header();
        void print_one_fft_header();
        void hello();

        //FUNCTIONS FOR ANALYSIS
        void initialize_rand();
	void get_data_length(unsigned long n);
	void initialize_variables();
        void initialize_data_set1();
	void initialize_data_set2();
        void delete_data_set1(){ delete [] data1; }
	void delete_data_set2(){ delete [] data2; }
        void initialize_train_set();
        void delete_train_set();
        void initialize_test_set();
        void delete_test_set();
        void set_shadow_dimensions();
        void dfft(float data[],unsigned long n,int isign);
        void sizer_ft(float data[], unsigned long nn, int isign);
        float fsob(float **b_train, float *c_train,
                   float **b_test, float *c_test);
        void sortarray(float array[], unsigned long n);
        void limits();
        void find_confid();
        void find_coeff();
        float sdev(float *array, float mean, unsigned long n);
        float sumcols(float **a, unsigned long norows, unsigned long colnumber);
        void make_train_set();
        void make_test_set();
        void interpret_data();
        void init_start_sample();
        bool is_sample();
        void next_sample();
        void put_final_analysis();

	void make_dechirped_file(char* dechirpFactor);
	unsigned long open_dechirped_file1();
	unsigned long open_dechirped_file2();
	unsigned long make_dechirped_file(unsigned long nChunkSize, FILE* &infile, 
					  FILE* &outfile);
	void close_dechirped_file1();
	void close_dechirped_file2();


      private:
	int unusual_signal;
        bool signal;
	double Factor;
	char *D_factor;
	string DFactor;
        char *in_file1, *in_file2, *out_file1, *out_file2;
        string vers;
        string filename1, filename2;
        FILE *infile1, *infile2, *outfile1, *outfile2;
        ofstream fout, gout;
        unsigned long b;
        unsigned long rows, cols;
        unsigned long start_sample;
        unsigned long end_sample;
        unsigned long data_len;
        float rc;
        float confid;
        float *data1, *data2;
        float **btn;
        float *ctn;
        float **btest;
        float *ctest;
	short int nBitsPerSample1, nBitsPerSample2;
};

Statpack1::Statpack1(){
     b=100;
     cols=8192;
     start_sample=0;
     rc=0.0;
     confid=0.0;
     signal = false;
     unusual_signal = 0;
     vers = "v1.33b beta";
}

void Statpack1::hello(){
     cout << "SETIEASY " << vers << " - Unusual Signal Detector" << endl;
     cout << "Copyright (c) 2000 University of Kentucky" << endl;
}

char* Statpack1::string_to_Cstring(string s1){
       char* s2 = new char[s1.length()+1];

       for(int i=0; i<s1.length(); i++)
          s2[i] = s1[i];

       s2[s1.length()] = '\0';
       return s2;
}

string Statpack1::Cstring_to_string(char* s1){
      string s2;

      for(int i=0; s1[i] != '\0'; i++)
         s2 += s1[i];
      return s2;
}

char* Statpack1::get_input_filename(int num){
      string temp;
      cout << "What is the name of the input file " << num << "?" << endl;
      cin >> temp;
      return string_to_Cstring(temp);
}

char* Statpack1::get_outData_filename(){
      string temp;
      cout << "What do you want to name the outData file?" << endl;
      cin >> temp;
      return string_to_Cstring(temp);
}

char* Statpack1::get_grafData_filename(){
      string temp;
      cout << "What do you want to name the grafData file?" << endl;
      cin >> temp;
      return string_to_Cstring(temp);
}

char* Statpack1::get_factor(){
      string temp;
      cout << "What do you want to use as the dechirp factor?" << endl;
      cin >> temp;
      return string_to_Cstring(temp);
}

void Statpack1::initialize_variables(){
      signal = false;
      unusual_signal = 0;
      start_sample=0;
}

unsigned long Statpack1::get_nChunkSize(FILE* &infile){
     /* Scan in the 4 byte integer */
     unsigned char ch;
     unsigned long nChunkSize = 0;

	for(int i=0; i<4; i++){
	   fscanf( infile, "%c", &ch);
	   nChunkSize = (unsigned long) (nChunkSize + ch*pow(2, i*8));
	}
	return nChunkSize;
}

void Statpack1::set_input_filename(char* fname1, char* fname2){
     in_file1 = fname1;
     filename1 = Cstring_to_string(in_file1);

     in_file2 = fname2;
     filename2 = Cstring_to_string(in_file2);
}

void Statpack1::open_input_file(){
     unsigned char ch;

     if( (infile1 = fopen( out_file1, "rb" )) == NULL ) {
       printf( "The source file %s was not opened\n", out_file1 );
       exit(1);
     }
     else{
       printf( "The source file %s was opened\n", out_file1 );
     }

     if( (infile2 = fopen( out_file2, "rb" )) == NULL ) {
       printf( "The source file %s was not opened\n", out_file2 );
       exit(1);    
     }
     else{
       printf( "The source file %s was opened\n", out_file2 );
     }

     for (int i=0; i<44; i++)
	  fscanf( infile1, "%c", &ch );

     for (int j=0; j<44; j++)
	  fscanf( infile2, "%c", &ch );
}

void Statpack1::close_input_file(){
     int Status;

     Status = fclose(infile1);
     if(Status == 0) printf("%s was closed\n", out_file1);
     else printf("%s did not close\n", out_file1);

     Status = fclose(infile2);
     if(Status == 0) printf("%s was closed\n", out_file2);
     else printf("%s did not close\n", out_file2);
}

void Statpack1::open_output_file(char *outData, char* graphData){
      char *out_data;
      char *graph_data;
      string temp1, temp2;
      string D = "D";
      string underscore = "_";

      temp1 = D;
      temp1 += DFactor;
      temp1 += underscore;
      temp1 += Cstring_to_string(outData);
      out_data = string_to_Cstring(temp1);

      temp2 = D;
      temp2 += DFactor;
      temp2 += underscore;
      temp2 += Cstring_to_string(graphData);
      graph_data = string_to_Cstring(temp2);

     fout.open(out_data);
     if(fout.fail()){
       cout << "Open outData file failed!" << endl;
       exit(1);
     }
     gout.open(graph_data);
     if(gout.fail()){
       cout << "Open graphData file failed!" << endl;
       exit(1);
     }
}

void Statpack1::put_output_file_header(){
     fout << "SETIEASY " << vers << " analysis of " << filename1 << " and "
     << filename2 << endl;
}

void Statpack1::put_output_file_footer(time_t *TIME){
     fout << ctime(TIME);
}

void Statpack1::put_dechirp_factor_header(char* dechirp_factor, time_t *TIME){
     fout << ctime(TIME);
     fout << "@~~~~~ Dechirp Factor: " << dechirp_factor << " ~~~~~@" << endl;
}

void Statpack1::initialize_rand(){
     cout << "Initializing random number generator" << endl;
     srand((unsigned int) time(NULL));
}

void Statpack1::get_data_length(unsigned long n){
     data_len = (unsigned long) (pow(2,n));
}

void Statpack1::print_bit(){
    if (nBitsPerSample1 == 8) {
	cout << in_file1 << " and " << in_file2 << " are 8-bit wav files." << endl;
	fout << in_file1 << " and " << in_file2 << " are 8-bit wav files." << endl;
    }
    else if (nBitsPerSample1 == 16) {
	cout << in_file1 << " and " << in_file2 << " are 16-bit wav files." << endl;
	fout << in_file1 << " and " << in_file2 << " are 16-bit wav files." << endl;
    }
}

void Statpack1::initialize_data_set1(){
    float sampleVal;
    unsigned char ch;

    cout << "Reading data set 1" << endl;
    data1 = new float[(2*data_len)+1];

    if (nBitsPerSample1 == 8){
        for(unsigned long index = 0; (feof(outfile1) == 0) && (index < data_len);
	    index++) {

            /* Read in data array */
            fscanf( infile1, "%1c", &ch );
            data1[index] = ((float)ch-128)/128;
        }
    }
    if (nBitsPerSample1 == 16){
        for(unsigned long index = 0; (feof(outfile1) == 0) && (index < data_len);
	    index++) {

            /* Read in the low byte */
            fscanf( infile1, "%1c", &ch );
            sampleVal = (float)ch;

            /* Read in the high byte. */
            fscanf( infile1, "%1c", &ch );
            sampleVal = sampleVal + 256 * (float) ch;
            data1[index] = (sampleVal-32768)/32768;
        }
    }
}

void Statpack1::initialize_data_set2(){
    float sampleVal;
    unsigned char ch;

    cout << "Reading data set 2" << endl;
    data2 = new float[(2*data_len)+1];

    if (nBitsPerSample2 == 8){
        for(unsigned long index = 0; (feof(outfile2) == 0) && (index < data_len);
	    index++) {

            /* Read in data array */
            fscanf( infile2, "%1c", &ch );
            data2[index] = ((float)ch-128)/128;
        }
    }
    if (nBitsPerSample2 == 16){
        for(unsigned long index = 0; (feof(outfile2) == 0) && (index < data_len);
	    index++) {

            /* Read in the low byte */
            fscanf( infile2, "%1c", &ch );
            sampleVal = (float)ch;

            /* Read in the high byte. */
            fscanf( infile2, "%1c", &ch );
            sampleVal = sampleVal + 256 * (float) ch;
            data2[index] = (sampleVal-32768)/32768;
        }
    }
}

void Statpack1::print_fftblock_header(){
    int mins;
    float sec;
    int sixty = 60;
    int samps_per_min = 480000;
    float fsamps_per_min = 480000.0;

   mins = data_len/samps_per_min;
   sec = ((data_len/fsamps_per_min)-float(data_len/samps_per_min))*sixty;

    cout << "Looking for signals " << mins << " mins and " << sec
	 << " secs long." << endl;
    fout << "Signals for " << mins << " mins and " << sec
	 << " secs long." << endl;
}

void Statpack1::print_one_fft_header(){
    int mins;
    float sec;
    int sixty = 60;
    int samps_per_min = 480000;
    float fsamps_per_min = 480000.0;

    mins = start_sample /samps_per_min;
    sec = ((start_sample/fsamps_per_min)-float(start_sample/samps_per_min))*sixty;

    cout << "End with " << mins << " mins and " << sec << " secs" << endl;
    fout << "End with " << mins << " mins and " << sec << " secs" << endl;
}

void Statpack1::initialize_train_set(){
     btn = new float*[b];
     for(unsigned long i=0; i<b; ++i) *(btn+i) = new float[cols+1];
     ctn = new float[cols+1];
}

void Statpack1::delete_train_set(){
     for(unsigned long i=0; i<b; ++i) delete [] *(btn+i);
     delete [] btn;
     delete [] ctn;
}

void Statpack1::initialize_test_set(){
     btest = new float*[b];
     for(unsigned long i=0; i<b; ++i) *(btest+i) = new float[cols+1];
     ctest = new float[cols+1];
}

void Statpack1::delete_test_set(){
     for(unsigned long i=0; i<b; ++i) delete [] *(btest+i);
     delete [] btest;
     delete [] ctest;
}

void Statpack1::set_shadow_dimensions(){
     rows = data_len/cols;
}

void Statpack1::dfft(float data[],unsigned long n,int isign){
     cout << "Calculating fast fourier transform\n";
     unsigned long i=0,init_one=0,init_too=0,init_temp=0,closer_init=0;
     unsigned long np_complete=0,start_flag=0;
     float end_flag=0.5,z_stack=0,pop_1=0,y_stack=0,dummy_y=0,counter_3y=0;
     double work=0,work_init=0,w_progress=0,pop_i=0,wtemp=0,angles=0;
     angles=3.141592653589793/(double) (n>>1);
     if(isign == 1){
       z_stack = -0.5;
       sizer_ft(data,n>>1,1);
     }
     else{
       z_stack = 0.5;
       angles = -angles;
     }
     wtemp = sin(0.5*angles);
     w_progress = -2.0*wtemp*wtemp;
     pop_i = sin(angles);
     work = 1.0+w_progress;
     work_init = pop_i;
     np_complete = n+3;
     for(i=2;i<=(n>>2);i++){
        closer_init=1+(init_temp=np_complete-(init_too=(1+(init_one=i+i-1))));
        pop_1=end_flag*(data[init_one]+data[init_temp]);
        y_stack=end_flag*(data[init_too]-data[closer_init]);
        dummy_y= -z_stack*(data[init_too]+data[closer_init]);
        counter_3y=z_stack*(data[init_one]-data[init_temp]);
        data[init_one]=pop_1+work*dummy_y-work_init*counter_3y;
        data[init_too]=y_stack+work*counter_3y+work_init*dummy_y;
        data[init_temp]=pop_1-work*dummy_y+work_init*counter_3y;
        data[closer_init]= -y_stack+work*counter_3y+work_init*dummy_y;
        work=(wtemp=work)*w_progress-work_init*pop_i+work;
        work_init=work_init*w_progress+wtemp*pop_i+work_init;
     }
     if(isign == 1){
       data[1]=((pop_1=data[1])+data[2]);
       data[2]=pop_1-data[2];
       float dtemp=data[2];
       for(start_flag=3;start_flag<=n;start_flag++) data[start_flag-1]=data[start_flag];
       data[start_flag-1]=dtemp;
     }
     else{
       data[1]=end_flag*((pop_1=data[1])+data[2]);
       data[2]=end_flag*(pop_1-data[2]);
       sizer_ft(data,n>>1,-1);
     }
}
void Statpack1::sizer_ft(float data[],unsigned long nn,int isign){
     unsigned long n=0,g_max=0,m=0,j=0,istep=0,i=0;
     float work=0,work_init=0,wtemp=0,w_progress=0,pop_i=0,angles=0;
     float tempr=0,tempi=0;
     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;
     }
     g_max=2;
     while(n>g_max){
       istep = g_max << 1;
       angles = isign*(6.2831853071796/g_max);
       wtemp = sin(0.5*angles);
       w_progress = -2.0*wtemp*wtemp;
       pop_i = sin(angles);
       work = 1.0;
       work_init = 0.0;
       for(m=1;m<g_max;m+=2){
         for(i=m;i<=n;i+=istep){
           j=i+g_max;
           tempr=work*data[j]-work_init*data[j+1];
           tempi=work*data[j+1]+work_init*data[j];
           data[j]=data[i]-tempr;
           data[j+1]=data[i+1]-tempi;
           data[i] += tempr;
           data[i+1] += tempi;
         }
         work=(wtemp=work)*w_progress-work_init*pop_i+work;
         work_init=work_init*w_progress+wtemp*pop_i+work_init;
       }
       g_max=istep;
     }
}

float Statpack1::fsob(float **b_train,float *c_train,
      float **b_test,float *c_test){
     unsigned long j,n,k;
     unsigned long m = b*2;
     float r=0.0;
     float xsum=0.0,ysum=0.0;
     float xysum=0.0;
     float xsqsum=0.0,ysqsum=0.0;
     float nnnum=0.0,dddenom=0.0;
     unsigned long d = cols;
     unsigned long cdf_len = b*2;
     float ntest = sqrt(b);
     float ntn = sqrt(b);
     float *ecdf;
     float *tcdf;
     float droot = 4;
     ecdf = new float[cdf_len+1];
     tcdf = new float[cdf_len+1];

     for(n=0;n<b;++n){
       for(k=0;k<d;++k){
          ecdf[n] += pow(fabs(b_train[n][k]-c_train[k]), droot)*ntn;
          ecdf[n+b] += pow(fabs(b_test[n][k]-c_train[k]), droot)*ntest;
       }
       ecdf[n] = sqrt(sqrt(ecdf[n]));
       ecdf[n+b] = sqrt(sqrt(ecdf[n+b]));
       tcdf[n]=ecdf[n];
       tcdf[n+b]=ecdf[n];
     }

     sortarray(ecdf-1, cdf_len);
     sortarray(tcdf-1, cdf_len);

     unsigned long templen = ((m-(m/10))-((m/10)+1))+1;

     for(j=0; j<templen; ++j){
        xsum += ecdf[j];
        ysum += tcdf[j];
        xysum += ecdf[j]*tcdf[j];
        xsqsum += ecdf[j]*ecdf[j];
        ysqsum += tcdf[j]*tcdf[j];
     }

     nnnum=(templen*xysum)-(xsum*ysum);
     dddenom=sqrt((templen*xsqsum)-(xsum*xsum));
     dddenom=dddenom*sqrt((templen*ysqsum)-(ysum*ysum));
     r=nnnum/dddenom;

     delete [] ecdf;
     delete [] tcdf;

     return r;
}

void Statpack1::sortarray(float array[], unsigned long n){
     float d_nodes;
     unsigned long i, indx_node, j, D_length;

     if(n<2) return;
     D_length=(n>>1)+1;
     indx_node=n;
     for(;;){
        if(D_length>1){
          d_nodes=array[--D_length];
        }
        else{
          d_nodes=array[indx_node];
          array[indx_node]=array[1];
          if(--indx_node == 1){
            array[1]=d_nodes;
            break;
          }
        }
        i=D_length;
        j=D_length+D_length;
        while(j<=indx_node){
           if(j < indx_node && array[j] < array[j+1]) j++;
           if(d_nodes < array[j]){
             array[i]=array[j];
             i=j;
             j <<= 1;
           }else break;
        }
        array[i]=d_nodes;
     }
}

void Statpack1::limits(){
     float sd=0.0;
     float r=0.0, mean=0.0;
     unsigned long ncors=b;
     unsigned long i,j,m,n,p,s,pick;
     float *cors;
     float **check;
     float *chtest;

     cors = new float[ncors+1];
     check = new float*[b];
     for(i=0;i<b;++i) *(check+i)=new float[cols+1];
     chtest = new float[cols+1];

     cout << "Preparing for calculation of confidence limits:" << endl;

     for(n=0; n<ncors; ++n){
       cout << "Limit " << n << " of " << ncors << ": ";
       for(j=0; j<b; ++j){
         for(m=0; m<cols; ++m){
           pick = rand() % b;
           check[j][m] = btn[pick][m];
         }
       }

       for(p=0; p<cols; ++p)
          chtest[p]=sumcols(check,b,p)/b;
       r = fsob(btn,ctn,check,chtest);
       cors[n] = r;
       cout << cors[n] << endl;
     }

     for(i=0; i<b; ++i) delete [] *(check+i);
     delete [] check;
     delete [] chtest;

     for(s=0; s<ncors; ++s) mean += cors[s];
     mean = mean/ncors;
     sd = sdev(cors,mean,ncors);
     confid = mean-(10.577*sd);

     cout << "Confidence limit: " << confid << endl;

     delete [] cors;
}

void Statpack1::find_confid(){
     dfft(data1-1, data_len, 1);
     make_train_set();
     delete_data_set1();
     limits();
}

void Statpack1::find_coeff(){
     dfft(data2-1, data_len, 1);
     make_test_set();
     delete_data_set2();
     rc = fsob(btn,ctn,btest,ctest);
}

float Statpack1::sdev(float *array,float mean, unsigned long n){
     cout << "Calculating standard deviation" << endl;
     float standdev=0.0,sum=0.0,s=0.0;
     unsigned long i;
     for(i=0;i<n;++i){
       s=(array[i]-mean)*(array[i]-mean);
       sum += s;
     }
     sum = sum/(n-1);
     standdev = sqrt(sum);
     return standdev;
}

float Statpack1::sumcols(float **a,unsigned long norows,
     unsigned long colnumber){
     unsigned long i;
     float sum;
     sum = 0.0;
     for(i=0; i<norows; i++) sum += a[i][colnumber];
     return sum;
}

void Statpack1::make_train_set(){
     cout << "Making replicates for training set" << endl;
     unsigned long i,j,k,pick,bsamp;

     for(bsamp=0; bsamp<b; ++bsamp){
        for(k=0; k<cols; ++k){
              btn[bsamp][k] = 0.0;
        }
        for(i=0; i<rows; ++i){
           pick = rand() % rows;
           for(j=0; j<cols; ++j){
              btn[bsamp][j] += fabs(data1[(pick*cols)+j]);
           }
        }
        for(k=0;k<cols;++k){
              btn[bsamp][k] = btn[bsamp][k]/rows;
        }
     }
     for(i=0;i<cols;++i) ctn[i]=sumcols(btn,b,i)/b;
}

void Statpack1::make_test_set(){
     unsigned long i,j,k,pick,bsamp;

     cout << "Making replicates for test set\n";

     for(bsamp=0; bsamp<b; ++bsamp){
        for(k=0; k<cols; ++k){
              btest[bsamp][k] = 0.0;
        }
        for(i=0; i<rows; ++i){
           pick = rand() % rows;
           for(j=0; j<cols; ++j){
              btest[bsamp][j] += fabs(data2[(pick*cols)+j]);
           }
        }
        for(k=0; k<cols; ++k){
              btest[bsamp][k] = btest[bsamp][k]/rows;
        }
     }
     for(i=0; i<cols; ++i) ctest[i]=sumcols(btest,b,i)/b;
}

void Statpack1::interpret_data(){
     cout << "rc = " << rc << "  ";
     cout << "confid = " << confid << endl;
     fout << "rc = " << rc << "  ";
     fout << "confid = " << confid << endl;
     gout << rc << endl;
     gout << confid << endl;

     if(rc < confid){
       fout << "Unusual signal detected - report to group via http://hendrix.pharm.uky.edu/seti/datafiles.html\n";
       cout << "Unusual signal detected - report to group via http://hendrix.pharm.uky.edu/seti/datafiles.html\n";
       signal = true;
       unusual_signal++;
     }
     else{
       fout << "No signals detected" << endl;
       cout << "No signals detected" << endl;
     }
}

void Statpack1::init_start_sample(){
     start_sample = data_len;
}

bool Statpack1::is_sample(){
     return ( (start_sample+data_len) < end_sample );
}

void Statpack1::next_sample(){
     start_sample += data_len;
}

void Statpack1::put_final_analysis(){
     if(signal){
	cout << "NUMBER OF UNUSUAL SIGNALS DETECTED BY SETIEASY " << vers << ": " 
	     << unusual_signal << endl << endl;
	fout << "NUMBER OF UNUSUAL SIGNALS DETECTED BY SETIEASY " << vers << ": " 
	     << unusual_signal << endl << endl;
     }
     else {
	cout << "NO UNUSUAL SIGNAL DETECTED BY SETIEASY " << vers << endl << endl;
	fout << "NO UNUSUAL SIGNAL DETECTED BY SETIEASY " << vers << endl << endl;
     }
}


unsigned long Statpack1::open_dechirped_file1(){
      char prefix[4];
      char fileFormat[8];
      char ckID[4];
      unsigned long nChunkSize, nChunkSize1;
      unsigned char ch;
      short wFormatTag;
      short nChannels;
      unsigned long nSamplesPerSecond;
      unsigned long nAvgBytesPerSecond;
      short nBlockAlign;

      string temp1;
      string D = "D";
      string underscore = "_";

/* Open source for binary read (will fail if file does not exist) */
      if( (infile1  = fopen(in_file1, "rb") ) == NULL ) {
	printf( "The source file %s was not opened\n", in_file1 );
	exit(1);
      }
      else 
	printf( "The source file %s was opened\n", in_file1 );

      temp1 = D;
      temp1 += DFactor;
      temp1 += underscore;
      temp1 += filename1;
      out_file1 = string_to_Cstring(temp1);

/* Open output for write */
      if( (outfile1 = fopen( out_file1, "w" )) == NULL ) {
	printf( "The output file %s was not opened\n", out_file1 );
	exit(1);
      }
      else
	printf( "The output file %s was opened\n", out_file1 );

     fscanf( infile1, "%4c", prefix );
     fprintf( outfile1, "%4c", prefix );
     fscanf( infile1, "%4c", &nChunkSize );
     fprintf( outfile1, "%4c", nChunkSize );
     fscanf( infile1, "%8c", fileFormat );
     fprintf( outfile1, "%8c", fileFormat );
     fscanf( infile1, "%4c", &nChunkSize );
     fprintf( outfile1, "%4c", nChunkSize );
     fscanf( infile1, "%2c", &wFormatTag );
     fprintf( outfile1, "%2c", wFormatTag );
     fscanf( infile1, "%2c", &nChannels );
     fprintf( outfile1, "%2c", nChannels );
     fscanf( infile1, "%4c", &nSamplesPerSecond );
     fprintf( outfile1, "%4c", nSamplesPerSecond );
     fscanf( infile1, "%4c", &nAvgBytesPerSecond );
     fprintf( outfile1, "%4c", nAvgBytesPerSecond );
     fscanf( infile1, "%2c", &nBlockAlign );
     fprintf( outfile1, "%2c", nBlockAlign );
     fscanf( infile1, "%1c", &ch );
     fprintf( outfile1, "%1c", ch );
     nBitsPerSample1 = (short int)ch;
     fscanf( infile1, "%1c", &ch );
     fprintf( outfile1, "%1c", ch );
     fscanf( infile1, "%4s", &ckID );
     fprintf( outfile1, "%4s", ckID );
     nChunkSize1 = get_nChunkSize(infile1);

     return nChunkSize1;
}

unsigned long Statpack1::open_dechirped_file2(){
      char prefix[4];
      char fileFormat[8];
      char ckID[4];
      unsigned long nChunkSize, nChunkSize2;
      unsigned char ch;
      short wFormatTag;
      short nChannels;
      unsigned long nSamplesPerSecond;
      unsigned long nAvgBytesPerSecond;
      short nBlockAlign;

      string temp2;
      string D = "D";
      string underscore = "_";

/* Open source for binary read (will fail if file does not exist) */
      if( (infile2  = fopen(in_file2, "rb") ) == NULL ) {
	printf( "The source file %s was not opened\n", in_file2 );
	exit(1);
      }
      else 
	printf( "The source file %s was opened\n", in_file2 );

      temp2 = D;
      temp2 += DFactor;
      temp2 += underscore;
      temp2 += filename2;
      out_file2 = string_to_Cstring(temp2);

/* Open output for write */
     if( (outfile2 = fopen( out_file2, "w" )) == NULL ) {
	printf( "The output file %s was not opened\n", out_file2 );
	exit(1);
     }
     else
	printf( "The output file %s was opened\n", out_file2 );

     fscanf( infile2, "%4c", prefix );
     fprintf( outfile2, "%4c", prefix );
     fscanf( infile2, "%4c", &nChunkSize );
     fprintf( outfile2, "%4c", nChunkSize );
     fscanf( infile2, "%8c", fileFormat );
     fprintf( outfile2, "%8c", fileFormat );
     fscanf( infile2, "%4c", &nChunkSize );
     fprintf( outfile2, "%4c", nChunkSize );
     fscanf( infile2, "%2c", &wFormatTag );
     fprintf( outfile2, "%2c", wFormatTag );
     fscanf( infile2, "%2c", &nChannels );
     fprintf( outfile2, "%2c", nChannels );
     fscanf( infile2, "%4c", &nSamplesPerSecond );
     fprintf( outfile2, "%4c", nSamplesPerSecond );
     fscanf( infile2, "%4c", &nAvgBytesPerSecond );
     fprintf( outfile2, "%4c", nAvgBytesPerSecond );
     fscanf( infile2, "%2c", &nBlockAlign );
     fprintf( outfile2, "%2c", nBlockAlign );
     fscanf( infile2, "%1c", &ch );
     fprintf( outfile2, "%1c", ch );
     nBitsPerSample2 = (short int)ch;
     fscanf( infile2, "%1c", &ch );
     fprintf( outfile2, "%1c", ch );
     fscanf( infile2, "%4s", &ckID );
     fprintf( outfile2, "%4s", ckID );
     nChunkSize2 = get_nChunkSize(infile2);

     return nChunkSize2;
}

void Statpack1::make_dechirped_file(char* dechirpFactor){
     unsigned long end_sample1, end_sample2;
     unsigned long nChunkSize1, nChunkSize2;
     char *stopstring;

     Factor = strtod(dechirpFactor, &stopstring);
     if (Factor>1.0 || Factor<-1.0) {
	 cerr << "Dechirp factor must be in the range of [-1, 1]" << endl;
	 exit(1);
     }

     D_factor = dechirpFactor;
     DFactor = Cstring_to_string(dechirpFactor);

     nChunkSize1 = open_dechirped_file1();
     nChunkSize2 = open_dechirped_file2();

//To make sure both files have same nChunksize and nBitsPerSample
   if ( nBitsPerSample1 != nBitsPerSample2 ){
	cerr << "Error, the two input files are not in the same bit!" << endl;
	exit(1);
     }

   if ( nChunkSize1 != nChunkSize2 ){
	cerr << "WARNING: The two input files are not the same length!" << endl;
   	cout << nChunkSize1 << ", " << nChunkSize2 << endl;
	exit(1);
   }

     cout << "Dechirping " << filename1 << endl;
     end_sample1 = make_dechirped_file(nChunkSize1, infile1, outfile1);
     close_dechirped_file1();

     cout << "Dechirping " << filename2 << endl;
     end_sample2 = make_dechirped_file(nChunkSize2, infile2, outfile2);
     close_dechirped_file2();

   if ( end_sample1 != end_sample2 ){
	cerr << "WARNING: Something gone wrong with the algorithm!" << endl;
   	cout << end_sample1 << ", " << end_sample2 << endl;
	exit(1);
   }
   else 
	end_sample = end_sample1;
}

unsigned long Statpack1::make_dechirped_file(unsigned long nChunkSize, FILE* &infile, 
					     FILE* &outfile){
     unsigned long byteCount, oldbyteCount;
     unsigned char ch, ch1, ch2;
     double slideFactor;


     oldbyteCount = 0;

     if (nBitsPerSample1 == 8){
	for(unsigned long i=0; i<nChunkSize; i++) {
	    slideFactor = (i/(double)nChunkSize)*Factor;
	    byteCount = (unsigned long) (i*(1+slideFactor));
	    fscanf( infile, "%1c", &ch );
	    if (byteCount>oldbyteCount) {
		do {
		   fprintf(outfile,"%c",ch);
		   oldbyteCount = oldbyteCount+1;
		} while (byteCount > oldbyteCount);
	    }
	}
     }
     else if (nBitsPerSample1 == 16){
	for(unsigned long i=0; i<(nChunkSize/2); i++) {
	    slideFactor = (i/(double)nChunkSize)*Factor;
	    byteCount = (unsigned long) (i*(1+slideFactor));
	    fscanf( infile, "%1c", &ch1 );
	    fscanf( infile, "%1c", &ch2 );
	    if (byteCount>oldbyteCount) {
		do {
		   fprintf(outfile,"%c",ch1);
		   fprintf(outfile,"%c",ch2);
		   oldbyteCount = oldbyteCount+1;
		} while (byteCount > oldbyteCount);
	    }
	}
     }

     return oldbyteCount;
}

void Statpack1::close_dechirped_file1(){
     int Status;

     Status = fclose(infile1);
     if(Status == 0) printf("%s was closed\n", in_file1);
     else printf("%s did not close\n", in_file1);

     Status = fclose(outfile1);
     if(Status == 0) printf("%s was closed\n", out_file1);
     else printf("%s did not close\n", out_file1);
}

void Statpack1::close_dechirped_file2(){
     int Status;

     Status = fclose(infile2);
     if(Status == 0) printf("%s was closed\n", in_file2);
     else printf("%s did not close\n", in_file2);

     Status = fclose(outfile2);
     if(Status == 0) printf("%s was closed\n", out_file2);
     else printf("%s did not close\n", out_file2);
}

void analyze(Statpack1&, char*, char*, char*, char*, int, char**);

void main(int argc, char* argv[]){
     time_t *start_time, *end_time;
     int x;
     int y=0;
     int no_of_factors = 0;
     int no_of_commands = argc;
     char *dechirpFactor[300];
     char *input_file1, *input_file2;
     char *outData_file;
     char *grafData_file;
     Statpack1 seti;

     seti.hello();

  if (no_of_commands <= 6){
     switch(no_of_commands){
       case 6: dechirpFactor[0] = argv[5];
	       no_of_factors = 1;
       case 5: grafData_file = argv[4];
       case 4: outData_file = argv[3];
       case 3: input_file2 = argv[2];
       case 2: input_file1 = argv[1];
     }
     switch(no_of_commands){
       case 1: input_file1 = seti.get_input_filename(1);
       case 2: input_file2 = seti.get_input_filename(2);
       case 3: outData_file = seti.get_outData_filename();
       case 4: grafData_file = seti.get_grafData_filename();
       case 5: dechirpFactor[0] = seti.get_factor();
	       no_of_factors = 1;
       case 6:
	analyze(seti, input_file1, input_file2, outData_file, grafData_file, 
	no_of_factors, dechirpFactor);
     }
  }
  else{
        input_file1 = argv[1];
        input_file2 = argv[2];
        outData_file = argv[3];
        grafData_file = argv[4];
        for ( x=5; x<no_of_commands; x++){
	    dechirpFactor[y] = argv[x];
	    no_of_factors++;
	    y++;
        }

	analyze(seti, input_file1, input_file2, outData_file, grafData_file, 
	no_of_factors, dechirpFactor);
  }
     cout << "@~~~ Program terminated ~~~~@" << endl;
}

void analyze(Statpack1& seti, char* input_file1, char* input_file2,
             char* outData_file, char* grafData_file, int num_factors, char** factor){
     time_t *TIME;
     unsigned long n;
     char* dechirp_factor;

     seti.initialize_rand();
     seti.set_input_filename(input_file1, input_file2);

  for(int x=0; x<num_factors; x++){
     dechirp_factor = factor[x];

     seti.make_dechirped_file(dechirp_factor);
     seti.open_output_file(outData_file, grafData_file);
     seti.put_output_file_header();
     seti.print_bit();

     time(TIME);
     seti.put_dechirp_factor_header(dechirp_factor, TIME);
     seti.open_input_file();

     seti.initialize_variables();

     for(n=16; n<=21; n++){
	seti.get_data_length(n);

        for( seti.init_start_sample(); seti.is_sample(); seti.next_sample() ){
	   seti.initialize_data_set1();
	   seti.print_fftblock_header();
           seti.initialize_train_set();
	   seti.set_shadow_dimensions();
           seti.find_confid();
           seti.print_one_fft_header();

           seti.initialize_data_set2();
           seti.initialize_test_set();
           seti.find_coeff();
           seti.delete_train_set();
           seti.delete_test_set();
           seti.interpret_data();
        }

        seti.close_input_file();
	seti.open_input_file();
     }

     seti.put_final_analysis();
     seti.close_output_file();
  }
}