# 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;       # sample rate
my $fc          = 2600;        # whistle pitch
my $bps         = 100;         # data speed
my $risetime    = 0.09;        # rise/fall slide duration
my $riseband    = 0.7 * $fc;   # bandwidth thereof
my $padding     = 0.5;         # silence before and after, sec
my $f_am        = 20;          # AM noise frequency
my $vol_am      = 0.18;        # AM noise modulation index
my $sig_vol     = dB(0);       # overall signal gain
my $noise_vol   = dB(12);      # overall noise gain after bandpass
my $noise_shift = 1.0;         # noise shift from harmonic (1 = none)
my $lp_alfa_fm  = 0.004;       # pre-fm lopass bandwidth
my $lp_alfa_noise = 0.009;     # noise bandpass bandwidth
my $fshift      = 0.077 * $fc; # FSK shift

# relative harmonic powers
my @harmonics       = (dB(0), dB(-17), dB(-33), dB(-28), dB(-31), dB(-44));
my @noise_harmonics = (dB(0), dB(-1),  dB(-2),  dB(-3),  dB(-4),  dB(-6));

$data           = "\xAA\xA7" . pack("C", length($data)) . $data;
my $dur         = length($data) * 8 / $bps;

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

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

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

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 = dB(0) - ($t-$dur-$risetime) / $risetime;
    $harmonics[0] = ($t-$dur-$risetime > $risetime * .33 ? 0 :
      dB(0) - (($t-$dur-$risetime) * .33) / ($risetime * .33));

  } else {

    $f = $fc;
    $harmonics[0] = dB(0);

    my $bit_i = int(($t - $risetime) * $bps);
    my $byte  = substr($data, int($bit_i / 8), 1);
    my $bit   = ((unpack("C", $byte) >> ($bit_i % 8)) & 1);

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

    $vol = dB(0);
  }

  $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 * $noise_harmonics[$n] *
      cos(($n + 1) * ($ph_fm * $noise_shift));
  }

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

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

close $stream;

sub dB { 10 ** ($_[0] / 20); }