File Coverage

blib/lib/Math/FFT.pm

Criterion	Covered	Total	%
statement	334	370	90.2
branch	77	164	46.9
condition	18	36	50.0
subroutine	31	31	100.0
pod	23	23	100.0
total	483	624	77.4

package Math::FFT;

use strict;
use warnings;

use 5.008;

use vars qw(@ISA);
require DynaLoader;

@ISA = qw(DynaLoader);

# Items to export into callers namespace by default. Note: do not export
# names by default without a very good reason. Use EXPORT_OK instead.
# Do not simply export all your public functions/methods/constants.

our $VERSION = '1.33';

bootstrap Math::FFT $VERSION;

# Preloaded methods go here.

sub new {
    my ($class, $data) = @_;
    die 'Must call constructor with an array reference for the data'
        unless ref($data) eq 'ARRAY';
    $data->[0] ||= 0;    # keep warnings happy
    my $n = @$data;
    my $nip = int(3 + sqrt($n));
    my $nw = int(2 + 5*$n/4);
    my $ip = pack("i$nip", ());
    my $w = pack("d$nw", ());
    return bless {
        type => '',
        mean => '',
        coeff => '',
        n   => $n,
        data => $data,
        ip => \$ip,
        w => \$w,
    }, $class;
}

# clone method to copy the ip and w arrays for data of equal size
sub clone {
    my ($self, $data) = @_;
    die 'Must call clone with an array reference for the data'
        unless ref($data) eq 'ARRAY';
    $data->[0] ||= 0;    # keep warnings happy
    my $n = @$data;
    die "Cannot clone data of unequal sizes" unless $n == $self->{n};
    my $class = ref($self);
    return bless {
        type => '',
        coeff => '',
        mean => '',
        n   => $self->{n},
        data => $data,
        ip => $self->{ip},
        w => $self->{w},
    }, $class;
}

# Complex Discrete Fourier Transform
sub cdft {
    my $self = shift;
    my $n = $self->{n};
    die "data size ($n) must be an integer power of 2" unless _check_n($n);
    my $data = [ @{$self->{data}} ];
    _cdft($n, 1, $data, $self->{ip}, $self->{w});
    $self->{type} = 'cdft';
    $self->{coeff} = $data;
    return $data;
}

# Inverse Complex Discrete Fourier Transform
sub invcdft {
    my $self = shift;
    my $data;
    my $n = $self->{n};
    if (my $arg = shift) {
        if (ref($arg) ne 'ARRAY')
        {
            die 'Must pass an array reference to invcdft';
        }
        if ($n != @$arg)
        {
            die "Size of data set must be $n";
        }
        $data = [ @$arg ];
    }
    else {
        die 'Must invert data created with cdft'
            unless $self->{type} eq 'cdft';
        $data =  [ @{$self->{coeff}} ];
    }
    _cdft($n, -1, $data, $self->{ip}, $self->{w});
    $_ *= 2.0/$n for (@$data);
    return $data;
}

# Real Discrete Fourier Transform
sub rdft {
    my $self = shift;
    my $n = $self->{n};
    if (!_check_n($n))
    {
        die "data size ($n) must be an integer power of 2";
    }
    my $data = [ @{$self->{data}} ];
    _rdft($n, 1, $data, $self->{ip}, $self->{w});
    $self->{type} = 'rdft';
    $self->{coeff} = $data;
    return $data;
}

# Inverse Real Discrete Fourier Transform
sub invrdft {
    my $self = shift;
    my $data;
    my $n = $self->{n};
    if (my $arg = shift) {
        die 'Must pass an array reference to invrdft'
            unless ref($arg) eq 'ARRAY';
        die "Size of data set must be $n"
            unless $n == @$arg;
        $data = [ @$arg ];
    }
    else {
        die 'Must invert data created with rdft'
            unless $self->{type} eq 'rdft';
        $data =  [ @{$self->{coeff}} ];
    }
    _rdft($n, -1, $data, $self->{ip}, $self->{w});
    $_ *= 2.0/$n for (@$data);
    return $data;
}

# Discrete Cosine Transform
sub ddct {
    my $self = shift;
    my $n = $self->{n};
    die "data size ($n) must be an integer power of 2" unless _check_n($n);
    my $data = [ @{$self->{data}} ];
    _ddct($n, -1, $data, $self->{ip}, $self->{w});
    $self->{type} = 'ddct';
    $self->{coeff} = $data;
    return $data;
}

