Message-Id: <3.0.1.16.20030616153401.34b73854@earthlink.net>
X-Sender: ethros AT earthlink DOT net
X-Mailer: Windows Eudora Light Version 3.0.1 (16)
Date: Mon, 16 Jun 2003 15:34:01 -0400
To: djgpp AT delorie DOT com
From: Ethan Rosenberg <ethros AT earthlink DOT net>
Subject: Program
Mime-Version: 1.0
Content-Type: multipart/mixed; boundary="=====================_1055806441==_"

--=====================_1055806441==_
Content-Type: text/plain; charset="us-ascii"

Dear mail list -

My apologies for sending this again, but 1]I have not received a reply
(though it might be in the mail-digest that will arrive tonight) and 2]I
have made some changes, and now an other error appears.

I am trying to run the attached program using djgpp and RHIDE.

My config.sys file contains FILES=50.

I changed stdio.h so that the max number of files is 50.  I made some other
changes, which I think I have reversed.  Just in case I did not, I am
appending the current version of stdio.h.

In reply to my original post, I received the following from a friend:

**********
To your string (char[]) constants within the function (stuff2), try 
adding the following modifiers before the "char" keyword: static and 
const. It may work with just the "static" keyword. If not, I would make 
the char constants global. You can then import them through the function 
interface, or if you don't like using globals directly, set them up to 
the locals but define them as char * locally. If you have any questions 
give me a call in the morning (after 9am).

Explanation: The reason compiler is not being friendly to your code is 
that a string of 39 characters is being allocated on your stack 
every time you call stuff2. It could be that the stack ran out of memory 
when it hit your code. If you make variables in a function static, they 
will only be allocated once but on the heap. This should not be used for 
recursive functions where the values may change in the process of the 
recursion. However, constant string variables are a good case of static 
when used within a function (outside of functions, static tells the 
compiler that the variable name is restricted to within the file).

***************
  
I followed the above advice, and the program ran u UNTIL the line in
stuff2.c which said "jj = 0;"  At that point, the program exited with an
exit code of 255.

I have now removed the references to string6 and string7, and the
associated write statements.

I had added the write statements, and the two files to debug stuff2.c.
With the program exiting as it does, I can't even debug it manually.

If it is truly a stack/heap problem, how do I increase the size of the
stack/heap?

I am totally perplexed.

Can you help??

MUCH THANKS IN ADVANCE.

Ethan Rosenberg
--=====================_1055806441==_
Content-Type: text/plain; charset="us-ascii"

/* Copyright (C) 1998 DJ Delorie, see COPYING.DJ for details */
/* Copyright (C) 1995 DJ Delorie, see COPYING.DJ for details */

#ifndef __dj_include_stdio_h_
#define __dj_include_stdio_h_

