Spectrograms with the FFTW FFT-Spectrograph library and GNUplot on Linux
FFTW (fftw.org) is a very nice library to perform FFT (Fast Fourier Transforms).
This page presents a quick hack to perform a FFT over a wav-file.
The C-program reads the data from the wav-file, performs the FFT and writes the raw data and
the FFT values to "data.raw" and "data.dat". The perl script invokes GNUplot - GNUplot uses
"data.dat" to plot the spectrum as PS (postscript) and PNG image file. This is a very convenient
way to perform FFT periodically, i.e. as a spectrum analyzer cron job or other means of batch processing data.
The wav-file used in this example is available here (187kB).
(sound-sample source: MPIfR LOPES Project)
Compile fftux.c with compile.sh. Download the wav-sample and save
it into the directory with fftux.c. Run "fftux" - this will produce data.raw and data.dat. Run
"generate-plots.pl" - this will produce "data.ps" and "data.png".
If all went well, you get this PNG-image with the power spectrum: (the original au-file had a sampling rate
of only 4000/s, hence the limitation of the plot to 2000Hz)
fftux.c:
// FFTW demo with wav file
// (C) 2005 www.captain.at
#include <stdio.h>
#include <stdlib.h>
#include <fftw3.h>
static double *rdata = NULL, *idata = NULL;
static fftw_plan rplan, iplan;
static int last_fft_size = 0;
double A[8192];
int n = 8192;
int i,m,k;
char erg[35];
char erg2[35];
double absval;
int is16bitfile = 1;
int main() {
double cc;
int v, val, val2;
char inc[n];
FILE *infile;
FILE *file;
FILE *file2;
int samplerate = 44100;
int bandwidth = samplerate / 2;
int nhalf = n / 2;
double correction = (double)samplerate / (double)n;
infile = fopen("theevent-zoom.wav", "r");
file = fopen("data.dat", "w");
file2 = fopen("data.raw", "w");
// skip wav file header
fgets(inc,37,infile);
fgets(inc,9,infile);
// read data and fill "A"-array
for (v=0;v<n;v++) {
val = fgetc(infile);
if (is16bitfile) {
val2 = fgetc(infile);
if (val2 > 127) { val2 = val2 - 255; }
A[v] = 256*val2 + val;
} else {
A[v] = val;
}
sprintf(erg2, "%d %f\n", v, A[v]);
fputs(erg2, file2);
}
// prepare fft with fftw
rdata = (double *)fftw_malloc(n * sizeof(double));
idata = (double *)fftw_malloc(n * sizeof(double));
// create the fftw plan
rplan = fftw_plan_r2r_1d(n, rdata, idata, FFTW_R2HC, FFTW_FORWARD);
// we have no imaginary data, so clear idata
memset((void *)idata, 0, n * sizeof(double));
// fill rdata with actual data
for (i = 0; i < n; i++) { rdata[i] = A[i]; }
// make fft
fftw_execute(rplan);
// post-process FFT data: make absolute values, and calculate
// real frequency of each power line in the spectrum
m = 0;
for (i = 0; i < (n-2); i++) {
absval = sqrt(idata[i] * idata[i]);
cc = (double)m * correction;
sprintf(erg, "%f %f\n", cc, absval);
fputs(erg, file);
m++;
}
// housekeeping
fclose(file);
fclose(file2);
fclose(infile);
fftw_destroy_plan(rplan);
fftw_free(rdata);
fftw_free(idata);
return 1;
}
compile.sh: (shell script to compile fftux.c)
#!/bin/sh
gcc -o fftux fftux.c -lfftw3 -lm
generate-plots.pl: (perl script to generate plots)
#!/usr/bin/perl
# Generate postscript and png plot with GNUplot
# (C) 2005 www.captain.at
# set custom font path
$ENV{GDFONTPATH} = "/usr/share/fonts/truetype/ttf-bitstream-vera/";
# GNUPLOT POSTSCRIPT
open (GNUPLOT, "|gnuplot");
print GNUPLOT <<EOPLOT;
set term post color "Courier" 12
set output "data.ps"
set size 1 ,1
set nokey
set data style line
set xlabel "frequency [Hz]" font "Courier,14"
set xrange [0:2000]
set yrange [0:4000000]
set title "FFT" font "Courier,14"
set grid xtics ytics
set xtics 100
plot "data.dat" using 1:2 w lines 1
EOPLOT
close(GNUPLOT);
# GNUPLOT PNG
open (GNUPLOT, "|gnuplot");
print GNUPLOT <<EOPLOT;
set term png small xFFFFFF font "VeraMono" 8
set output "data.png"
set size 1 ,1
set nokey
set data style line
set xlabel "frequency [Hz]" font "VeraMono,10"
set xrange [0:2000]
set yrange [0:4000000]
set title "FFT" font "VeraMono,10"
set grid xtics ytics
set xtics 100
plot "data.dat" using 1:2 w lines 1
EOPLOT
close(GNUPLOT);
Last-Modified: Sat, 04 Feb 2006 16:03:16 GMT