# Inverse Discrete Cosine Transform
sub invddct {
    my $self = shift;
    my $data;
    my $n = $self->{n};
    if (my $arg = shift) {
        die 'Must pass an array reference to invddct'
            unless ref($arg) eq 'ARRAY';
        die "Size of data set must be $n"
            unless $n == @$arg;
        $data = [ @$arg ];
    }
    else {
        die 'Must invert data created with ddct'
            unless $self->{type} eq 'ddct';
        $data =  [ @{$self->{coeff}} ];
    }
    $data->[0] *= 0.5;
    _ddct($n, 1, $data, $self->{ip}, $self->{w});
    $_ *= 2.0/$n for (@$data);
    return $data;
}

# Discrete Sine Transform
sub ddst {
    my $self = shift;
    my $n = $self->{n};
    die "data size ($n) must be an integer power of 2" unless _check_n($n);
    my $data = [ @{$self->{data}} ];
    _ddst($n, -1, $data, $self->{ip}, $self->{w});
    $self->{type} = 'ddst';
    $self->{coeff} = $data;
    return $data;
}

# Inverse Discrete Sine Transform
sub invddst {
    my $self = shift;
    my $data;
    my $n = $self->{n};
    if (my $arg = shift) {
        die 'Must pass an array reference to invddst'
            unless ref($arg) eq 'ARRAY';
        die "Size of data set must be $n"
            unless $n == @$arg;
        $data = [ @$arg ];
    }
    else {
        die 'Must invert data created with ddst'
            unless $self->{type} eq 'ddst';
        $data =  [ @{$self->{coeff}} ];
    }
    $data->[0] *= 0.5;
    _ddst($n, 1, $data, $self->{ip}, $self->{w});
    $_ *= 2.0/$n for (@$data);
    return $data;
}

# Cosine Transform of RDFT (Real Symmetric DFT)
sub dfct {
    my $self = shift;
    my $np1 = $self->{n};
    my $n = $np1 - 1;
    die "data size ($n) must be an integer power of 2" unless _check_n($n);
    my $nt = int(2 + $n/2);
    my $t = [];
    my $data = [ @{$self->{data}} ];
    pdfct($nt, $n, $data, $t, $self->{ip}, $self->{w});
    $self->{type} = 'dfct';
    $self->{coeff} = $data;
    return $data;
}

# Inverse Cosine Transform of RDFT (Real Symmetric DFT)
sub invdfct {
    my $self = shift;
    my $data;
    my $np1 = $self->{n};
    my $n = $np1 - 1;
    if (my $arg = shift) {
        die 'Must pass an array reference to invdfct'
            unless ref($arg) eq 'ARRAY';
        die "Size of data set must be $n"
            unless $np1 == @$data;
        $data = [ @$arg ];
    }
    else {
        die 'Must invert data created with dfct'
            unless $self->{type} eq 'dfct';
        $data =  [ @{$self->{coeff}} ];
    }
    my $nt = int(2 + $n/2);
    my $t = [];
    $data->[0] *= 0.5;
    $data->[$n] *= 0.5;
    pdfct($nt, $n, $data, $t, $self->{ip}, $self->{w});
    $data->[0] *= 0.5;
    $data->[$n] *= 0.5;
    $_ *= 2.0/$n for (@$data);
    return $data;
}

# Sine Transform of RDFT (Real Anti-symmetric DFT)
sub dfst {
    my $self = shift;
    my $n = $self->{n};
    die "data size ($n) must be an integer power of 2" unless _check_n($n);
    my $data = [ @{$self->{data}} ];
    my $nt = int(2 + $n/2);
    my $t = [];
    pdfst($nt, $n, $data, $t, $self->{ip}, $self->{w});
    $self->{type} = 'dfst';
    $self->{coeff} = $data;
    return $data;
}

# Inverse Sine Transform of RDFT (Real Anti-symmetric DFT)
sub invdfst {
    my $self = shift;
    my $n = $self->{n};
    my $data;
    if (my $arg = shift) {
        die 'Must pass an array reference to invdfst'
            unless ref($arg) eq 'ARRAY';
        die "Size of data set must be $n"
            unless $n == @$arg;
        $data = [ @$arg ];
    }
    else {
        die 'Must invert data created with dfst'
            unless $self->{type} eq 'dfst';
        $data =  [ @{$self->{coeff}} ];
    }
    my $nt = int(2 + $n/2);
    my $t = [];
    pdfst($nt, $n, $data, $t, $self->{ip}, $self->{w});
    $_ *= 2.0/$n for (@$data);
    return $data;
}

# check if $n is a power of 2
sub _check_n {
    my ($n) = @_;

    return scalar($n == int($n) and $n > 0 and (! ($n & ($n-1))));
}

