/*

Randomize phase-spectrum of EBS file
Copyright (C) 1993/94 Ralf Schamburger, Institut fuer Physiologie und
     		      Biokybernetik (IPB), Universitaet Erlangen,
		      Deutschland

Permission to use, copy, modify, distribute, and sell this software
and its documentation for any purpose is hereby granted without fee,
provided that the above copyright notice appear in all copies and
related publications and that both that copyright notice and this
permission notice appear in supporting documentation. The authors make
no representations about the suitability of this software for any
purpose. It is provided "as is" without express or implied warranty.

Send your comments, suggestions or bug reports to
        ftpebs@uni-erlangen.de
 
Or mail to
        Institut fuer Physiologie und Biokybernetik
        Markus Prosch
        Universitaetsstrasse 17
        D-91054 Erlangen
        Germany

$Id: rphase.c,v 1.13 1994/02/28 17:01:41 msprosch Exp $

*/


#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#include <string.h>
#include <limits.h>

#include "ebs.h"
#include "ebshead.h" 

#ifdef MSDOS
#error This program needs more memory than DOS provides
#endif

#define MIN(a,b) (((a)<(b))?(a):(b))

#ifndef SEEK_CUR
#define SEEK_CUR 1
#endif

#ifndef RAND_MAX
#define RAND_MAX INT_MAX
#endif

char copyright[]= 
"Randomize phase-spectrum of EBS file\n"
"Copyright (C) 1993/94 Ralf Schamburger, Institut fuer Physiologie und\n"
"                      Biokybernetik (IPB), Universitaet Erlangen, Deutschland\n\n"
"Permission to use, copy, modify, distribute, and sell this software\n"
"and its documentation for any purpose is hereby granted without fee,\n"
"provided that the above copyright notice appear in all copies and\n"
"related publications and that both that copyright notice and this\n"
"permission notice appear in supporting documentation. The authors make\n"
"no representations about the suitability of this software for any\n"
"purpose. It is provided \"as is\" without express or implied warranty.\n\n"
"Send your comments, suggestions or bug reports to\n"
"        ftpebs@uni-erlangen.de\n\n"
"Or mail to\n"
"        Institut fuer Physiologie und Biokybernetik\n"
"        Markus Prosch\n"
"        Universitaetsstrasse 17\n"
"        D-91054 Erlangen\n"
"        Germany\n\n";

short v = 0; /* not verbose is default */

/* Random numbers between 0 and 2pi */
#define frandpi() ((double) ((2.0*M_PI*rand()) / (RAND_MAX)))


/* four1 and realft is taken from the book
         Numerical Recipes in C             */

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

/* fast fourier transform algorithm */ 
void 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;
        }
}

/* fft of a single real function uses four1 */
void realft(float data[], unsigned long n, int isign)

{
  unsigned long i,i1,i2,i3,i4,np3;
  float c1=0.5,c2,h1r,h1i,h2r,h2i;
  double wr,wi,wpr,wpi,wtemp,theta;
  void four1();
  
  theta=3.141592653589793/(double) (n>>1);
  if (isign == 1) {
    c2 = -0.5;
    four1(data,n>>1,1);
  } else {
    c2=0.5;
    theta = -theta;
  }
  wtemp=sin(0.5*theta);
  wpr = -2.0*wtemp*wtemp;
  wpi=sin(theta);
  wr=1.0+wpr;
  wi=wpi;
  np3=n+3;
  for (i=2;i<=(n>>2);i++) {
    i4=1+(i3=np3-(i2=1+(i1=i+i-1)));
    h1r=c1*(data[i1]+data[i3]);
    h1i=c1*(data[i2]-data[i4]);
    h2r = -c2*(data[i2]+data[i4]);
    h2i=c2*(data[i1]-data[i3]);
    data[i1]=h1r+wr*h2r-wi*h2i;
    data[i2]=h1i+wr*h2i+wi*h2r;
    data[i3]=h1r-wr*h2r+wi*h2i;
    data[i4] = -h1i+wr*h2i+wi*h2r;
    wr=(wtemp=wr)*wpr-wi*wpi+wr;
    wi=wi*wpr+wtemp*wpi+wi;
  }
  if (isign == 1) {
    data[1] = (h1r=data[1])+data[2];
    data[2] = h1r-data[2];
  } else {
    data[1]=c1*((h1r=data[1])+data[2]);
    data[2]=c1*(h1r-data[2]);
    four1(data,n>>1,-1);
  }
}