#ifdef __cplusplus
extern "C" {
#endif

#ifndef __dj_ENFORCE_ANSI_FREESTANDING

#include <sys/version.h>
#include <sys/djtypes.h>
  
#define _IOFBF    	00001
#define _IONBF    	00002
#define _IOLBF    	00004

/* Some programs think they know better... */
#undef NULL

#define BUFSIZ		16384
#define EOF		(-1)
#define FILENAME_MAX	260
#define FOPEN_MAX	50
#define L_tmpnam	260
#define NULL		0
#define TMP_MAX		999999

#define SEEK_SET	0
#define SEEK_CUR	1
#define SEEK_END	2

__DJ_va_list
#undef __DJ_va_list
#define __DJ_va_list
__DJ_size_t
#undef __DJ_size_t
#define __DJ_size_t

/* Note that the definitions of these fields are NOT guaranteed!  They
   may change with any release without notice!  The fact that they
   are here at all is to comply with ANSI specifictions. */
   
typedef struct {
  int   _cnt;
  char *_ptr;
  char *_base;
  int   _bufsiz;
  int   _flag;
  int   _file;
  char *_name_to_remove;
  int   _fillsize;
} FILE;

typedef unsigned long		fpos_t;

extern FILE __dj_stdin, __dj_stdout, __dj_stderr;
#define stdin	(&__dj_stdin)
#define stdout	(&__dj_stdout)
#define stderr	(&__dj_stderr)

void	clearerr(FILE *_stream);
int	fclose(FILE *_stream);
int	feof(FILE *_stream);
int	ferror(FILE *_stream);
int	fflush(FILE *_stream);
int	fgetc(FILE *_stream);
int	fgetpos(FILE *_stream, fpos_t *_pos);
char *	fgets(char *_s, int _n, FILE *_stream);
FILE *	fopen(const char *_filename, const char *_mode);
int	fprintf(FILE *_stream, const char *_format, ...);
int	fputc(int _c, FILE *_stream);
int	fputs(const char *_s, FILE *_stream);
size_t	fread(void *_ptr, size_t _size, size_t _nelem, FILE *_stream);
FILE *	freopen(const char *_filename, const char *_mode, FILE *_stream);
int	fscanf(FILE *_stream, const char *_format, ...);
int	fseek(FILE *_stream, long _offset, int _mode);
int	fsetpos(FILE *_stream, const fpos_t *_pos);
long	ftell(FILE *_stream);
size_t	fwrite(const void *_ptr, size_t _size, size_t _nelem, FILE *_stream);
int	getc(FILE *_stream);
int	getchar(void);
char *	gets(char *_s);
void	perror(const char *_s);
int	printf(const char *_format, ...);
int	putc(int _c, FILE *_stream);
int	putchar(int _c);
int	puts(const char *_s);
int	remove(const char *_filename);
int	rename(const char *_old, const char *_new);
void	rewind(FILE *_stream);
int	scanf(const char *_format, ...);
void	setbuf(FILE *_stream, char *_buf);
int	setvbuf(FILE *_stream, char *_buf, int _mode, size_t _size);
int	sprintf(char *_s, const char *_format, ...);
int	sscanf(const char *_s, const char *_format, ...);
FILE *	tmpfile(void);
char *	tmpnam(char *_s);
int	ungetc(int _c, FILE *_stream);
int	vfprintf(FILE *_stream, const char *_format, va_list _ap);
int	vprintf(const char *_format, va_list _ap);
int	vsprintf(char *_s, const char *_format, va_list _ap);

#ifndef __STRICT_ANSI__

#define L_ctermid
#define L_cusrid
/* #define STREAM_MAX	20 - DOS can change this */

int	fileno(FILE *_stream);
FILE *	fdopen(int _fildes, const char *_type);
int	pclose(FILE *_pf);
FILE *	popen(const char *_command, const char *_mode);

#ifndef _POSIX_SOURCE

extern FILE __dj_stdprn, __dj_stdaux;
#define stdprn	(&__dj_stdprn)
#define stdaux	(&__dj_stdaux)

#define P_tmpdir "c:/"

void	_djstat_describe_lossage(FILE *_to_where);
int	_doprnt(const char *_fmt, va_list _args, FILE *_f);
int	_doscan(FILE *_f, const char *_fmt, void **_argp);
int	_doscan_low(FILE *, int (*)(FILE *_get), int (*_unget)(int, FILE *), const char *_fmt, void **_argp);
int	fpurge(FILE *_f);
int	getw(FILE *_f);
int	mkstemp(char *_template);
char *	mktemp(char *_template);
int	putw(int _v, FILE *_f);
void	setbuffer(FILE *_f, void *_buf, int _size);
void	setlinebuf(FILE *_f);
char *	tempnam(const char *_dir, const char *_prefix);
int	_rename(const char *_old, const char *_new);	/* Simple (no directory) */
int	vfscanf(FILE *_stream, const char *_format, va_list _ap);
int	vscanf(const char *_format, va_list _ap);
int	vsscanf(const char *_s, const char *_format, va_list _ap);

#endif /* !_POSIX_SOURCE */
#endif /* !__STRICT_ANSI__ */
#endif /* !__dj_ENFORCE_ANSI_FREESTANDING */

#ifndef __dj_ENFORCE_FUNCTION_CALLS
#endif /* !__dj_ENFORCE_FUNCTION_CALLS */

#ifdef __cplusplus
}
#endif

