Analyzing WAV Files with Perl





Analyzing WAV Files with Perl

WAV files created with a microphone and a sound card can be analyzed with a perl program. A description of this process and supporting code are presented.

I have built several electronic measurement instruments that take advantage of an ordinary PC sound capability. The analog to digital conversion done by a sound card provides an easy way to convert real-world signals into computerized measurements.

These instruments create an audio tone, which is recorded and then analyzed with a perl program. I created two types of instruments, AM and FM. In the AM circuits, a tone is either on or off, and the information extracted by the program is the times when the tone transitions from between the on and off states. The FM circuits create a tone which has a pitch which is proportional to the measured quantity. This pitch is tracked and a waveform of pitch versus time is extracted.

The code uses a module called Freqext, named for frequency extraction. The code is simple, procedural, and somewhat tuned to the specific measurements that it was written to extract. Perhaps I will clean it up someday. If I were to do it again, I would look at using some PDL modules that deal with WAV files directly.

Tracking Frequency by Finding Zero Crossings

This example program measures the time between zero crossings to extract the frequency of the signal. This turned out to not be the best approach for the FM frequency extraction. A moving average is used to smooth the results.

#!/usr/bin/perl

use strict;
use Freqext;
$|++;

my $datadir= "/home/toma/perl/examples/wav/data/index.html";
my $datafile= "engine2.wav";

my $filename= "$datadir/$datafile";
print "Extracting file $filename\n";

Freqext::trackz($filename, 10, "/perl/wav/datadir/eng2freq.txt");

Finding the Frequency of a Fixed Tone in a WAV File

This example program uses the FFT to find the frequency of a tone in a WAV file.

#!/usr/bin/perl

use strict;
use Freqext;

my $datadir= "/home/toma/perl/examples/wav/data/index.html";
my @files= qw ( YFTape.wav b1s1.wav b1s2.wav b1s3.wav
b1s4.wav b1s5.wav b1s6.wav b1s7.wav b2s1.wav b2s2.wav b2s3.wav
b2s4.wav b2s5.wav b2s6.wav b2s7.wav big1.wav big10.wav big2.wav
big3.wav big4.wav big5.wav big6.wav big7.wav big8.wav big9.wav
noload.wav small1.wav small2.wav small3.wav small4.wav small5.wav
small6.wav small7.wav );

open(REPORT, "/perl/wav/gtreport.txt") or die "Can't write report";
for (@files)
{
  my $filename= "$datadir/$_";
  print "Extracting file $filename\n";
  my $f= Freqext::stationary($filename);
  print REPORT "$filename $f\n";
}
close REPORT;

Using the FFT to Track a Moving Frequency

This program is able to track the frequency of a signal with the FFT. It tries various FFT sizes on the data files. This worked well.

#!/usr/bin/perl

use strict;
use Freqext;

my $datadir= "/home/toma/perl/examples/wav/data/index.html";
my @files= qw ( 
  engine1.wav engine2.wav engine1.wav engine2.wav engine1.wav
  engine2.wav engine1.wav engine2.wav
);

my @fftsize= qw (
  2048 2048 1024 1024 512 512 256 256
);

my @outfile= qw (
  eng1_2048 eng2_2048 eng1_1024 eng2_1024 eng1_512 eng2_512
  eng1_256 eng2_256
);

for (my $i=0; $i<8; $i++)
{
  for (@files)
  {
    Freqext::moving($datadir."/index.html".$files[$i], $fftsize[$i], $outfile[$i]);
  }
}

Example Output Data

This is the result of FFT frequency versus time extraction.

This is a result of AM detection (code not shown).

Frequency Extraction Module

#
# Module to extract the frequency of a signal in a wav file
# Tue Jan  1 15:47:02 PST 2002
# Copyright 2002 Tom Anderson, All Rights Reserved
#

package Freqext;

use strict;
use warnings;
use diagnostics;

use Audio::Wav;
use PDL;
use PDL::FFT;

