Message-Id: <3.0.1.16.20030615182018.34b7942a@earthlink.net>
X-Sender: ethros AT earthlink DOT net
X-Mailer: Windows Eudora Light Version 3.0.1 (16)
Date: Sun, 15 Jun 2003 18:20:18 -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="=====================_1055730018==_"

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

Dear mail list -

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

I've indicated where the debugger says it dies.  The debugger says
"segmemntation error", and further reading of the screen is an indication
of 'General Protection Fault".  

My config.sys file contains FILES=50.

I am totally perplexed.

Can you help??

MUCH THANKS IN ADVANCE.

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


#include "fft.h"


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 Inuput 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.
 *  This program uses only real functions, wso it dpes NOT calculate
 *  the hypotenuse of a triangle
 */
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[])
{
    int i,limit,j,jj;
    complx temp[FFT_LEN];
    complx temper6[FFT_LEN];
    complx temper7[FFT_LEN];        
/*>>>>>>the next line is where it seems to die <<<<<<<<<*/
    char string6[39] = "D:\\OBST\\wave1.dat";   
    char string7[39] = "D:\\OBST\\wave2.dat";
    /*char string3[32] ="D:\\OBST\\fft.out";  */
    FILE *fptr6, *fptr7; 
    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;
}

--=====================_1055730018==_
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);

--=====================_1055730018==_--