#endif /* !__dj_include_stdio_h_ */

--=====================_1055806441==_
Content-Type: text/plain; charset="us-ascii"


#include "fft.h"

/*char string6[] = "D:\\OBST\\wavein.dat";   
char string7[] = "D:\\OBST\\waveout.dat";*/
/*static char *string6;
static char *string7; */         

static FLOAT	 	*SineTable, *CosineTable, *WinTable;

# define		SINE(x)		SineTable[(x)]
# define		COSINE(x)       CosineTable[(x)]
# define		WINDOW(x)	WinTable[(x)]


/*
 *  This PERMUTE could be defined as a table lookup (Sampson did this),
 *  but I found that to be significantly slower than just computing it
 *  inline.  Of course, this is w/ GCC, and not MS-DOS Turbo C.
 *  I'll leave it as a macro; you can tweak with it.
 */

#define PERMUTE(x, y)	reverse((x), (y))

/* Number of samples in one "frame" */
#define SAMPLES     8192
#define REAL(x)		wave[(x)].re
#define IMAG(x)		wave[(x)].im
#define ALPHA		0.54	/* ala hanning */
#define FFT_LEN     16384
#define STEP        320

void main(void)
{
FILE *fptr1, *fptr2, *fptr3, *fptr4, *fptr5;
int i,j,k,l,m,ii,jj,kk,num_read,z,len,inum,icycles;
float num,cycles;
char string1[32];/* ="d:\\obst\\delta.dat";*/
char string2[32];/* ="D:\\OBST\\w4.bin";        */
char string3[32] ="D:\\OBST\\fft.out";
char string4[32] ="d:\\obst\\psd.dat";
char string5[32];/* ="d:\\obst\\psdlot.dat";*/
FLOAT  data[SAMPLES],sum1,sum2,avg,index,freqavg, fft_out[FFT_LEN/2],temp;
FLOAT  temp2[FFT_LEN],e_low,e_high,e_tot,delta,min,delta_avg,std_dev;
FLOAT  avg_del[30],minplus;
complx wave[FFT_LEN];
long pos_byte,pos;
struct jha_data
{                                                                 
    short fhr;
    short toco;
    short emg;
};
struct jha_data data_in[SAMPLES];

    printf("Enter Input Data File [including path]:  ");
    gets(string2);
    printf("Enter DELTA Output Data File [including path]:  ");
    gets(string1);
    printf("Enter LOTUS Input Data File [including path]:   ");
    gets(string5);
    
    for (i=0;i<SAMPLES;i++)
        data[i]=0.0;
    for (i=0;i<31;i++)
        avg_del[i]=0.0;
    delta = 0.0;
    num_read = 0;
    sum1 = 0.0;
    ii = 0;            
    kk = 0;               
    minplus = 0.0;
    sum1 = sizeof(struct jha_data);
    
    if( (fptr1 = fopen(string2, "rb")) == NULL)
    {
        printf("\x7");
        perror("Error in opening source file\n");
        fclose(fptr1);
        exit(1);
    }
    else
    {
    printf("The source file WAS opened\n");
    }

/*    if( (fptr2 = fopen(string3, "wb")) == NULL)
    {
        printf("\x7");
        perror("Error in opening output file\n");
        fclose(fptr2);
        exit(1);
    }
    else
    {
        printf("The output file WAS opened\n");
    }*/

    if( (fptr3 = fopen(string4, "wb")) == NULL)
    {
        printf("\x7");
        perror("Error in opening STD DEV  file\n");
        fclose(fptr3);
        exit(1);
    }
    else
    {
        printf("The STD DEV file WAS opened\n");
    }

    if( (fptr4 = fopen(string1, "ab")) == NULL)
    {
        printf("\x7");
        perror("Error in opening DELTA file\n");
        fclose(fptr4);
        exit(1);
    }

    else
    {
        printf("The DELTA file WAS opened\n");
    }

    if( (fptr5 = fopen(string5, "a")) == NULL)
    {
        printf("\x7");
        perror("Error in opening LOTUS INPUT file\n");
        fclose(fptr5);
        exit(1);
    }
    else
    {
        printf("The LOTUS INPUT file WAS opened\n");
    }

    fft_init(13);

    while(!feof(fptr1))
    {
       /* if(kk > 180)printf("%d\n",ii);
        if(ii == 17 && kk>180)
           z=1;*/
/* Use this when importing from DADiSP fileas with ONE data stream        */

        pos_byte = ftell(fptr1);
        pos = pos_byte/sizeof(FLOAT);
        min = (FLOAT)pos/1920.0;
        num_read = fread(data,sizeof(FLOAT),SAMPLES,fptr1);
        if(num_read < SAMPLES)
            break;
        fseek(fptr1,-(SAMPLES-STEP)*sizeof(FLOAT),SEEK_CUR);
        for (i = 0;i < SAMPLES; i++)
        {
            wave[i].re = data[i];/*(float)data_in[i].emg;*/
            wave[i].im = 0.0;
        }                                      
/* Use this when importing from DAQWARE fileas with three data streams        */
        
/*        pos_byte = ftell(fptr1);
        pos = pos_byte/sizeof(struct jha_data);
        min = (FLOAT)pos/1920.0;
        num_read = fread(data_in,sizeof(struct jha_data),SAMPLES,fptr1);
        if(num_read < SAMPLES)
            break;
        fseek(fptr1,-(SAMPLES-STEP)*sizeof(struct jha_data),SEEK_CUR);
        for (i = 0;i < SAMPLES; i++)
        {
            wave[i].re = 0.0;/*(float)data_in[i].emg;*/
/*            wave[i].im = 0.0;
        }*/
        scale(wave, 13);
        window(wave,13);
        stuff2(wave);
        zero_cmplx(wave);
        
        fft(wave,14);
        for (i = 0;i < FFT_LEN; i++)
            fft_out[i] = amp(i, wave, 14);

/*/        fwrite(fft_out,sizeof(float),FFT_LEN/2,fptr2);*/
        for(i=0;i<FFT_LEN/2;i++)
            temp2[i] = fft_out[i]*fft_out[i]*2;
        fwrite(temp2,sizeof(float),FFT_LEN/2,fptr3);      
        e_low = trap_integ(temp2,8,18,0.2,0.45);
/*        printf("\n E_Low = %f\n ",e_low);*/
        e_high = trap_integ(temp2,33,123,0.45,3.0);
/*        printf("E_High = %f\n ",e_high);   */
        e_tot = trap_integ(temp2,8,123,0.2,3.0);
/*        printf("E_Tot = %f\n ",e_tot);       */
        delta = -14.9317*(e_low/e_tot) + 137.681*(e_high/e_tot);
       if(ii <30)
       {
           avg_del[ii] = delta;
           ii++;
           kk++;
       }
       if(ii==1)break;
  /*     if(ii == 30)                                           
         {
              minplus = min + (1.0/6.0);
              printf("DELTA = %f  min = %f  minplus = %f\n ",delta,min,minplus);         
       }       */
       if(ii == 30)
      {
          for(jj = 0;jj < 30; jj++)
              sum1 += avg_del[jj];
          delta_avg = sum1/30.0;
          std_dev = std(avg_del,30);
          if(kk == 180)
            index = 0.0;
/*          fprintf(fptr3,"%.2f  %.4f  %.4f %d  %d  %.4f\n",\
          minplus,delta_avg,std_dev);/*fptr3*/   
          /*printf(/"%.2f  %.4f  %.4f %d  %d  %.4f\n",min,delta_avg,std_dev,ii,kk,delta);*/             
          ii = 0;
      }
        fwrite(&delta,sizeof(FLOAT),1,fptr4);
        fprintf(fptr5,"%.2f  %.4f\n",min,delta);
    }
    fclose(fptr1);
/*    fclose(fptr2);   */
    fclose(fptr3);   
    fclose(fptr4);
    fclose(fptr5);
    free(SineTable);
    free(CosineTable);
    free(WinTable);
    printf("\nSuccessful Completion!!!");
}