/* calculates the biggest number 2^y but equal or smaller then x */
unsigned long pot2(unsigned long x)
{
  unsigned long res;
  res=pow(2.0,floor(log10((double) x)/log10(2.0)));
  return res;
}

/* reads attributes in inputfile and writes to outputfile */
void readwriteattr(unsigned long attr_len, FILE *fi, FILE *fo)
{
  #define BUFLEN 1024
  long buffer[BUFLEN], mini;
  int sizer, sizew;

  while (attr_len > 0) {
    mini = MIN(BUFLEN,attr_len);
    sizer = fread(buffer, sizeof(long), mini, fi);
    if (sizer > 0) { 
        sizew = fwrite(buffer, sizeof(long), sizer, fo);
	if (sizew != sizer) {
	  fprintf(stderr, "Error: Problems by reading and writing of header.\n");
	  fprintf(stderr, "Program terminates!\n");
	  fclose(fi);
	  fclose(fo);
	  exit(1);
	}
	attr_len -= sizer;
      }
    else {
      fprintf(stderr, "Error: Can't read from input file.\n");
      fprintf(stderr, "Program terminates!\n");
      fclose(fi);
      fclose(fo);
      exit(1);
    }
  }
}

/* copies header from inputfile to outputfile and adds processing-history */
/* void kopfhead(FILE *fin, FILE *fout, char *processmess, unsigned long *n, unsigned long *cm)
{
  unsigned long length;
  unsigned long channels;
  unsigned long magic1, magic2, id;
  unsigned long tag;
  unsigned long attribute_length, attr_help;
  long a,b,a_old,b_old;
  char string[256], string_old[256];
  int fndhis;

  int i,j,k,l,m;

  fndhis = 0; /* found history = 1 */

  /* read fixed header */
/*  magic1 = fgeti32(fin);
  magic2 = fgeti32(fin);
  id = fgeti32(fin);
  channels = fgeti32(fin);
  length = fgeti32(fin);
  if (feof(fin)) {
    fprintf(stderr, "Input file is not an EBS file, too short.\n");
    exit(1);
  }
  if ((magic1 & 0xffffff00) != (EBS_MAGIC1 & 0xffffff00)) {
    fprintf(stderr, "Input file is not an EBS file.\n");
    exit(1);
  }
  if ((magic1 != EBS_MAGIC1) || (magic2 != EBS_MAGIC2)) {
    fprintf(stderr, "Input file looks like an EBS file that has been\n"
                    "corrupted during file transfer.\n");
    exit(1);
  }

  /*  if (id == CIB_16) {
      if (v == 1)
	fprintf(stderr, " EBS data format CIB-16 (channel order, 16-bit bigendian integer)\n"); }
  
  else {
    fprintf(stderr, "EBS data format is not CIB-16. Please use an converter.\n"
	   "Program terminates.");
    exit(1);
  }


  fputi32(magic1, fout);
  fputi32(magic2, fout);
  fputi32(id, fout);
  fputi32(channels, fout);
  fputi32(pot2(length), fout);  /* shorting to a max power of 2 */
  
  
  /* read and write variable header */
/*  while ((tag = fgeti32(fin)) != 0) {
    attribute_length = fgeti32(fin);
    if (tag == PROCESSING_HISTORY) {
      fndhis = 1;
      if (v == 1)
	fprintf(stderr, " Found process entry and adding my own one.\n");
      fputi32(tag, fout);
      attr_help = attribute_length + s4len(processmess);
      fputi32(attr_help, fout);
      readwriteattr(attribute_length, fin ,fout);
      fputs4(processmess, fout);
    }
    else { 
      fputi32(tag, fout);
      fputi32(attribute_length, fout);
      readwriteattr(attribute_length, fin ,fout);
    }
  }

  if (fndhis == 0) {
    if (v == 1)
      fprintf(stderr, " Didn't find process entry - writing a new one.\n");
    fputi32(PROCESSING_HISTORY, fout);
    attribute_length = s4len(processmess);
    fputi32(attribute_length, fout);
    fputs4(processmess, fout);
  }
  tag = 0;
  fputi32(tag, fout);
  /* header is now copied to fout and process-history is up-to-date */
  /* returning # of channels and samples per channel */

