#!/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 take FFT spectrum from instead of standard input
-r calculate RMS
-a calculate A-weighting
-h show column heading (optional, with -r or -a)
-l show as row heading (optional, with -r or -a)
-e extrapolate A curve for f>20kHz (optional, with -a)
-b use as equivalent noise bandwidth of FFT window
(mandatory with -r or -a)
-p rebin spectrum into N bins
-o write rebinned spectrum to instead of standard output
-c 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(){
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;
}
}