/*
 *  Bit reverser for unsigned ints
 *  Reverses 'bits' bits.
 */
const unsigned int reverse (unsigned int val, int bits)
{
  unsigned int retn = 0;

  while (bits--)
    {
      retn <<= 1;
      retn |= (val & 1);
      val >>= 1;
    }
  return (retn);
}

/*
 *  Here is the real work-horse.
 *  It's a generic FFT, so nothing is lost or approximated.
 *  The samples in wave[] should be in order, and they
 *  will be decimated when fft() returns.
 */
void fft (complx wave[], int bits)
{
  register int	loop, loop1, loop2;
  unsigned	i1;		/* going to right shift this */
  int		i2, i3, i4, y;
  FLOAT		a1, a2, b1, b2, z1, z2;
/*  extern float SineTable[], CosineTable[];*/
  i1 = FFT_LEN / 2;
  i2 = 1;

  /* perform the butterfly's */

  for (loop = 0; loop < bits; loop++)
    {
      i3 = 0;
      i4 = i1;

      for (loop1 = 0; loop1 < i2; loop1++)
	{
	  y  = PERMUTE(i3 / (int)i1, bits);
	  z1 = CosineTable[y];
	  z2 = -SineTable[y];

	  for (loop2 = i3; loop2 < i4; loop2++)
	    {
	      a1 = REAL(loop2);
	      a2 = IMAG(loop2);

	      b1 = z1 * REAL(loop2+i1) - z2 * IMAG(loop2+i1);
	      b2 = z2 * REAL(loop2+i1) + z1 * IMAG(loop2+i1);

	      REAL(loop2) = a1 + b1;
	      IMAG(loop2) = a2 + b2;

	      REAL(loop2+i1) = a1 - b1;
	      IMAG(loop2+i1) = a2 - b2;
	    }

	  i3 += (i1 << 1);
	  i4 += (i1 << 1);
	}

      i1 >>= 1;
      i2 <<= 1;
    }
}

