206 lines
4.3 KiB
C
206 lines
4.3 KiB
C
/*
|
|
* fft.c
|
|
*
|
|
* Original code by ulrich@Gaston.westfalen.de (Heinz Ulrich Stille).
|
|
*
|
|
* $Id: fft.c,v 1.5 1999/03/16 09:59:19 marcs Exp $
|
|
*
|
|
*/
|
|
#include <stdio.h>
|
|
#include <string.h>
|
|
#include <math.h>
|
|
#include "pvf.h"
|
|
|
|
static void fft (float *real, float *imag, int n)
|
|
{
|
|
float a;
|
|
float tr;
|
|
float ti;
|
|
float wr;
|
|
float wi;
|
|
int mr = 0;
|
|
int l;
|
|
int istep;
|
|
int m;
|
|
int i;
|
|
int j;
|
|
|
|
for (m = 1; m <= (n - 1); m++)
|
|
{
|
|
l = n / 2;
|
|
|
|
while ((mr + l) > (n - 1))
|
|
l /= 2;
|
|
|
|
mr = (mr % l) + l;
|
|
|
|
if (mr > m)
|
|
{
|
|
tr = real[m];
|
|
real[m] = real[mr];
|
|
real[mr] = tr;
|
|
ti = imag[m];
|
|
imag[m] = imag[mr];
|
|
imag[mr] = ti;
|
|
}
|
|
|
|
}
|
|
|
|
for (l = 1; l < n; l = istep)
|
|
{
|
|
istep = 2 * l;
|
|
|
|
for (m = 1; m <= l; m++)
|
|
{
|
|
a = -M_PI * (m - 1.0) / l;
|
|
wr = cos(a);
|
|
wi = sin(a);
|
|
i = m - 1;
|
|
|
|
do
|
|
{
|
|
j = i + l;
|
|
tr = wr * real[j] - wi * imag[j];
|
|
ti = wr * imag[j] + wi * real[j];
|
|
real[j] = real[i] - tr;
|
|
imag[j] = imag[i] - ti;
|
|
real[i] += tr;
|
|
imag[i] += ti;
|
|
i += istep;
|
|
}
|
|
while (i < n);
|
|
|
|
}
|
|
|
|
}
|
|
|
|
}
|
|
|
|
int pvffft (FILE *fd_in, pvf_header *header_in, int skip, int sample_size,
|
|
double threshold, int vgetty_pid, int display)
|
|
{
|
|
|
|
/*
|
|
* Read pvf data, transform it into amplitude data and
|
|
* compute an index for the distribution of energy over the
|
|
* frequencies.
|
|
*/
|
|
|
|
int i;
|
|
int data;
|
|
float val;
|
|
float sum = 0.0;
|
|
float max = 0.0;
|
|
float *real;
|
|
float *imag;
|
|
|
|
for (i = sample_size; i > 0; i >>= 1)
|
|
{
|
|
|
|
if (((i & 0x01) != 0) && (i != 1))
|
|
{
|
|
fprintf(stderr, "%s: sample size (%d) must be a power of 2\n",
|
|
"libpvf", sample_size);
|
|
return(0);
|
|
}
|
|
|
|
}
|
|
|
|
real = malloc((sample_size + 8) * sizeof(float));
|
|
imag = malloc((sample_size + 8) * sizeof(float));
|
|
|
|
if ((real == NULL) || (imag == NULL))
|
|
{
|
|
fprintf(stderr, "%s: not enough memory\n", "libpvf");
|
|
return(0);
|
|
}
|
|
|
|
for(i = 0; i < sample_size; i++)
|
|
{
|
|
real[i] = 0;
|
|
imag[i] = 0;
|
|
}
|
|
|
|
|
|
/*
|
|
* skip the first few seconds
|
|
*/
|
|
|
|
for (i = 0; i < skip; i++)
|
|
data = header_in->read_pvf_data(fd_in);
|
|
|
|
/*
|
|
* store a few seconds' worth of data
|
|
*/
|
|
|
|
for (i = 0; i < sample_size; i++)
|
|
{
|
|
|
|
real[i] = (float) header_in->read_pvf_data(fd_in);
|
|
if (feof(fd_in))
|
|
{
|
|
|
|
/*
|
|
* assume it isn't data
|
|
*/
|
|
|
|
fprintf(stderr, "%s: not enough samples available\n",
|
|
"libpvf");
|
|
return(0);
|
|
}
|
|
|
|
}
|
|
|
|
/*
|
|
* tell the calling process that it can stop writing
|
|
*/
|
|
|
|
if (vgetty_pid > 0)
|
|
kill(vgetty_pid, SIGPIPE);
|
|
|
|
/*
|
|
* note that it will get another sigpipe once this is
|
|
* really finished - a close won't do since other processes
|
|
* (i.e. the calling shell) will still have the fd opened.
|
|
*/
|
|
|
|
fft(real, imag, sample_size);
|
|
|
|
for (i = 10; i < sample_size; i++)
|
|
{
|
|
val = sqrt(real[i] * real[i] + imag[i] * imag[i]);
|
|
sum += val;
|
|
|
|
if (val > max)
|
|
max = val;
|
|
}
|
|
|
|
sum /= sample_size;
|
|
|
|
if (max > 1e-10)
|
|
sum /= max;
|
|
|
|
if (vgetty_pid > 0)
|
|
{
|
|
|
|
if (sum < threshold)
|
|
kill(vgetty_pid, SIGUSR2);
|
|
|
|
return(1);
|
|
}
|
|
|
|
if (display)
|
|
{
|
|
|
|
for (i = 10; i < (sample_size >> 1); i++)
|
|
printf("%f %f\n", ((double) i / sample_size) *
|
|
header_in->speed, sqrt(real[i] * real[i] +
|
|
imag[i] * imag[i]) / max);
|
|
|
|
return(1);
|
|
};
|
|
|
|
printf("%s: FFT level is %g\n", "libpvf", sum);
|
|
return(1);
|
|
}
|