sub trackz
{
  my ($input_file, $mov_avg_pts, $output_file) = @_;

  # First estimate the starting frequency

  my $wav = new Audio::Wav;
  my $data= $wav->read( $input_file );  
  my $maxind= $data->length();
  my $seconds= $data->length_seconds();
  my $timestep = $seconds/$maxind;

  my @channels;
  my $v;
  my $last_v=0;
  my $time=0;
  my @zcross;
  my $crosstime;

  for(my $i=0; $i<$maxind; $i++)
  {
    @channels= $data->read();
    $v= $channels[0];
    if ($v > 0 and $last_v <= 0)
    {
      $crosstime= $time-$last_v/($v-$last_v)*$timestep;
      push @zcross, $crosstime;
      # print $crosstime,"\n";
    }
    $last_v= $v;
    $time += $timestep;
  }

  my $i;
  my $total=0;
  my @freqt;
  for ($i=0; $i<$mov_avg_pts; $i++)
  {
    $total += $zcross[$i]
  }
  for ($i=$mov_avg_pts; $i< $#zcross; $i++)
  {
    push @freqt, $total/$mov_avg_pts;
    $total += $zcross[$i]-$zcross[$i-$mov_avg_pts];
  }

  # Now change from list of times of zero crossings to frequency versus time
  open(OUTFILE, ">".$output_file) or die "Can't write $output_file";
  for ($i=1; $i< $#freqt; $i++)
  {
    my $freq= 1/($freqt[$i]-$freqt[$i-1]);
    print OUTFILE $freqt[$i-1],' ',$freq,"\n";
  }
  close OUTFILE;
}

sub tracking
{
  my ($input_file, $chunksize, $output_file) = @_;

  # First estimate the starting frequency

  my $wav = new Audio::Wav;
  my $data= $wav->read( $input_file );  
  my $maxind= $data->length();
  my $seconds= $data->length_seconds();
  my $fsample= $maxind/$seconds;

  my $pointsleft= $maxind;
  my $chunk=0;
  my $break=0;

  # open(OUTFILE, ">".$output_file) or die "Can't write $output_file";

  my @waveform;
  my @channels;
  for(my $i=0; $i<$chunksize; $i++)
  {
    @channels= $data->read();
    push @waveform, $channels[0];
    $pointsleft--;
    # print $channels[0],"\n";
  }

  my $fbin= $fsample/$chunksize;
  my $freq= sequence($chunksize/2) * $fbin;

  # The signal is analytic (imaginary part=0)
  my $v= new PDL(@waveform);

  # Compute the stats on the first chunk of the signal 
  my $pmean= daverage($v);
  $v= $v - $pmean;

  my ($mean, $rms, $median, $min, $max) = stats($v);
  print "Mean: $mean\nRMS: $rms\nMin: $min\nMax: $max\n";

  my $vi= zeroes($chunksize);
  fft($v,$vi);

  # Since signal was analytic, throw out redundant part of fft
  reshape($v, $chunksize/2);
  reshape($vi, $chunksize/2);
  $v= $v**2+$vi**2;

  # This just finds the bin with the highest signal.  
  # Interpolation could improve the accuracy here.
  my $index= maximum_ind($v);
  my $freqest= $index * $fbin;

  print "Freq: $freqest\n";

  # my $n= new PDL(@waveform);
  # $n= ($n-$pmean)/(sqrt(2)*$mean);

  my $x;
  my $scale= 1/(sqrt(2)*$rms);

  my $timestep = 1 / $fsample;
  my $time     = 0;
  my $pi2      = 8 * atan2(1,1);
  my $omest    = $pi2 * $freqest;
  my $sum      = 0;

  my $result="";
  my $cntr=0;
  my $err;
  my $tracksrc;
  my $phase=0;

  # Read the rest of the waveform and track it
  for(my $i=$chunksize; $i<$maxind; $i++)
  {
    @channels = $data->read();
    $x        = ($channels[0] - $pmean) * $scale;
    $tracksrc = sin($phase);

    # Need an FIR filter here that is probably between 1 and 3 cycles long.
    # The ability to look ahead would be nice.
    # Use a nice raised-cosine window function.
    # Still need to analyze the loop gain.

    # Would also be nice to keep the FFT going to make a slower outside 
    # frequency loop to take care of locking on a wrong harmonic.

    $err      = -atan2($tracksrc,$x);
    $sum     += $err;
    $omest   += $sum * .01;
    # $sum      = $sum * .98;
    $freqest  = $omest / $pi2;
    $result  .= $time.' '.$x.' '.$tracksrc.' '.$freqest.' '.$phase.' '.$sum.' '.$err."\n";
    $phase   += $omest * $timestep;
    $phase   -= $pi2 if $phase > $pi2;
    $time    += $timestep;
    if ($cntr++ > 1000)
    {
      $cntr=0;
      print $result;
      $result= "";
    }
  }

  # print OUTFILE "\n";
  # close OUTFILE;
}
sub moving
{
  my ($input_file, $chunksize, $output_file) = @_;

  my $wav = new Audio::Wav;
  my $data= $wav->read( $input_file );  
  my $maxind= $data->length();
  my $seconds= $data->length_seconds();
  my $fsample= $maxind/$seconds;

  my $pointsleft= $maxind;
  my $chunk=0;
  my $break=0;

  open(OUTFILE, ">".$output_file) or die "Can't write $output_file";

  while ($break == 0)
  {
    if ($pointsleft >= $chunksize)
    {
      my @waveform;
      my @channels;
      for(my $i=0; $i<$chunksize; $i++)
      {
        @channels= $data->read();
        push @waveform, $channels[0];
        $pointsleft--;
      }

      my $fbin= $fsample/$chunksize;
      my $freq= sequence($chunksize/2) * $fbin;

      # The signal is analytic (imaginary part=0)
      my $v= new PDL(@waveform);
      my $vi= zeroes($chunksize);
      fft($v,$vi);

      # Since signal was analytic, throw out redundant part of fft
      reshape($v, $chunksize/2);
      reshape($vi, $chunksize/2);
      $v= $v**2+$vi**2;

      # This just finds the bin with the highest signal.  
      # Interpolation could improve the accuracy here.
      my $index= maximum_ind($v);
      my $freqest= $index * $fbin;
  
      my $chunktime= $chunk * $chunksize / $fsample;
      print OUTFILE "$chunktime $freqest\n";
      $chunk++;
    }
    else
    {
      $break=1;
    }
  }
  close OUTFILE;
}

sub stationary
{
  # Routine to return an estimate of the frequency of a stationary
  # signal in a wav file.
  my ($input_file) = @_;
  my $wav = new Audio::Wav;
  my $data= $wav->read( $input_file );  
  my $maxind= $data->length();
  my $seconds= $data->length_seconds();
  my $fsample= $maxind/$seconds;

  my @waveform;
  my @channels;

  # Use up to a closest power of 2 points
  my $npts=1;
  while ($npts < $maxind) { $npts = $npts * 2 }
  $npts= $npts/2;

  for(my $i=0; $i<$npts; $i++)
  {
    @channels= $data->read();
    push @waveform, $channels[0];
  }

  my $fbin= $fsample/$npts;
  my $freq= sequence($npts/2) * $fbin;

  # The signal is analytic (imaginary part=0)
  my $v= new PDL(@waveform);
  my $vi= zeroes($npts);
  fft($v,$vi);

  # Since signal was analytic, throw out redundant part of fft
  reshape($v, $npts/2);
  reshape($vi, $npts/2);
  $v= $v**2+$vi**2;

  # This just finds the bin with the highest signal.  
  # Interpolation could improve the accuracy here.
  my $index= maximum_ind($v);
  my $freqest= $index * $fbin;
  
  return $freqest;
}

1;

06/01/2002

By toma