#!/usr/bin/perl
# Compute an FFT plot of a wav file.

use Audio::Wav;
use Math::FFT;
use Data::Dumper;

(my $myname = $0) =~ s|(.*/)*||;        # strip path component from name
my $Usage = "Usage: $myname [samples_to_read] \n";
# Check one or more arguments:
die $Usage if ($#ARGV < 0);

$file = $ARGV[0];
print "Reading $file\n";

my $wav = new Audio::Wav;
my $read = $wav -> read($file);
my $base = $file;
$base =~ s/\.wav$//;

my $max = 0;
$max = $ARGV[1] if defined($ARGV[1]);

my $info = $read->details();
my $freq = $$info{sample_rate}; # Sample frequency (Hz)
my $time = $$info{length};      # Length of wav file in seconds

my @data = ();
my $datum;
my $n = 0;
while(@datum = ($read->read())) {
  push(@data, $datum[0]);
  $n++;
  last if $max && ($n > $max);
}
$n--;

print "Read $n samples\n";

my $total = 1;
$total *=2 while $total*2 <= $n/2;

$total = 65536 if $total > 65536;

print "Computing for $total samples\n";
#@dat = @data[0..$total-1];

my $fft = new Math::FFT(\@data);
#my %options = (window => "bartlett");
my %options = ();
my $segments = int(($#data + 1)/$total) + 1;
print "Processing $segments segments.\n";
$options{segments} = $segments;
$options{number} = $total/2;
$options{overlap} = 1;
my $spectrum = $fft->spctrm(%options);


log_scale($spectrum);

open(OUT, ">$base.dat") or die;
# Each unit in the spectrum is $freq/($#$spectrum+1)/2 Hz wide:
my $unit = $freq/($#$spectrum+1)/2;
foreach $n (0..$#$spectrum) {
  print OUT $unit*$n, " $$spectrum[$n]\n";
}
close(OUT);

open(PLOT, "|gnuplot") or die;
print PLOT <<EOF;
#set term post color "Courier" 12
#set output "$base.ps"
set term gif size 1024, 768
set output "$base.gif"
set size 1 ,1
set nokey
set style data lines
set xlabel "frequency [Hz]" font "Courier,14"
#set xrange [0:15]
#set yrange [0:100]
#set yrange [0:80]
#set xrange [490:500]

set xrange [10:20000]
set logscale x

#set xrange [2200:2400]
#set xrange [1100:1200]

set title "FFT" font "Courier,14"
set grid xtics ytics
set style fill solid 0.25 border
#set style fill solid 1
#plot "$base.dat" using 1:2 w lines 1
plot "$base.dat" using 1:2 with boxes
EOF
close(PLOT);
#system qq[xli "$base.gif" &];


# Convert spectrum data to relative dB values ###

sub log_scale($) {
  my ($spectrum) = @_;
  my $cutoff = 1;
  $_ = abs($_) for @$spectrum;
  #print "Spectrum elements = ", $#$spectrum+1, "\n";
  #print "cutoff = $cutoff\n";
  # Normalise results to dB values
  my $log10 = log(10);
  for (@$spectrum) {
    $_ = abs($_) / $cutoff;
    if ($_ < 1) {
      $_ = 0;
    } else {
      # Convert to dB values relative to the cutoff value:
      $_ = 20 * log($_)/$log10;
    }
  }
}