sub correl {
    my ($self, $other) = @_;
    my $n = $self->{n};
    my $d1 = $self->{type} ?
    ($self->{type} eq 'rdft' ? [ @{$self->{coeff}} ] :
        die 'correl must involve a real function' ) :
    $self->rdft &&  [ @{$self->{coeff}} ];
    my $d2 = [];
    if (ref($other) eq 'Math::FFT') {
        $d2 = $other->{type} ?
        ($other->{type} eq 'rdft' ? [ @{$other->{coeff}}] :
            die 'correl must involve a real function' ) :
        $other->rdft && [ @{$other->{coeff}}];
    }
    elsif (ref($other) eq 'ARRAY') {
        $d2 = [ @$other ];
        _rdft($n, 1, $d2, $self->{ip}, $self->{w});
    }
    else {
        die 'Must call correl with either a Math::FFT object or an array ref';
    }
    my $corr = [];
    _correl($n, $corr, $d1, $d2, $self->{ip}, $self->{w});
    return $corr;
}

sub convlv {
    my ($self, $r) = @_;
    die 'Must call convlv with an array reference for the response data'
        unless ref($r) eq 'ARRAY';
    my $respn = [ @$r ];
    my $m = @$respn;
    die 'size of response data must be an odd integer' unless $m % 2 == 1;
    my $n = $self->{n};
    my $d1 = $self->{type} ?
    ($self->{type} eq 'rdft' ? [ @{$self->{coeff}} ] :
        die 'correl must involve a real function' ) :
    $self->rdft &&  [ @{$self->{coeff}} ];
    for (my $i=1; $i<=($m-1)/2; $i++) {
        $respn->[$n-$i] = $respn->[$m-$i];
    }
    for (my $i=($m+3)/2; $i<=$n-($m-1)/2; $i++) {
        $respn->[$i-1] = 0.0;
    }
    my $convlv = [];
    _convlv($n, $convlv, $d1, $respn, $self->{ip}, $self->{w});
    return $convlv;
}

sub deconvlv {
    my ($self, $r) = @_;
    die 'Must call deconvlv with an array reference for the response data'
        unless ref($r) eq 'ARRAY';
    my $respn = [ @$r ];
    my $m = @$respn;
    die 'size of response data must be an odd integer' unless $m % 2 == 1;
    my $n = $self->{n};
    my $d1 = $self->{type} ?
    ($self->{type} eq 'rdft' ? [ @{$self->{coeff}} ] :
        die 'correl must involve a real function' ) :
    $self->rdft &&  [ @{$self->{coeff}} ];
    for (my $i=1; $i<=($m-1)/2; $i++) {
        $respn->[$n-$i] = $respn->[$m-$i];
    }
    for (my $i=($m+3)/2; $i<=$n-($m-1)/2; $i++) {
        $respn->[$i-1] = 0.0;
    }
    my $convlv = [];
    if (_deconvlv($n, $convlv, $d1, $respn, $self->{ip}, $self->{w}) != 0) {
        die "Singularity encountered for response in deconvlv";
    }
    return $convlv;
}