/*  *n = channels;
  *cm = length;
}
*/
/* reads one complete channel of the inputfile 
   but only the first lengthp2 numbers are saved in datai */
void readchannel(FILE *fi, unsigned long lengthp2, unsigned long lengthm, short *datai)
{
  unsigned long size, offset;
  int dum;

  size = fread(datai, sizeof(short), lengthp2, fi);
  if (size != lengthp2) {
    fprintf(stderr, "Error: Can't read channel!\n");
    fprintf(stderr, "Program terminates!\n");
    exit(1);
  }
  offset = (lengthm - lengthp2)*sizeof(short);
  dum = fseek(fi, offset, SEEK_CUR);  /* read remainder of channel */
  if (dum != 0) {
    fprintf(stderr, "Error: Can't read remainder of channel!\n");
    fprintf(stderr, "Program terminates!\n");
    exit(1);
  }
}

/* writes the samples from datai in the outputfile */
void writechannel(FILE *fo, unsigned long lengthp2, short *datai)
{
  unsigned long size;
  
  size = fwrite(datai, sizeof(short), lengthp2, fo);
  if (size != lengthp2) {
    fprintf(stderr, "Error: Can't write channel!\n");
    fprintf(stderr, "Program terminates!\n");
    exit(1);
  }
}

/* changes short to float */
void chasho2flo(short *di, float *da, unsigned long newm)
{
  unsigned long y;

  for (y = 0; y < newm; y++) {
	da[y+1] = (float)(di[y]); /* four1 works on 1.. and not on 0.. */
  }
}

/* changes float to short */
void chaflo2sho(short *di, float *da, unsigned long newm)
{
  unsigned long y;

  for (y = 0; y < newm; y++)
	di[y] = (short)((da[y+1])+0.5); /* round  */
}

/* changes phaseinformation in fft-ed data */
void ranphase(float *data, unsigned long len)
{
  unsigned long i;
  double betr, p;
  int bogen;
  float x, y;

  for (i = 2; i <= len/2-1; i++) {
    betr = sqrt((double)(data[2*i-1]*data[2*i-1] + data[2*i]*data[2*i]));
    p = frandpi();
    x = (float)(cos(p)*betr);
    y = (float)(sin(p)*betr);
    data[2*i-1] = x;
    data[2*i] = y;
  }
  /* because of order of realft */
  betr = (double)(data[2]); /* it's the real of len+1 and has no imag */
  p = frandpi();
  x = (float)(cos(p)*betr);
  data[2] = x;
}    

/* norms the samples points */
void norm(float *da, unsigned long len)
{
  unsigned long i;

  for (i = 1; i <= len; i++)
    da[i] = (float)(2.0*da[i]/(double)(len));  /* 2*.... because of realft */
}

/* if -w option set using data windowing */
void welch(float *da, unsigned long len)
{
  unsigned long i;
  double a, wj;
  for (i = 1; i <= len; i++) {
    a = (double)(((double)(i) - 0.5*len)/(0.5*len));
    wj = 1.0 - a*a;
    da[i] *= wj;
  }
}

/* norms samples points if -w option set */
void normwelch(float *da, unsigned long len, double wss)
{
  unsigned long i;

  for (i = 1; i <= len; i++)
    da[i] = 2.0*0.76*da[i]/wss; /* 2*... because of realft */
}

/* calculates energie of one channel */
double energie(short *di, unsigned long len)
{
  unsigned long i;
  double res;

  res = 0;
  for (i = 0; i < len; i++)
    res = res + di[i]*di[i];
return (double)(res);
}

