#!/usr/bin/perl
#
# spectrum.pl
#
# Copyright (C) 2007 Alberto Longhi
#
# This script expects, as input, a text file in 3 columns: frequency,
# FFT level left channel, FFT level right channel.  One such file can
# be created for example with CoolEdit Pro after a "Frequency
# analysis" by doing "Copy to clipboard" and pasting into a text file.
# It is expected that there are 2^n data lines at equidistant
# frequencies starting with zero and ending with the Nyquist
# frequency, that is half of the sampling frequency.  All lines not
# beginning with a digit are ignored.  It is also assumed, as is
# conventional in programs such as CoolEdit, that the FFT level is
# expressed in dB in such a way that the spectrum of a full-scale sine
# wave peaks at 0 dB.
#
# The script does the following:
#
# 1. Calculates the total RMS power of the track expressed in dB
#    relative to the total RMS of a full-scale sine wave.
#
# 2. Calculates the A-weighted level dBA (also relative to the total
#    RMS of a full-scale sine wave).
#
# 3. Rearranges the spectrum so to have data points uniform in log f:
#    this is convenient when plotting the spectrum in logscale, so the
#    graphic is not burdened with overdense points in the high end of
#    the spectrum; when the new binning is sparser the data is
#    averaged and the plot is smoother, where the new binning is
#    denser, a stepping shows up.
#
# Launch without arguments to see a help message.

# -------------------------------------------------------------------
# Subroutines
# -------------------------------------------------------------------

# Help messages
sub usage(){   
    print STDERR << "EOF";
Usage:

spectrum.pl -f filename options
spectrum.pl options < filename

 -f <file>   take FFT spectrum from <file> instead of standard input
 -r          calculate RMS
 -a          calculate A-weighting
 -h          show column heading (optional, with -r or -a)
 -l <label>  show <label> as row heading (optional, with -r or -a)
 -e          extrapolate A curve for f>20kHz (optional, with -a)
 -b <B>      use <B> as equivalent noise bandwidth of FFT window 
             (mandatory with -r or -a)
 -p <N>      rebin spectrum into N bins
 -o <file>   write rebinned spectrum to <file> instead of standard output
 -c <n>      what to show in the output (totals and spectrum)
             1: show only power average of L and R channel
             2: show left and right channel (default)
             3: show left, right and average
 -s          show sum, not average (with -c1 or -c3)
             use -sc 1 if the input has only one channel
EOF
    exit;
}

sub usage_B(){   
    print STDERR << "EOF";
You must specify an equivalent noise bandwidth (ENB).
The ENB depends on what window was used to calculate the FFT.
See the following table for common examples:

WINDOW                      ENB
------------------------    ------
no window or rectangular    1
Bartlett or triangular      1.33
Bartlett-Hann               1.46
Blackman                    1.73
Blackman-Harris             2.01
Blackman-Nuttall            1.98
Flat top                    3.77
Gauss                       (varies)
Hamming                     1.37
Hann                        1.50
Kaiser                      (varies)
Nuttall                     2.02
Sine                        1.24
EOF
}

# A-weighting function
sub w_a(){
    my $f=shift;
    my $f2=$f**2;
    my $a=12200**2;
    my $b=20.6**2;
    my $c=107.7**2;
    my $d=737.9**2;
    my $Ra=$a*$f2**2/(($f2+$b)*($f2+$a)*sqrt($f2+$c)*sqrt($f2+$d));
    return 2.0+20/log(10)*log($Ra);
}

# 10*Log() function
sub log10(){
    my $x=shift;
    return 10/log(10)*log($x);
}

# -------------------------------------------------------------------
# Get options
# -------------------------------------------------------------------

use Getopt::Std;
getopts("f:rahl:eb:p:o:c:s");

# Sanity checks
if(not defined $opt_r and not defined $opt_a and not defined $opt_p){
    print STDERR "Please specify at least one of -r, -a, -p.\n";
    &usage();
}
if(defined $opt_f){
    if(! -f $opt_f){
	print STDERR "File $opt_f not found.\n";
	&usage();
    } else {
	$DATA="< $opt_f";
    }
} else {
    $DATA=STDIN;
}
if(defined $opt_p){
    $N1=int($opt_p);
    if($N1<1){
	print STDERR "$opt_p is not a valid number of bins.\n";
	&usage();
    }
    if(defined $opt_o){
	$OUT="> $opt_o";
    } else {
	$OUT=STDOUT;
    }
    open OUT, $OUT;
}
if(defined $opt_r or defined $opt_a){
    if(!defined $opt_b){
	&usage_B();
	&usage();
    } elsif($opt_b<=0) {
	print STDERR "$opt_b is not a valid ENB.\n";
	&usage_B();
	&usage();
    } else {
	$B=$opt_b;
    }
}
if(defined $opt_c){
    if($opt_c==1 or $opt_c==2 or $opt_c==3){
	$col=($opt_c?$opt_c:2);
    } else {
	print STDERR "$opt_c is not a valid number of output columns.\n";
	&usage();
    }
}
$fact=($opt_s?1:0.5);

# -------------------------------------------------------------------
# Get input data
# -------------------------------------------------------------------

open FH, $DATA;
while(<FH>){
    next unless /^\d/;
    @a=split;
    push @f, $a[0];
    push @datL, $a[1];
    push @datR, $a[2];
}
close FH;

# number of bins
$N=$#f+1;