{

    my $PI2 = 2 * 4.0 * atan2(1,1);
sub spctrm {
    my ($self, %args) = @_;
    my %accept = map {$_ => 1} qw(window segments number overlap);
    for (keys %args) {
        die "`$_' is not a valid argument to spctrm" if not $accept{$_};
    }
    my $win_fun = $args{window};
    if ($win_fun and ref($win_fun) ne 'CODE') {
        my %accept = map {$_ => 1} qw(hamm hann welch bartlett);
        die "`$win_fun' is not a known window function in spctrm"
        if not $accept{$win_fun};
    }
    die 'Please specify a value for "segments" in spctrm()'
    if ($args{number} and ! $args{segments});
    my $n = $self->{n};
    my $d;
    my $n2 = 0;
    my $spctrm = [];
    my $win_sub = do {
        my $h = sub {
            my ($j, $n) = @_;
            return (1 - cos($PI2*$j/$n))/2;
        };

        +{
            'hamm' => $h,
            'hann' => $h,
            'welch' => sub {
                my ($j, $n) = @_;
                return 1 - 4*($j-$n/2)*($j-$n/2)/$n/$n;
            },
            'bartlett' => sub {
                my ($j, $n) = @_;
                return 1 - abs(2*($j-$n/2)/$n);
            },
        }
    };
    if (not $args{segments} or ($args{segments} == 1 and not $args{number})) {
        die "data size ($n) must be an integer power of 2" unless _check_n($n);
        if ($win_fun) {
            $d = [ @{$self->{data}}];
            $win_fun = $win_sub->{$win_fun} if ref($win_fun) ne 'CODE';
            for (my $j=0; $j<$n; $j++) {
                my $w = $win_fun->($j, $n);
                $d->[$j] *= $w;
                $n2 += $w * $w;
            }
            $n2 *= $n;
            _spctrm($n, $spctrm, $d, $self->{ip}, $self->{w}, $n2, 1);
        }
        else {
            $d = $self->{type} ?
            ($self->{type} eq 'rdft' ? $self->{coeff} :
                die 'correl must involve a real function' ) :
            $self->rdft && $self->{coeff};
            $n2 = $n*$n;
            _spctrm($n, $spctrm, $d, $self->{ip}, $self->{w}, $n2, 0);
        }
    }
    else {
        $d = [ @{$self->{data}}];
        my ($data, @w);
        my $k = $args{segments};
        my $m = $args{number};
        die 'Please specify a value for "number" in spctrm()'
        if ($k and ! $m);
        die "number ($m) must be an integer power of 2" unless _check_n($m);
        my $m2 = $m+$m;
        my $overlap = $args{overlap};
        my $N = $overlap ? ($k+1)*$m : 2*$k*$m;
        die "Need $N data points (data only has $n)" if $N > $n;
        if ($win_fun) {
            $win_fun = $win_sub->{$win_fun} if ref($win_fun) ne 'CODE';
            for (my $j=0; $j<$m2; $j++) {
                $w[$j] = $win_fun->($j, $m2);
                $n2 += $w[$j]*$w[$j];
            }
        }
        else {
            $n2 = $m2;
        }
        if ($overlap) {
            my @old =  splice(@$d, 0, $m);
            for (0..$k-1) {
                push @{$data->[$_]}, @old;
                my @new = splice(@$d, 0, $m);
                push @{$data->[$_]}, @new;
                @old = @new;
                if ($win_fun) {
                    my $j=0;
                    $data->[$_] = [ map {$w[$j++]*$_} @{$data->[$_]}];
                }
            }
        }
        else {
            for (0..$k-1) {
                push @{$data->[$_]}, splice(@$d, 0, $m2);
                if ($win_fun) {
                    my $j=0;
                    $data->[$_] = [ map {$w[$j++]*$_} @{$data->[$_]}];
                }
            }
        }
        my $tmp = [];
        my $nip = int(3 + sqrt($m2));
        my $nw = int(2 + 5*$m2/4);
        my $ip = pack("i$nip", ());
        my $w = pack("d$nw", ());
        _spctrm_bin($k, $m2, $spctrm, $data, \$ip, \$w, $n2, $tmp);
    }
    return $spctrm;
}
}

sub mean {
    my $self = shift;
    my $sum = 0;
    my ($n, $data);
    my $flag = 0;
    if ($data = shift) {
        die 'Must call with an array reference'
        unless ref($data) eq 'ARRAY';
        $n = @$data;
        $flag = 1;
    }
    else {
        $data = $self->{data};
        $n = $self->{n};
    }
    $sum += $_ for @$data;
    my $mean = $sum / $n;
    $self->{mean} = $mean unless $flag == 1;
    return $mean;
}

sub rms {
    my $self = shift;
    my $sum = 0;
    my ($n, $data);
    if ($data = shift) {
        die 'Must call with an array reference'
        unless ref($data) eq 'ARRAY';
        $n = @$data;
    }
    else {
        $data = $self->{data};
        $n = $self->{n};
    }
    $sum += $_*$_ for @$data;
    return sqrt($sum / $n);
}

sub stdev {
    my $self = shift;
    my ($n, $data, $mean);
    if ($data = shift) {
        die 'Must call with an array reference'
        unless ref($data) eq 'ARRAY';
        $n = @$data;
        $mean = $self->mean($data);
    }
    else {
        $data = $self->{data};
        $n = $self->{n};
        $mean = $self->{mean} || $self->mean;
    }
    die 'Cannot find the standard deviation with n = 1'
    if $n == 1;
    my $sum = 0;
    $sum += ($_ - $mean)*($_ - $mean) for @$data;
    return sqrt($sum / ($n-1));
}

sub range {
    my $self = shift;
    my ($n, $data);
    if ($data = shift) {
        die 'Must call with an array reference'
        unless ref($data) eq 'ARRAY';
        $n = @$data;
    }
    else {
        $data = $self->{data};
        $n = $self->{n};
    }
    my $min = $data->[0];
    my $max = $data->[0];
    for (@$data) {
        $min = $_ if $_ < $min;
        $max = $_ if $_ > $max;
    }
    return ($min, $max);
}

sub median {
    my $self = shift;
    my ($n, $data);
    if ($data = shift) {
        die 'Must call with an array reference'
            unless ref($data) eq 'ARRAY';
        $n = @$data;
    }
    else {
        $data = $self->{data};
        $n = $self->{n};
    }
    my @sorted = sort {$a <=> $b} @$data;
    return
    (
        ($n & 0x1)
        ? $sorted[($n-1)/2]
        : ($sorted[$n/2] + $sorted[$n/2-1])/2
    );
}

# Autoload methods go after =cut, and are processed by the autosplit program.


1;

__END__