/* how to use rphase */
void printusage(char *name)
{
  fprintf(stderr, "Randomize phase-spectrum of EBS file, V1.0\n"
  "(C) 1993/94, Ralf Schamburger, IPB, Erlangen\n\n"
  "Usage: %s [-vwe] [-t<max>] <inputfile> <outputfile>\n\n", name);
  fprintf(stderr, "options:\n"
  "\t-c           print out copyright notice\n"
  "\t-h           print out this message\n"
  "\t-v           print extra information\n"
  "\t-w           use Welch windowing before transformation\n"
  "\t-e           print energie for each channel before\n"
  "\t             and after transformation\n"
  "\t-t<max>      truncates at positive and negative amplitude\n"
  "\t             (useful in case of artefacts)\n\n");
}


void truncate(short *di, unsigned long len, short trunpo, short trunne)
{
  unsigned long i;

  for (i = 0; i < len; i++)
    if (di[i] > trunpo)
      di[i] = trunpo;
    else if (di[i] < trunne)
      di[i] = trunne;
}




/* *********************************************************************** */
/*                           ** main program **                            */
/***************************************************************************/
int main(int argc, char *argv[])
{
  FILE *fin, *fout;
  EBS_FILE *efin= NULL, *efout;
  char processmess[400];
  unsigned long n,m, newm, i, y;
  float *data;
  short *datai;
  double wss, wj, a, pv, pn;
  short wopt, eopt, fdi, fdo, topt, trun= 0;
  void *memput;
  ebs_fixedheader fix;

  pv = 0;    /* Energie before transformation */
  pn = 0;    /* Energie after transformation */
  wopt = 0;  /* set => -w option set (welch window) */
  eopt = 0;  /* set => -p option set (print energie) */
  fdi = -1;  /* input file not set */
  fdo = -1;  /* output file not set */
  topt = 0;  /* Truncate artefact */
  fin = stdin;

  /* the following tests if the commandline is correct */
  if ((argc == 1) || (argc > 7)) {
    printusage(argv[0]);
    exit(1);
  }

  for (y = 1; y < argc; y++) {
    if (argv[y][0] == '/' && argv[y][1] == '?' && argv[y][2] == 0) {
      printusage(argv[0]);
      exit(1);
    }
    else if (argv[y][0] == '-') 
      switch (argv[y][1]) {
      case 'v':
      case 'V':
	v = 1;
	break;
      case 't':
      case 'T':
	topt=1;
	trun = atoi(argv[y]+2);
	break;
      case 'w':
      case 'W':
        wopt = 1;
        break;
      case 'e':
      case 'E':
        eopt = 1;
        break;
      case '\0':
        if (fdi == 0) {
          printusage(argv[0]);
          exit(1);
        }
        else
	  fdi = 0;  /* input file = stdin */
        break;
      case 'c':
      case 'C':
	fprintf (stderr, copyright);
	exit (1);
      default:
        printusage(argv[0]);
        exit(1);
      }
    else {
      if (fdi == -1)
	fdi = y;
      else
	if (fdo == -1)
	  fdo = y;
	else {
	  printusage(argv[0]);
	  exit(1);
	}
    }
  }
  if (fdi == -1) {
    printusage(argv[0]);
    exit(1);
  }
  if (fdi == 0) {
    fprintf(stderr,"Reading from stdin in new version not yet implemented.\n");
    fprintf(stderr,"Use old version.\n");
    exit(1);
    /* fin = stdin; */
  }
  else
    if (fdi > 0) {
      efin = ebs_open(argv[fdi], 'r');
      if(efin == NULL) {
	fprintf(stderr,"Can't open input file: %s\n",argv[fdi]);
	exit(1);
      }
      if (efin->fixedheader.encoding_id != EBS_CIB_16) {
	fprintf(stderr, "EBS data format is not CIB-16. Please use an");
	fprintf(stderr, " converter.\nProgram terminates.\n");
	ebs_close(efin);
	exit(1);
      }
    }
  if (fdo == -1) {
    fprintf(stderr,"Writting to stdout in new version not yet implemented.\n");
    fprintf(stderr,"Use old version.\n");
    exit(1);
    /* fout = stdout; */
  }
  else {
    efout = ebs_open(argv[fdo], 'w');
    if(efout == NULL) {
      fprintf(stderr,"Can't open output file: %s\n", argv[fdo]);
      ebs_close(efin);
      exit(1);
    }
  }
  /* commandline is correct */
  

  /* built string for processing-history */
  sprintf(processmess, "RPHASE removed phase information:\n ");
  for (y = 0; y < argc; y++) {
    strcat(processmess," ");
    strcat(processmess,argv[y]);
  }
  if ((wopt == 1) && (v == 1))
    fprintf(stderr, " Using Welch window on input data.\n");
  if ((eopt == 1) && (v == 1))
    fprintf(stderr, " Printing energie for each channel.\n");
  if ((topt == 1) && (v == 1))
    fprintf(stderr, " Truncating samples at: %d\n", trun);

  ebs_copyheader(efout, efin);
  memput = ebs_appattr(efout, EBS_PROCESSING_HISTORY, s4len(processmess));
  puts4(processmess, &memput);
  n = efin->fixedheader.channels;
  m = efin->fixedheader.sampleslow;
  ebs_getfixedheader(efout, &fix);
  fix.sampleslow = pot2(m);
  ebs_putfixedheader(efout, &fix);
  ebs_flush(efout);
  fin = ebs_stream(efin);
  fout = ebs_stream(efout);
  ebs_rewind(efin);
  ebs_rewind(efout);
  /*  kopfhead(fin,fout,processmess, &n, &m); */
  /* reads header and writes it to fout with changed processing history
     n: number of channels    m: samples per channel                    */


  if (v == 1)
    fprintf(stderr, " Channels: %d   Samples per channel: %d\n", n ,m);
  newm = pot2(m); /* newm = 2^x =< m  because of using fft */
  if (v == 1)
    fprintf(stderr, " Shorting samples per channel to: %d\n", newm);

  /* getting memory */
  data = (float *)malloc((newm+1)*sizeof(float));
  datai = (short *)malloc((newm+1)*sizeof(short));
  if ((data == NULL) || (datai == NULL)) {
    fprintf(stderr,"Can't get the MEMORY I need!\n");
    exit(1);
  }
 
  /* calculates correcting factor if -w option is set */
  wss = 0.0;
  if (wopt == 1) { /* -w option set */
    for (y = 1; y <= newm; y++) {
          a = (double)(((double)(y) - 0.5*newm)/(0.5*newm));
	  wj = 1.0 - a*a;
	  wss += wj*wj;
	}
  }

  if (v == 1)
    fprintf(stderr, "\n");


  /* the following works over all channels */
  for(i=1; i<=n; i++) {

    if (v == 1)
      fprintf(stderr, "Reading channel:%3d", i); fflush(stderr);
    readchannel(fin, newm, m, datai); /* read chan. and remainder of the ch.*/

    if (topt == 1) 
      truncate(datai, newm, trun, -trun);


    if ((eopt == 1) && (v == 1)) { /* -e option set */
      pv = energie(datai,newm);
    }

    chasho2flo(datai, data, newm);  /* change short to float */

    if (wopt == 1) { /* -w option set */
      welch(data, newm);
    }

    if (v == 1) {
      fprintf(stderr," - transforming");
      fflush(stderr);
    }
    realft(data, newm, 1); 

    if (v == 1) {
      fprintf(stderr, ", changing phase information,\n");
      fflush(stderr);
    }
    ranphase(data, newm);

    if (v == 1) {
      fprintf(stderr, "                      retransforming");
      fflush(stderr);
    }
    realft(data, newm, -1);

    if (wopt!=1) { /* -w option not set */
      norm(data, newm);  /* norm data */
    }

    if (wopt==1) { /* -w option set - norming after Welch */
      normwelch(data, newm, wss);
    }

    chaflo2sho(datai, data, newm);  /* change float to short */

    if (v == 1)
      fprintf(stderr, " and writing back.\n");
    writechannel(fout, newm, datai);

    if ((eopt == 1) && (v == 1)) {  /* energie opt set */
      pn = energie(datai,newm);
      fprintf(stderr, "        Energie before: %9.3E  and after: %9.3E", pv, pn);
      fprintf(stderr, "  transf.\n");
    }

  }

  free(data);
  free(datai);

  ebs_setendof_data_block(efout);
  ebs_close(efin);
  ebs_close(efout);

  if (v == 1)
    fprintf(stderr, "Program successful.\n");

  return 0;
}