fft_dump.c

#include <rfftw.h>
#include <math.h>
#include <stdio.h>
/* the out format is: */
/* r0, r1, r2, ..., rn/2, i(n+1)/2-1, ..., i2, i1 */

void fft(int n, fftw_real in[n], fftw_real power_spectrum[n/2+1])
{
    fftw_real out[n];
    rfftw_plan p;
    int k = 0;

    p = rfftw_create_plan(n, FFTW_REAL_TO_COMPLEX, FFTW_ESTIMATE);
    rfftw_one(p, in, out);

    for (k = 1; k < (n+1)/2; ++k)  /* (k < n/2 rounded up) */
    {
/*         power_spectrum[k] = (out[k]*out[k] + out[n-k]*out[n-k]); */
        power_spectrum[k] = 10*log(out[k]*out[k] + out[n-k]*out[n-k])/log(10.0);
    }
    if (n % 2 == 0) /* n is even */
    {
        power_spectrum[n/2] = out[n/2]*out[n/2];  /* nyquist freq. */
    }
    fftw_destroy_plan(p);  
}

#define BYTES_PER_SAMPLE (2)
#define SAMPLE_RATE (44100)
#define N (4096*2)
#define PI (3.1415926)

#define FREQ_MIN (150)
#define FREQ_MAX (2000)

void int16_to_real(int num, char *in, double *out)
{
    int i = 0;
    for (i = 0; i < num; i++)
    {
      out[i] = in[i*2] + (in[i*2+1] << 8);
    }
}

void write_sine_int16(FILE *out, int n, double freq, int sample_rate, double power)
{
    int i = 0;
    double x = 0;
    int i_x = 0;
    double step = ((double)(freq))/sample_rate;
    for (i = 0; i < n; i++)
    {
      x = (sin(i*step*2*PI))*(50*power); //((1<<14)-1);
        i_x = (int)round(x);
        fputc((char)(i_x >> 8), out);
        fputc((char)(i_x & 0xFF), out);
    }
}

int main(const int argc, const char * const argv [])
{
    FILE *pcm_file = NULL;
    char data[N*BYTES_PER_SAMPLE] = {0};
    fftw_real real_data[N] = {0}, power_spectrum[N/2+1] = {0};
    
    int k = 0;
    int pos = 0;
    double max_freq = 0, max_mag = 0, freq = 0;//, prev_freq = 0;
    
    if (argc != 2) {
        fprintf(stderr, "Usage: %s <filename>\n", argv[0]);
        return -1;
    }
    
    pcm_file = fopen(argv[1], "rb");
    if (NULL == pcm_file)
    {
        return -1;
    }

    pos = 0;
    //    prev_freq = 0;
    while (!feof(pcm_file)) {
        pos += fread(data, sizeof(char), N*BYTES_PER_SAMPLE, pcm_file);
        int16_to_real(N, data, real_data);
        fft(N, real_data, power_spectrum);

        max_mag = 0;
        max_freq = 10;
        for (k=0;k<N/2; k++) {
            freq = ((double)k)/N*SAMPLE_RATE;
/*             printf("%f ", freq); */
/*             printf("%f\n", power_spectrum[k]); */
            if ((freq < FREQ_MAX) && (freq > FREQ_MIN) && (power_spectrum[k] > 20)) {
                if ((max_mag == 0) || (power_spectrum[k] > max_mag)) {
                    max_freq = freq;
                    max_mag = power_spectrum[k];
                }
            }
        }
/*         printf("%d %f\n", pos, max_freq); */
        write_sine_int16(stdout, N, max_freq, SAMPLE_RATE, max_mag/100);
        //        prev_freq = max_freq;
/*         break; */
    }
                  
    return 0;
}