identiglyph.pl

#!/usr/bin/env perl

use bytes;
use warnings;
use strict;
use POSIX qw(ceil);
use Math::FFT;
use constant { FFTLEN => 2048, SAMPRATE => 44100 };
use constant BINHZ => SAMPRATE / FFTLEN;
use constant BUFZERO => ((0.0) x FFTLEN);

my @samp;
my @fftbuf;
my @symbols;
my $matrix = gen_matrix();
my ($prevsymstart, $numsym) = (-1, 0);
my %symstats;

@ARGV == 1 or print "usage: $0 filename.wav\n" and exit 1;

open IN, "-|:raw", "sox $ARGV[0] -t raw -" or die("error invoking sox: $!");
push @samp, unpack("s*", $_) while sysread(IN, $_, 1024);
close IN;
print "Read ${\ scalar(@samp)} samples\n";
rss_used();
find_symbols();
rss_used();
print "Processed $numsym symbols\n";
print "Symbol statistics:\n";
print "  $_: $symstats{$_} times\n" for (sort { length($a) <=> length($b) || $a cmp $b } keys %symstats);
print "\n";

open OUT, ">", "identiglyph.out" or die("error opening identiglyph.out: $!");
print OUT join('', @symbols);
close OUT;
exit 0;

sub find_symbols
{
  use constant { THRESH => ceil(0.10 * 32768), MINRUNLEN => 40 };
  my ($inrun, $runlen) = (0, 0);

  for (my $i = 0; $i < @samp; $i++) {
    if (abs($samp[$i]) < THRESH) {
      if (!$inrun) {
        $inrun = $runlen = 1;   # begin a run of zeros
      } else {
        $runlen++;              # continue a run
      }
    } else {
      if ($inrun) {             # end of run -- process if long enough
        do_symbol($i, $runlen) if($runlen >= MINRUNLEN);
        $inrun = 0;
      }
    }
  }
  # file probably ended with a run of zeros, so get the last symbol
  do_symbol($#samp+1, $runlen) if($inrun && $runlen >= MINRUNLEN);
}

sub do_symbol
{
  # we're called when at index $i in @samp and we've just seen the end of a $runlen run of zeros
  my ($i, $runlen) = @_;

  # assume file starts with a run of zeros
  $prevsymstart = $i, return if ($prevsymstart == -1);

  my ($symstart, $symend) = ($prevsymstart, $i - $runlen - 1);
  my $len = $symend - $symstart + 1;

  # open(OUT, "|-:raw", "sox -t raw -r 44100 -s -w -c 1 - symbol_wavs/symbol_" . sprintf("%08d", $symstart) . ".wav") or die("error invoking sox: $!");
  # print OUT pack("s*", @samp[$symstart .. $symend]);
  # close(OUT);

  my $bufstart = ceil((FFTLEN - $len)/2);
  $numsym++;
  $prevsymstart = $i;
  print "Warning: symbol start=$symstart end=$symend len=$len probably incorrect, skipping.\n", return if ($len > 750 || $len < 700);

  # FFT buffer is wider than the sample, so zero it out and then copy the sample into the center,
  # converting from signed integers to floating point along the way
  @fftbuf = BUFZERO;
  # $fftbuf[$bufstart + $_] = $samp[$symstart + $_] / 32768.0 for (0..$len-1);  # faster?
  @fftbuf[$bufstart..$bufstart+$len-1] = map { $_ / 32768.0 } @samp[$symstart..$symend];
  my $fft = Math::FFT->new(\@fftbuf);
  my $spectrum = $fft->spctrm(window => "hamm");
  # print "symbol $numsym:\n";
  # printf "           DC: %8f\n", $$spectrum[0];
  # my $h = BINHZ / 2;
  # printf("  %8.2f Hz: %8f\n", $h, $$spectrum[$_]), $h+=BINHZ for (1..$#$spectrum);
  print "Symbol $numsym top bins: ";
  # printf "%.2f Hz, power=%8f  ", (($_-1) * BINHZ + BINHZ/2), $$spectrum[$_] for (sort { $$spectrum[$b] <=> $$spectrum[$a] } (0..$#$spectrum))[0..1];
  # my @topbins = sort { $a <=> $b } (sort { $$spectrum[$b] <=> $$spectrum[$a] } (0..$#$spectrum))[0..2];

  my @bins = ((sort { $$spectrum[$b] <=> $$spectrum[$a] } (28..42))[0], (sort { $$spectrum[$b] <=> $$spectrum[$a] } (50..65))[0]);
  printf "%d/%.2f Hz  ", $_, (($_-1) * BINHZ + BINHZ/2) for @bins;
  $symstats{$_}++ for @bins;
  $symstats{join("x", @bins)}++;
  push @symbols, $matrix->{$bins[0]}{$bins[1]};
  print "\n";
}

sub gen_matrix
{
  my %matrix;
  my $i = 1;

  # these are bin numbers for N=2048, they need scaling if N changes
  # corresponding bin frequencies are (635, 700, 764, 851) and (1109, 1217, 1346) in Hz
  # actual measured frequencies are   (636, 709, 785, 862) and (1109, 1228, 1357) in Hz
  for my $r (30, 33, 36) {
    for my $c (52, 57, 63) {
      $matrix{$r}{$c} = $i++;
    }
  }
  $matrix{40}{57} = 0;

  return \%matrix;
}

sub rss_used
{
  open(STAT, "<", "/proc/$$/stat") and print STDERR sprintf("Resident memory: %.2f MB\n", (((split(/ /, <STAT>))[23]) * 4096 / 1024 / 1024)) and close(STAT);
}