gamma-ticks.pl

#! /usr/local/bin/perl -w

# Generate a ticking sound where the intervals between ticks are
# distributed following a Gamma distribution.

# Produces raw 48k 16-bit signed single channel audio.  So pipe this to:
# sox -t raw -r 48k -c 1 -e signed -b 16 -L - -t wav ticks.wav

# Command-line parameters are:
# -r <number>: the average tick rate in ticks per second (defaults to 2)
# -s <number>: the INVERSE shape parameter of the Gamma distribution
# -d <number>: the duration in seconds of audio produced
# -0: whether to include a tick at time t=0

use strict;
use warnings;
use Math::Trig qw(:pi);
use POSIX qw(floor tgamma);
use Getopt::Std;

my %opts;
getopts("r:s:d:0", \%opts);

my $samplrate = 48000;
my $inverseshape = ($opts{s} // 1) + 0;
my $rate = ($opts{r} // 2) + 0;
my $scaleparm = (1.0/$rate)*$inverseshape;
my $duration = ($opts{d} // 30) + 0;
my $zerotick = $opts{0};
my $tickfrq = 220;

sub johnk_generate_gamma_small {
    # Returns a gamma-distributed random variable with shape $alpha
    # and scale 1 (expected value $alpha and variance $alpha), for
    # 0<$alpha<1.
    my $alpha = shift;
    die unless $alpha>0 && $alpha<1;
    my $beta = 1-$alpha;
    my ($v1, $v2);
    do {
	my $u1 = rand;
	my $u2 = rand;
	$v1 = $u1 ** (1/$alpha);
	$v2 = $u2 ** (1/$beta);
    } while ( ! ($v1+$v2 <= 1 ) );
    my $z = $v1/($v1+$v2);
    my $u = rand;
    my $e = -log($u);
    return $e * $z;
}

sub generate_gamma {
    # Returns a gamma-distributed random variable with shape $alpha
    # and scale 1 (expected value $alpha and variance $alpha).
    my $alpha = shift;
    die unless $alpha>=0;
    my $n = floor($alpha);
    my $resalpha = $alpha - $n;
    my $x = 0;
    for ( my $i=0 ; $i<$n ; $i++ ) {
	my $u = rand;
	$x += -log($u);
    }
    $x += johnk_generate_gamma_small($resalpha) if $resalpha > 0;
    return $x;
}

print STDERR "rate=$rate\n";
print STDERR "scale=$scaleparm\n";

sub distr {
    # Returns the time until next tick (in seconds).
    return (1.0/$rate) if $scaleparm == 0;
    return (generate_gamma(1/$inverseshape))*$scaleparm;
}

sub tickfunc {
    # The tick wavefunction (returns a pair with first value being the
    # function value and second being false if this tick can be
    # considered to have fully died out).
    my $i = shift;
    my $x = $i/$samplrate;
    my $xx = $x * $tickfrq * 2*pi;
    return (0, 0) if $xx>15;
    return ($xx * exp(1-$xx), 1);
}

my @ticktimes;  # Starting times pipeline of currently active ticks (in samples).
my $nexttick;   # Starting time of next pending tick.
my $lasttick;   # Starting time of last tick (time base).
my $tickcount;  # Tick count added to pipeline.
my $remj = 0;   # Tick count removed from pipeline.

$lasttick = 0;
if ( $zerotick ) {
    push @ticktimes, $lasttick;
    $tickcount++;
}
my $maxsize = scalar(@ticktimes);
my $overflowed = 0;
for ( my $i=0 ; $i<$samplrate*$duration ; $i++ ) {
    if ( ! defined($nexttick) ) {
	# Decide when next tick will occur.
	$nexttick = $lasttick + $samplrate*distr;
    }
    if ( $i >= $nexttick ) {
	# Start the next tick.
	$lasttick = $nexttick;
	push @ticktimes, $lasttick;
	$tickcount++;
	$nexttick = undef;
    }
    $maxsize = scalar(@ticktimes) if scalar(@ticktimes) > $maxsize;
    my $v = 0;  # The value to output
    my $remvticks = 0;
    for ( my $j=0 ; $j<scalar(@ticktimes) ; $j++ ) {
	my ($w, $ct) = tickfunc($i-$ticktimes[$j]);
	$w = -$w if ($j+$remj)%2;  # Every other tick is flipped in direction!
	$v += $w;
	if ( !$ct && $j == $remvticks ) {
	    $remvticks++;
	}
    }
    if ( $remvticks ) {
	# Remove died out ticks
	splice @ticktimes, 0, $remvticks;
	$remj += $remvticks;
    }
    # Convert to 16-bit audio:
    my $wrd = floor($v*15000 + 0.5);
    if ( $wrd >= 32767 ) {
	$overflowed = 1;
	$wrd = 32767;
    } elsif ( $wrd <= -32768 ) {
	$overflowed = 1;
	$wrd = -32768;
    }
    print pack("s<",$wrd);
}
print STDERR "total tick count: $tickcount\n";
print STDERR "ticks completely removed: $remj\n";
print STDERR "max ticks in pipeline: $maxsize\n";
print STDERR "WARNING: overflowed\n" if $overflowed;

weibull-ticks.pl

#! /usr/local/bin/perl -w

# Generate a ticking sound where the intervals between ticks are
# distributed following a Weibull distribution.

# Produces raw 48k 16-bit signed single channel audio.  So pipe this to:
# sox -t raw -r 48k -c 1 -e signed -b 16 -L - -t wav ticks.wav

# Command-line parameters are:
# -r <number>: the average tick rate in ticks per second (defaults to 2)
# -s <number>: the INVERSE shape parameter of the Weibull distribution
# -d <number>: the duration in seconds of audio produced
# -0: whether to include a tick at time t=0

use strict;
use warnings;
use Math::Trig qw(:pi);
use POSIX qw(floor tgamma);
use Getopt::Std;

my %opts;
getopts("r:s:d:0", \%opts);

my $samplrate = 48000;
my $inverseshape = ($opts{s} // 1) + 0;
my $rate = ($opts{r} // 2) + 0;
my $scaleparm = (1.0/$rate)/tgamma(1+$inverseshape);
my $duration = ($opts{d} // 30) + 0;
my $zerotick = $opts{0};
my $tickfrq = 220;

print STDERR "rate=$rate\n";
print STDERR "scale=$scaleparm\n";

sub distr {
    # Returns the time until next tick (in seconds).
    my $t = rand;
    my $w = -log($t);
    return ($w**$inverseshape)*$scaleparm;
}

sub tickfunc {
    # The tick wavefunction (returns a pair with first value being the
    # function value and second being false if this tick can be
    # considered to have fully died out).
    my $i = shift;
    my $x = $i/$samplrate;
    my $xx = $x * $tickfrq * 2*pi;
    return (0, 0) if $xx>15;
    return ($xx * exp(1-$xx), 1);
}

my @ticktimes;  # Starting times pipeline of currently active ticks (in samples).
my $nexttick;   # Starting time of next pending tick.
my $lasttick;   # Starting time of last tick (time base).
my $tickcount;  # Tick count added to pipeline.
my $remj = 0;   # Tick count removed from pipeline.

$lasttick = 0;
if ( $zerotick ) {
    push @ticktimes, $lasttick;
    $tickcount++;
}
my $maxsize = scalar(@ticktimes);
my $overflowed = 0;
for ( my $i=0 ; $i<$samplrate*$duration ; $i++ ) {
    if ( ! defined($nexttick) ) {
	# Decide when next tick will occur.
	$nexttick = $lasttick + $samplrate*distr;
    }
    if ( $i >= $nexttick ) {
	# Start the next tick.
	$lasttick = $nexttick;
	push @ticktimes, $lasttick;
	$tickcount++;
	$nexttick = undef;
    }
    $maxsize = scalar(@ticktimes) if scalar(@ticktimes) > $maxsize;
    my $v = 0;  # The value to output
    my $remvticks = 0;
    for ( my $j=0 ; $j<scalar(@ticktimes) ; $j++ ) {
	my ($w, $ct) = tickfunc($i-$ticktimes[$j]);
	$w = -$w if ($j+$remj)%2;  # Every other tick is flipped in direction!
	$v += $w;
	if ( !$ct && $j == $remvticks ) {
	    $remvticks++;
	}
    }
    if ( $remvticks ) {
	# Remove died out ticks
	splice @ticktimes, 0, $remvticks;
	$remj += $remvticks;
    }
    # Convert to 16-bit audio:
    my $wrd = floor($v*15000 + 0.5);
    if ( $wrd >= 32767 ) {
	$overflowed = 1;
	$wrd = 32767;
    } elsif ( $wrd <= -32768 ) {
	$overflowed = 1;
	$wrd = -32768;
    }
    print pack("s<",$wrd);
}
print STDERR "total tick count: $tickcount\n";
print STDERR "ticks completely removed: $remj\n";
print STDERR "max ticks in pipeline: $maxsize\n";
print STDERR "WARNING: overflowed\n" if $overflowed;

See https://twitter.com/gro_tsen/status/1503704920624619527 for context about this.