/*
 *  Put the samples back in order after the FFT scrambles them.
 *  (Decimation-in-frequecy)
 *  Untested and unused because I haven't done IFFT yet.
 */
void reorder (complx wave[], int bits)
{
  int i, j;
  complx temp;

  for (i = 0; i < SAMPLES; i++)
    {
      j = PERMUTE (i, bits);
      if (i > j)		/* So we only do each pair once */
	{
	  temp = wave[i];
	  wave[i] = wave[j];
	  wave[j] = temp;
	}
    }
}

/*
 *  Scale sampled values.
 *  Do this *before* the fft.
 */
void scale (complx wave[], int bits)
{
  int i;

  for (i = 0; i < SAMPLES; i++)
  {
    wave[i].re /= SAMPLES;
    wave[i].im /= SAMPLES;
  }
}

/*
 *  Calculate the hypotenuse of a right triangle.  The function in math.h
 *  assumes the the input and return values are double.  This uses FLOAT
 */
FLOAT hypot2(FLOAT a, FLOAT b)
{
  FLOAT c,asq,bsq;
  double aa,bb,cc;
  asq = a*a;
  bsq = b*b;
  aa = (double)asq;
  bb = (double)bsq;
  cc = sqrt(aa+bb);
  c = (FLOAT)sqrt((double)(a*a) + (double)(b*b));
  return c;
}


/*
 *  Calculate amplitude of component n in the decimated wave[] array.
 */

FLOAT amp (int n, complx wave[], int bits)
{
  n = PERMUTE (n, bits);
  return (hypot2 (REAL(n), IMAG(n)));

/*  return (REAL(n), IMAG(n)); */
}


/*
 *  Calculate phase of component n in the decimated wave[] array.
 */
FLOAT phase (int n, complx wave[], int bits)
{
  n = PERMUTE (n, bits);
  if (REAL(n) != 0.0)
    return (atan (IMAG(n) / REAL(n)));
  else
    return (0.0);
}


/*
 *  Initializer for FFT routines.  Currently only sets up tables.
 *  - Generates scaled lookup tables for sin() and cos()
 *  - Fills a table for the Hamming window function
 */
