http://www.cmlab.csie.ntu.edu.tw/DSPCourse/reference/Sliding%20DFT.pdf

goertzel.pl

#!/usr/bin/env perl
use 5.010;
use strict;
use warnings;

use Carp qw(croak);
use Fcntl qw(:DEFAULT);
use List::MoreUtils qw(pairwise);
use Math::Complex;

my $k = 4;
my $N = 16;

my @input = (read_data(@ARGV))[0 .. 999];

# my @output = dft(@input);

my @output = sliding_dft(@input);

# my @y_real = sliding_goertzel(map { Re($_) } @input);
# my @y_imag = sliding_goertzel(map { Im($_) } @input);
# my @output = pairwise { sqrt($a ** 2 + $b ** 2) } @y_real, @y_imag;

# my @output = sliding_goertzel(@input);

say for @output;

sub read_data {
    my ($filename) = @_;
    my ($buf, @y);
    sysopen(my $fh, $filename, O_RDONLY) || croak "Can't open $filename: $!";
    while (sysread($fh, $buf, 2)) {
        my ($i, $q) = map { ($_ - 127) * (1 / 128) } unpack('CC', $buf);
        push @y, cplx($i, $q);
    }
    close $fh;
    return @y;
}

sub dft {
    my @x = @_;
    my @y;

    for my $n ($N .. $#x) {
        my @samples = @x[($n - $N + 1) .. $n];
        my @freqs;
        $#freqs = $N - 1;
        for my $i (0 .. $#freqs) {
            $freqs[$i] = 0;
            for my $j (0 .. $#samples) {
                my $a = -2 * pi * $i * $j / $N;
                $freqs[$i] += $samples[$j] * cplx(cos($a), sin($a));
            }
        }
        push @y, sqrt(Re($freqs[$k]) ** 2 + Im($freqs[$k]) ** 2);
    }

    return @y;
}

sub sliding_dft {
    my @x = @_;
    my @y;

    my @samples = @x[0 .. ($N - 1)];
    my @freqs;
    $#freqs = $N - 1;
    for my $i (0 .. $#freqs) {
        $freqs[$i] = 0;
        for my $j (0 .. $#samples) {
            my $a = -2 * pi * $i * $j / $N;
            $freqs[$i] += $samples[$j] * cplx(cos($a), sin($a));
        }
    }
    push @y, sqrt(Re($freqs[$k]) ** 2 + Im($freqs[$k]) ** 2);

    my @coeffs;
    for my $i (0 .. $#freqs) {
        $coeffs[$i] = exp(i * 2 * pi * $i / $N);
    }

    for my $n ($N .. $#x) {
        # for my $i (0 .. ($N - 1)) {
        #     $freqs[$i] = ($freqs[$i] - $x[$n - $N] + $x[$n]) * exp(i * 2 * pi * $i / $N);
        # }
        $freqs[$k] = ($freqs[$k] - $x[$n - $N] + $x[$n]) * $coeffs[$k];
        push @y, sqrt(Re($freqs[$k]) ** 2 + Im($freqs[$k]) ** 2);
    }

    return @y;
}

sub goertzel {
    my @x = @_;
    my @y;

    my $coeff = 2 * cos(2 * pi * $k / $N);

    for my $n ($N .. $#x) {
        my ($q0, $q1, $q2) = (0, 0, 0);
        my @samples = @x[($n - $N + 1) .. $n];
        for my $sample (@samples) {
            $q0 = $coeff * $q1 - $q2 + $sample;
            $q2 = $q1; $q1 = $q0;
        }
        push @y, sqrt($q1 ** 2 + $q2 ** 2 - $q1 * $q2 * $coeff);
    }

    return @y;
}

sub sliding_goertzel {
    my @x = @_;
    my @y;

    my $coeff1 = 2 * cos(2 * pi * $k / $N);
    my $coeff2 = exp(i * 2 * pi * $k / $N);

    my ($q0, $q1, $q2) = (0, 0, 0);
    my $y0;
    for my $n ($N .. $#x) {
        $q0 = $coeff1 * $q1 - $q2 + $x[$n] - $x[$n - $N];
        $y0 = $q0 - $coeff2 * $q1;
        $q2 = $q1; $q1 = $q0;

        push @y, sqrt(Re($y0) ** 2 + Im($y0) ** 2);
    }

    return @y;
}

plot.png

short2.dat