# Nyquist frequency
$F=$f[$#f];

# -------------------------------------------------------------------
# Total RMS and dBA
# -------------------------------------------------------------------

if($opt_r or $opt_a){

    # initialization
    $PL=$PR=$PLa=$PRa=0;

    # loop
    foreach $j (1..$N-1){
	$PL+=10**($datL[$j]/10);
	$PR+=10**($datR[$j]/10);
	if($f<=20000 or $opt_e){
	    $w=&w_a($f[$j]);
	    $PLa+=10**(($datL[$j]+$w)/10);
	    $PRa+=10**(($datR[$j]+$w)/10);
	}
    }

    # convert to dB
    $dbL=&log10($PL/$B);
    $dbR=&log10($PR/$B);
    $dbLR=&log10(($PL+$PR)*$fact/$B);
    $dbaL=&log10($PLa/$B);
    $dbaR=&log10($PRa/$B);
    $dbaLR=&log10(($PLa+$PRa)*$fact/$B);

    # show results
    $str="%9s";
    $flt="%9.2f";
    if($col==1 and $opt_r and not $opt_a){
	printf STDERR "$str\n", 
	"RMS L+R" if $opt_h;
	printf STDERR "$flt    %s\n",
	$dbLR, $opt_l;
    } elsif($col==1 and not $opt_r and $opt_a){
	printf STDERR "$str\n", 
	"dBA L+R" if $opt_h;
	printf STDERR "$flt    %s\n",
	$dbaLR, $opt_l;
    } elsif($col==1 and $opt_r and $opt_a){
	printf STDERR "$str $str\n", 
	"RMS L+R", "dBA L+R" if $opt_h;
	printf STDERR "$flt $flt    %s\n",
	$dbLR, $dbaLR, $opt_l;
    } elsif($col==2 and $opt_r and not $opt_a){
	printf STDERR "$str $str\n", 
	"RMS L", "RMS R" if $opt_h;
	printf STDERR "$flt $flt    %s\n",
	$dbL, $dbR, $opt_l;
    } elsif($col==2 and not $opt_r and $opt_a){
	printf STDERR "$str $str\n", 
	"dBA L", "dBA R" if $opt_h;
	printf STDERR "$flt $flt    %s\n",
	$dbaL, $dbaR, $opt_l;
    } elsif($col==2 and $opt_r and $opt_a){
	printf STDERR "$str $str $str $str\n", 
	"RMS L", "RMS R", "dBA L", "dBA R" if $opt_h;
	printf STDERR "$flt $flt $flt $flt    %s\n",
	$dbL, $dbR, $dbaL, $dbaR, $opt_l;
    } elsif($col==3 and $opt_r and not $opt_a){
	printf STDERR "$str $str $str\n", 
	"RMS L", "RMS R", "RMS L+R" if $opt_h;
	printf STDERR "$flt $flt $flt    %s\n",
	$dbL, $dbR, $dbLR, $opt_l;
    } elsif($col==3 and not $opt_r and $opt_a){
	printf STDERR "$str $str $str\n", 
	"dBA L", "dBA R", "dBA L+R" if $opt_h;
	printf STDERR "$flt $flt $flt    %s\n",
	$dbaL, $dbaR, $dbaLR, $opt_l;
    } elsif($col==3 and $opt_r and $opt_a){
	printf STDERR "$str $str $str | $str $str $str\n", 
	"RMS L", "RMS R", "RMS L+R", "dBA L", "dBA R", "dBA L+R" if $opt_h;
	printf STDERR "$flt $flt $flt | $flt $flt $flt    %s\n",
	$dbL, $dbR, $dbLR, $dbaL, $dbaR, $dbaLR, $opt_l;
    }
}

# -------------------------------------------------------------------
# Rebin spectrum
# -------------------------------------------------------------------

if($opt_p){

    # Define all frequency bin limits
    $Df=$F/$N;
    for($j=0;$j<=$N;$j++){
	$flim[$j]=($j+0.5)*$Df;
    }
    for($i=0;$i<=$N1;$i++){
	$flim1[$i]=($Df/2)*(1+2*$F/$Df)**($i/$N1);
    }

    # Loop through data
    $jlast=1;
    $PLlast=$PRlast=0;
    for $i (1..$N1){

	# initializations for this pass
	$f1=sqrt($flim1[$i-1]*$flim1[$i]);
	$n=($flim1[$i]-$flim1[$i-1])/$Df;
	$j=$jlast;

	# process data
	if($flim[$j]>$flim1[$i]){
	    $PL=$n * 10**($datL[$j]/10);
	    $PR=$n * 10**($datR[$j]/10);
	} else {
	    $PL=$PLlast;
	    $PR=$PRlast;
	    $j++;
	    while($flim[$j]<$flim1[$i]){
		$PL+=10**($datL[$j]/10);
		$PR+=10**($datR[$j]/10);
		$j++;
	    }
	    if($j<$N){
		$fraz_in = ($flim1[$i]-$flim[$j-1])/$Df;
		$PL += $fraz_in * 10**($datL[$j]/10);
		$PR += $fraz_in * 10**($datR[$j]/10);
	    }
	}

	# print one line of data
	printf OUT "%-.2f\t", $f1;
	if($col==1){
	    printf OUT "%-.5f\n",
	    &log10(($PL+$PR)*$fact/$n);
	} elsif($col==2) {
	    printf OUT "%-.5f\t%-.5f\n",
	    &log10($PL/$n),
	    &log10($PR/$n);
	} elsif($col==3){
	    printf OUT "%-.5f\t%-.5f\t%-.5f\n",
	    &log10($PL/$n),
	    &log10($PR/$n),
	    &log10(($PL+$PR)*$fact/$n);
	}

	# prepare for next pass
	$fraz_out = ($flim[$j]-$flim1[$i])/$Df;
	$PLlast = $fraz_out * 10**($datL[$j]/10);
	$PRlast = $fraz_out * 10**($datR[$j]/10);
	$jlast = $j;
    }
}