void fft_init (int bits)
{
  int i;

  const FLOAT  	TWOPIoN   = (atan(1.0) * 8.0) / (FLOAT)SAMPLES;
  const FLOAT 	TWOPIoNm1 = (atan(1.0) * 8.0) / (FLOAT)(SAMPLES - 1);
  const FLOAT   TWOPISAM  = 2.0*PI/(FLOAT)SAMPLES;

  SineTable   = (FLOAT *)malloc (sizeof(FLOAT) * SAMPLES);
  CosineTable = (FLOAT *)malloc (sizeof(FLOAT) * SAMPLES);
  WinTable    = (FLOAT *)malloc (sizeof(FLOAT) * SAMPLES);

  for (i=0; i < SAMPLES; i++)
    {
      SineTable[i]   = sin((FLOAT) i * TWOPIoN);
      CosineTable[i] = cos((FLOAT) i * TWOPIoN);
      /*
       * Generalized Hanning window function.
       * Set ALPHA to 0.54 for a Hamming window. (Good idea)
       */
      WinTable[i] = 0.5 - 0.5*(cos(TWOPISAM*(FLOAT)i));

    }
}


/*
 *  Apply some windowing function to the samples.
 */
void window (complx wave[], int bits)
{
  int i;

  for (i = 0; i < SAMPLES; i++)
    {
      REAL(i) *= WinTable[i];
      IMAG(i) *= WinTable[i];
    }
}

/*
 *  Undo the window function to restore full magnitude.
 *  Only works if there are NO zeros in WinTable!  (like hanning's)
 */
void unwindow (complx wave[], int bits)
{
  int i;


  for (i = 0; i < SAMPLES; i++)
    {
      REAL(i) /= WinTable[i];
      IMAG(i) /= WinTable[i];
    }
}

/*  Pad the end of the input with zeros */

void stuff(complx wave[])
{
    int i;

    for(i = SAMPLES; i < FFT_LEN; i++)
    {
      REAL(i) = 0.0;
      IMAG(i) = 0.0;
    }
}
/*  Pad the begining and end of the input with zeros */

void stuff2(complx wave[])
{
/*    FILE *fptr6, *fptr7; */
    int i,limit,j,jj;
    complx temp[FFT_LEN];
    complx temper6[FFT_LEN];
    complx temper7[FFT_LEN];        
/*    extern char string6;/* = "D:\\OBST\\wave1.dat";   */
/*    extern char string7;/* = "D:\\OBST\\wave2.dat";*/
    /*char string3[32] ="D:\\OBST\\fft.out";  */              
/*  SineTable   = (FLOAT *)malloc (sizeof(FLOAT) * SAMPLES);*/
                                                              
 /*   string6 = (char *)malloc(sizeof(char) * 34);
    string7 = (char *)malloc(sizeof(char) * 34);
    string6 = "D:\\OBST\\wavein.dat";
    string7 = "D:\\OBST\\waveout.dat";*/
    
    jj=0;    
/*    if( (fptr6 = fopen(string6, "wb")) == NULL)
    {
        printf("\x7");
        perror("Error in opening WAVE IN file\n");
        fclose(fptr6);
        exit(1);
    }
    else
    {
        printf("The WAVE IN file WAS opened\n");
    }

    if( (fptr7 = fopen(string7, "wb")) == NULL)
    {
        printf("\x7");
        perror("Error in opening WAVE OUT  file\n");
        fclose(fptr7);
        exit(1);
    }
    else
    {
        printf("The WAVE OUT file WAS opened\n");
    }*/

    limit = (FFT_LEN - SAMPLES)/2;
  
    
    for(i = 0; i < FFT_LEN; i++)
    {
      temp[i].re = 0.0;
      temp[i].im = 0.0;
    }
    for(i = 0; i < FFT_LEN; i++)
    {
      temper6[i].re = 0.0;
      temper6[i].im = 0.0;
    }
    for(i = 0; i < FFT_LEN; i++)
    {
      temper7[i].re = 0.0;
      temper7[i].im = 0.0;
    }
    for(i = 0; i < FFT_LEN; i++)
    {
      temper6[i].re = wave[i].re;
      temper6[i].im = wave[i].im;
    }
/*    while(!feof(fptr6))
       fwrite(temper6,sizeof(float),1,fptr6);*/
    
    for(i = 0; i < limit; i++)
    {
      temp[i].re = 0.0;
      temp[i].im = 0.0;
    }
    
    for(i = limit, j = 0; i < limit+SAMPLES, j < SAMPLES; i++,j++)
    {
      temp[i].re = wave[j].re;
      temp[i].im = wave[j].im;
    }
    for(i = limit+SAMPLES; i < FFT_LEN; i++)
    {
      temp[i].re = 0.0;
      temp[i].im = 0.0;
    }
    for(i = 0; i < FFT_LEN; i++)
    {
      wave[i].re = temp[i].re;
      wave[i].im = temp[i].im;
    }
    for(i = 0; i < FFT_LEN; i++)
    {
      temper7[i].re = wave[i].re;
      temper7[i].im = wave[i].im;
    }
/*    while(!feof(fptr7))
       fwrite(temper7,sizeof(float),1,fptr7);  
    fclose(fptr6);
    fclose(fptr7);                             */
}

