# Oona Räisänen's pea whistle encoder
# http://www.windytan.com/2015/10/pea-whistle-steganography.html

use warnings;
use strict;

my $data        = $ARGV[0] // "OHAI!";

my $fs          = 44100;
my $fc          = 2600;
my $bps         = 100;
my $risetime    = 0.09;
my $padding     = 0.5;
my $f_am        = 20;
my $vol_am      = 0.18;
my $sig_vol     = 1;
my $noise_vol   = 4;
my $noise_shift = 1.0;
my $lp_alfa_fm  = 0.004;
my $lp_alfa_noise = 0.009;
my $fshift      = 0.0769 * $fc;
my $riseband    = 0.7 * $fc;
$data           = "\xAA\xA7" . pack("C", length($data)) . $data;
my $dur         = length($data) * 8 / $bps;

my @harmonics = (1, 0.141, 0.022, 0.0398, 0.0282, 0.0063);
my @harmonics_noise = (1, 0.9, 0.8, 0.7, 0.6, 0.5);

open my $stream, '|-', "sox -t .s16 -c 1 -r $fs - whistle.wav"
  or die $!;

print $stream pack("s", 0) x ($padding * $fs);

my $t = 0;
my $prev_noise = 0;
my $prev_f = $fc - $riseband;
my $f, my $f_filt, my $vol;
my $ph_am, my $ph_fm;

for (my $t = 0; $t < $risetime*2 + $dur; $t += 1/$fs) {

  if ($t < $risetime) {

    $f = $fc - $riseband + ($t/$risetime * $riseband);
    $vol = $t/$risetime;
    $harmonics[0] = ($t < $risetime * .66 ? 0 :
      ($t-$risetime * .66) / ($risetime * .66));

  } elsif ($t > $risetime + $dur) {

    $f = $fc - ($t-$dur-$risetime) / $risetime * $riseband;
    $vol = 1 - ($t-$dur-$risetime) / $risetime;
    $harmonics[0] = ($t-$dur-$risetime > $risetime * .33 ? 0 :
      1 - (($t-$dur-$risetime) * .33) / ($risetime * .33));

  } else {

    $f = $fc;
    $harmonics[0] = 1;

    my $bitind = int(($t - $risetime) * $bps);
    my $byte   = substr($data, int($bitind / 8), 1);
    my $bit    = ((ord($byte) >> ($bitind % 8)) & 1);

    $f = ($bit - .5) * $fshift + $fc;

    $vol = 1;
  }

  $f_filt = $lp_alfa_fm * $f + (1 - $lp_alfa_fm) * $prev_f;
  $prev_f = $f_filt;

  my $noise = rand() - .5;
  my $noise_filt = $lp_alfa_noise * $noise + (1 - $lp_alfa_noise) *
    $prev_noise;
  $prev_noise = $noise_filt;

  $ph_fm += 2 * 3.14159 * $f_filt / $fs;
  $ph_am += 2 * 3.14159 * $f_am / $fs;

  my $mix = 0;
  for my $n (0..$#harmonics) {
    my $sig = $harmonics[$n] * cos(($n + 1) * $ph_fm);
    $mix += $sig_vol * $sig * (1 + $vol_am * cos($ph_am));
    $mix += $noise_vol * $noise_filt * $harmonics_noise[$n] *
      cos(($n + 1) * ($ph_fm * $noise_shift));
  }

  print $stream pack("s", 0.2 * $vol * $mix * 0x7FFF);
}

print $stream pack("s", 0) x ($padding * $fs);

close $stream;