File Coverage

blib/lib/Math/WalshTransform.pm

Criterion	Covered	Total	%
statement	82	288	28.4
branch	15	66	22.7
condition			n/a
subroutine	14	18	77.7
pod	12	16	75.0
total	123	388	31.7

# Math::WalshTransform.pm
#########################################################################
#        This Perl module is Copyright (c) 2002, Peter J Billam         #
#               c/o P J B Computing, www.pjb.com.au                     #
#                                                                       #
#     This module is free software; you can redistribute it and/or      #
#            modify it under the same terms as Perl itself.             #
#########################################################################

package Math::WalshTransform;
no strict;
$VERSION = '1.17';
# gives a -w warning, but I'm afraid $VERSION .= ''; would confuse CPAN
require DynaLoader;
require Exporter;
@ISA = qw(Exporter DynaLoader);
@EXPORT = qw(fht fhtinv fwt fwtinv biggest logical_convolution
 logical_autocorrelation power_spectrum walsh2hadamard hadamard2walsh);
@EXPORT_OK = qw(sublist distance size average normalise product arr2txt);
%EXPORT_TAGS = (ALL => [@EXPORT,@EXPORT_OK]);
bootstrap Math::WalshTransform $VERSION;

$PP = 0;

sub fht {
	if (! $PP) { return &xs_fht(@_); }
	my @mr = @_;
	my $k = 1;
	my $n = scalar @mr;
	my $l = $n;
	my $i; my $nl; my $nk; my $j;
	while () {
		$i = $[-1; $l = $l/2;
		for ($nl=1; $nl<= $l; $nl++) {
			for ($nk=1; $nk<= $k; $nk++) {
				$i++; $j = $i+$k;
				$mr[$i] = ($mr[$i] + $mr[$j])/2;
				$mr[$j] =  $mr[$i] - $mr[$j];
			}
			$i = $j;
		}
		$k = 2*$k;
		last if $k >= $n;
	}
	if ($k == $n) {
		return @mr;
	} else {
		warn "Math::WalshTransform::fht \$n = $n but must be power of 2\n";
		return ();
	}
}
sub fhtinv {
	if (! $PP) { return &xs_fhtinv(@_); }
	my @mr = @_;
	my $k = 1;
	my $n = scalar @mr;
	my $l = $n;
	my $i; my $nl; my $nk; my $j;
	while () {
		$i = $[-1; $l = $l/2;
		for ($nl=1; $nl<= $l; $nl++) {
			for ($nk=1; $nk<= $k; $nk++) {
				$i++; $j = $i+$k;
				$mr[$i] = $mr[$i] + $mr[$j];
				$mr[$j] = $mr[$i] - 2*$mr[$j];
			}
			$i = $j;
		}
		$k = 2*$k;
		last if $k >= $n;
	}
	if ($k == $n) {
		return @mr;
	} else {
		warn "Math::WalshTransform::fhtinv \$n = $n but must be power of 2\n";
		return ();
	}
}
sub fwt { # might be easier to Hadamard transform and shuffle the results
	if (! $PP) { return &xs_fwt(@_); }
	my @mr = @_;
	my $n = scalar @mr; my @nr; $#nr=$#mr;
	my $k; my $l; my $i; my $nl; my $nk; my $kp1;

	my $m = 0;  # log2($n)
	my $tmp = 1; while () { last if $tmp >= $n; $tmp<<=1; $m++; }
	my $alternate = $m & 1;

	if ($alternate) {
		for ($k=$[; $k<$n-$[; $k+=2) {
			$kp1 = $k+1;
			$mr[$k]   = ($mr[$k] + $mr[$kp1])/2;
			$mr[$kp1] =  $mr[$k] - $mr[$kp1];
		}
	} else {
		for ($k=$[; $k<$n-$[; $k+=2) {
			$kp1 = $k+1;
			$nr[$k]   = ($mr[$k] + $mr[$kp1])/2;
			$nr[$kp1] =  $nr[$k] - $mr[$kp1];
		}
	}

	$k = 1; my $nh = $n/2;
	while () {
		my $kh = $k; $k = $k+$k; $kp1 = $k+1; last if $kp1>$n;
		$nh = $nh/2; $l = $[; $i = $[; $alternate = !$alternate;
		for ($nl=1; $nl<=$nh; $nl++) {
			for ($nk=1; $nk<=$kh; $nk++) {
				if ($alternate) {
					$mr[$l]   = ($nr[$i]   + $nr[$i+$k])/2;
					$mr[$l+1] =  $mr[$l]   - $nr[$i+$k];
					$mr[$l+2] = ($nr[$i+1] - $nr[$i+$kp1])/2;
					$mr[$l+3] =  $mr[$l+2] + $nr[$i+$kp1];
				} else {
					$nr[$l]   = ($mr[$i]   + $mr[$i+$k])/2;
					$nr[$l+1] =  $nr[$l]   - $mr[$i+$k];
					$nr[$l+2] = ($mr[$i+1] - $mr[$i+$kp1])/2;
					$nr[$l+3] =  $nr[$l+2] + $mr[$i+$kp1];
				}
				$l = $l+4; $i = $i+2;
			}
			$i = $i+$k;
		}
	}
	return @mr;
}
sub fwtinv {
	if (! $PP) { return &xs_fwtinv(@_); }
	my @mr = @_; my $n = scalar @mr; my @nr; $#nr=$#mr;
	my $k; my $l; my $i; my $nl; my $nk; my $kp1;

	my $m = 0;  # log2($n)
	my $tmp = 1; while () { last if $tmp >= $n; $tmp<<=1; $m++; }
	my $alternate = $m & 1;

	if ($alternate) {
		for ($k=$[; $k<$n-$[; $k+=2) {
			$kp1 = $k+1;
			$mr[$k]   =  $mr[$k] + $mr[$kp1];
			$mr[$kp1] =  $mr[$k] - $mr[$kp1] - $mr[$kp1];
		}
	} else {
		for ($k=$[; $k<$n-$[; $k+=2) {
			$kp1 = $k+1;
			$nr[$k]   =  $mr[$k] + $mr[$kp1];
			$nr[$kp1] =  $mr[$k] - $mr[$kp1];
		}
	}

	$k = 1; my $nh = $n/2;
	while () {
		my $kh = $k; $k = $k+$k; $kp1 = $k+1; last if $kp1>$n;
		$nh = $nh/2; $l = $[; $i = $[; $alternate = !$alternate;
		for ($nl=1; $nl<=$nh; $nl++) {
			for ($nk=1; $nk<=$kh; $nk++) {
				if ($alternate) {
					$mr[$l]   =  $nr[$i]   + $nr[$i+$k];
					$mr[$l+1] =  $nr[$i]   - $nr[$i+$k];
					$mr[$l+2] =  $nr[$i+1] - $nr[$i+$kp1];
					$mr[$l+3] =  $nr[$i+1] + $nr[$i+$kp1];
				} else {
					$nr[$l]   =  $mr[$i]   + $mr[$i+$k];
					$nr[$l+1] =  $mr[$i]   - $mr[$i+$k];
					$nr[$l+2] =  $mr[$i+1] - $mr[$i+$kp1];
					$nr[$l+3] =  $mr[$i+1] + $mr[$i+$kp1];
				}
				$l = $l+4; $i = $i+2;
			}
			$i = $i+$k;
		}
	}
	return @mr;
}

sub logical_convolution { my ($xref, $yref) = @_;
	if (ref $xref ne 'ARRAY') { warn
	"Math::WalshTransform::logical_convolution 1st arg must be array ref\n";
		return undef;
	} elsif (ref $yref ne 'ARRAY') { warn
	"Math::WalshTransform::logical_convolution 2nd arg must be array ref\n";
		return undef;
	}
	my @Fx = &fwt(@$xref);  my @Fy = &fwt(@$yref);
	# my @Fz; foreach ($[ .. $#Fx) { $Fz[$_] = $Fx[$_] * $Fy[$_]; } return @Fz;
	return &fwtinv(&product(\@Fx, \@Fy));
}

sub old_logical_convolution { my ($xref, $yref) = @_;
	if (ref $xref ne 'ARRAY') { warn
	"Math::WalshTransform::logical_convolution 1st arg must be array ref\n";
		return undef;
	} elsif (ref $yref ne 'ARRAY') { warn
	"Math::WalshTransform::logical_convolution 2nd arg must be array ref\n";
		return undef;
	}
	local $[ = 0;
	my @x = @$xref; my @y = @$yref;
	my $n = scalar @x;
	my @z; $#z=$#x;
	my $j; my $k; my $sum;
	for ($k=$[; $k<=$#x; $k++) {
		$sum = 0.0;
		for ($j=$[; $j<=$#x; $j++) { $sum += $x[$j^$k] * $y[$j]; }
		$z[$k] = $sum/$n;
	}
	return @z;
}
sub logical_autocorrelation {
	&logical_convolution(\@_,\@_);
}
sub power_spectrum {
	&fwt( &logical_convolution(\@_,\@_) );
}
sub walsh2hadamard {
	my @h; $#h = $#_;
	my @jw2jh = &jw2jh(scalar @_);
	my $i; for ($i=$[; $i<=$#_; $i++) { $h[$jw2jh[$i]] = $_[$i]; }
	return @h;
}
sub hadamard2walsh {
	my @w; $#w = $#_;
	my @jw2jh = &jw2jh(scalar @_);
	my $i; for ($i=$[; $i<=$#_; $i++) { $w[$i] = $_[$jw2jh[$i]]; }
	return @w;
}

# ---------------------- EXPORT_OK stuff ---------------------------

sub biggest { my $k = shift @_; my @weeded = @_;
	my $smallest;
	if ($k <= 0) {
		my $tot = 0.0; foreach (@weeded) { $tot += abs $_; }
		$smallest = $tot / scalar @weeded;
	} else {
		my @sorted = sort { abs $b <=> abs $a } @weeded;
		$smallest = abs $sorted[$[-1+$k];
	}
	foreach (@weeded) { if (abs $_ < $smallest) { $_=0.0; } }
	return @weeded;
}
sub sublist { my ($aref, $offset, $length) = @_;
	if (ref $aref ne 'ARRAY') {
	warn "Math::WalshTransform::sublist 1st arg must be array ref\n"; return ();
	}
	my $first;
	if ($offset<0) { $first = $#{$aref}+$offset+1; } else { $first=$offset; }
	my $last;
	if (! defined $length) { $last = $#{$aref};
	} elsif ($length < 0) { $last = $#{$aref} + $length;
	} else { $last = $first + $length - 1;
	}
	# warn "offset=$offset length=$length first=$first last=$last\n";
	my @sublist = (); my $i;
	for ($i=$first; $i<=$last; $i++) { push @sublist, ${$aref}[$i]; }
	return @sublist;
}
# sub distance { my ($xref, $yref) = @_;  # Euclidian metric
# 	if (ref $xref ne 'ARRAY') {
# 		warn "Math::WalshTransform::distance 1st arg must be array ref\n";
# 		return undef;
# 	} elsif (ref $yref ne 'ARRAY') {
# 		warn "Math::WalshTransform::distance 2nd arg must be array ref\n";
# 		return undef;
# 	}
# 	my $distance = 0.0; my $i; my $diff;
# 	for ($i=$[; $i<= $#$xref; $i++) {
# 		$diff = ${$xref}[$i] - ${$yref}[$i];
# 		$distance += $diff * $diff;
# 	}
# 	return sqrt $distance;
# }
# sub size {my $sum=0.0; foreach (@_) {$sum+=$_*$_;} return sqrt $sum;}
# sub normalise {
#	my $size = &size(@_);
#	my @normalised = ();
#	foreach (@_) { push @normalised, $_/$size; }
#	return @normalised;
#}
sub average {
	my $i = $[; my $j; my @sum;
	foreach (@_) {
		if (ref $_ ne 'ARRAY') {
			warn "Math::WalshTransform::average argument $i must be array ref\n";
			return undef;
		}
		for ($j=$[; $j<=$#$_; $j++) { $sum[$j] += ${$_}[$j]; }
		$i++;
	}
	foreach (@sum) { $_ /= ($i-$[); }
	return @sum;
}
sub arr2txt { # neat printing of arrays for debug use
	my @txt; foreach (@_) { push @txt, sprintf('%g',$_); }
	return join (' ',@txt)."\n";
}

# ---------------------- infrastructure ---------------------------

sub jw2jh { my $n = shift;
	if ($n < 16) {
		if ($n == 8) { return (0,4,6,2,3,7,5,1); }
		if ($n == 4) { return (0,2,3,1); }
		if ($n == 2) { return (0,1); }
		if ($n == 1) { return (0); }
		warn "jw2jh: n=$n must be power of 2\n"; return ();
	}
	my @half; foreach (&jw2jh($n/2)) { push @half, $_<<1; }
	my @whole = @half; foreach (reverse @half) { push @whole, $_|1; }
	return @whole;
}

my $flag = 0;
sub gaussn {	my $standdev = $_[$[];
	# returns normal distribution around 0.0 by the Box-Muller rules
	if (! $flag) {
		$a = sqrt(-2.0 * log(rand));
		$b = 6.28318531 * rand;
		$flag = 1;
		return ($standdev * $a * sin($b));
	} else {
		$flag = 0;
		return ($standdev * $a * cos($b));
	}
}
1;

__END__