/* Insure that the complex component of the wave is zero */

void zero_cmplx(complx wave[])
{
    int i;
    for (i = 0;i < FFT_LEN; i++)
    {
        wave[i].im = 0.0;
    }
}

/* Integrate the PSD function, using the trapaziodal rule */

FLOAT trap_integ(FLOAT temp2[],int start_pt,int end_pt,FLOAT start_val,FLOAT end_val)
{
    FLOAT sum;
    int i;

    sum = 0.0;
    sum += temp2[start_pt];
    for(i = start_pt+1; i < end_pt; i++)
        sum += 2*temp2[i];
    sum += temp2[end_pt];
    sum = (end_val-start_val)/(2*((float)end_pt-(float)start_pt))*sum;
    return sum;
}

/*
 * Calculate the std deviation of the delta calculation
 */

FLOAT std(FLOAT avg_del[], int nn)
{              

    int i,jj;
    FLOAT sum1,sumx2,std_dev,sum2;

    sum1  = 0.0;
    sumx2 = 0.0;
    sum2 = 0.0;

    for(jj = 0;jj < nn; jj++)
       sum1 += avg_del[jj];

    for(i = 0; i<nn; i++)
        sumx2 += avg_del[i]*avg_del[i];
    std_dev = (FLOAT)sqrt(((double)sumx2 -\
             (double)sum1*(double)sum1/(double)nn)/((double)nn-1));

    return std_dev;
}

--=====================_1055806441==_
Content-Type: text/plain; charset="us-ascii"

/* fft.h: Defs for fft routines
 * [Part of simple-fft-1.5.tar.Z (ver 1.5)]*/

#include <stdio.h>
#include <math.h>
/*#include <libc.h>*/

/*
 * If you want double-precision floats used everywhere,
 * define this to be "double"
 */
#define FLOAT 		float

/*
 * Since every math.h seems to have their own idea of what the
 * elements in struct complex are called, we define our own.
 */
struct _complx
{
  FLOAT re;
  FLOAT im;
};
typedef struct _complx complx;

/*
 * Functions provided by fft.c
 */
void		fft (complx wavetrum[], int bits);
void		ifft (complx wavetrum[], int bits);
void		scale (complx wavetrum[], int bits);
void		window (complx wavetrum[], int bits);
void		unwindow (complx wavetrum[], int bits);
void		fft_init (int bits);
FLOAT		phase (int n, complx wavetrum[], int bits);
FLOAT		amp (int n, complx wavetrum[], int bits);
void        reorder (complx wave[], int bits);
const unsigned int reverse (unsigned int val, int bits);
FLOAT   hypot2(FLOAT a, FLOAT b);
void    stuff(complx wave[]);
void    zero_cmplx(complx wave[]);
void    stuff2(complx wave[]);
FLOAT   trap_integ(FLOAT temp2[],int start_pt,int end_pt,FLOAT start_val,FLOAT end_val);
FLOAT   std(FLOAT avg_del[], int nn);

--=====================_1055806441==_--