Mail Archives: djgpp/2003/06/15/18:24:29
--=====================_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==_--
- Raw text -