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.36';

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__