view release on metacpan or  search on metacpan

lua/WalshTransform.lua  view on Meta::CPAN

-- ----------------------------------------------------------------- --
--      This Lua5 module is Copyright (c) 2010, Peter J Billam       --
--                        www.pjb.com.au                             --
--                                                                   --
--   This module is free software; you can redistribute it and/or    --
--          modify it under the same terms as Lua5 itself.           --
-- ----------------------------------------------------------------- --
local M = {} -- public interface
M.Version = '1.17'
M.VersionDate = '12aug2010'

-- Example usage:
-- local WT = require 'WalshTransform'
-- H = WT.fht(a)

--------------------- infrastructure ----------------------
local function warn(str)
    io.stderr:write(str,'\n')
end
local function die(str)
    io.stderr:write(str,'\n')
    os.exit(1)
end
local function jw2jh(n)
	if (n < 16) then
		if n == 8 then return {0,4,6,2,3,7,5,1}; end
		if n == 4 then return {0,2,3,1}; end
		if n == 2 then return {0,1}; end
		if n == 1 then return {0}; end
		warn "jw2jh: n=$n must be power of 2\n"; return false;
	end
	local half_n = n / 2;
	local half = jw2jh(half_n)
	local whole = {}
	local i = 1; while i <= half_n do
		whole[i] = 2 * half[i]
		i = i + 1
	end
	i = half_n; local j = half_n + 1
	while i >= 1 do
		whole[j] = 1 + whole[i]
		i = i - 1
		j = j + 1
	end
	return whole
end
-------------------------------------------------------

function M.fht(a)
	local i; local mr = {}; for i=1,#a do mr[i] = a[i] end
	local k = 1
	local n = #mr
	local l = n
	local nl; local nk; local j;
	while true do
		i = 0; l = l/2;
		for nl=1,l do
			for nk=1,k do
				i = i + 1; j = i + k;
				mr[i] = (mr[i] + mr[j])/2;
				mr[j] =  mr[i] - mr[j];
			end
			i = j;
		end
		k = 2*k;
		if k >= n then break end
	end
	if k == n then
		return mr
	else
		warn("WalshTransform::fht \$n = $n but must be power of 2\n")
		return false
	end

lua/WalshTransform.lua  view on Meta::CPAN

function M.size(a)
	local sum=0.0; local i; for i=1,#a do sum = sum + a[i]*a[i] end
	return math.sqrt(sum)
end

function M.product(a,b)
	if type(a) ~= 'table' then
		warn("WalshTransform.product 1st arg must be a table\n")
		return nil
	end
	if type(b) ~= 'table' then
		warn("WalshTransform.product 2nd arg must be a table\n")
		return nil
	end
	if #a ~= #b then
		warn("WalshTransform.product the 2 args must be the same size\n")
		return nil
	end
	local c = {}; local i;
	for i=1,#a do c[i] = a[i]*b[i] end
	return c
end

function M.normalise(a)
	if type(a) ~= 'table' then
		warn("WalshTransform.product 1st arg must be a table\n")
		return nil
	end
	local s = M.size(a);  local normalised = {};  local i
	for i = 1, #a do normalised[i] = a[i]/s end
	return normalised
end

function M.average (...)
	local i = 1; local j; local sum = {}; local n = 0
	for i,a in ipairs({...}) do
		if type(a) ~= 'table' then
			warn("WalshTransform.average arg must be a table\n")
			return nil
		end
		for j=1, #a do sum[j] = (sum[j] or 0) + a[j] end
		n = n + 1
	end
	for j=1,#sum do sum[j] = sum[j] / n end
	return sum
end

return M

--[[

=pod

=head1 NAME

WalshTransform.lua - Fast Hadamard and Walsh Transforms

=head1 SYNOPSIS

 local WT = require 'WalshTransform' -- name doesn't matter of course
 local f = {1.618, 2.718, 3.142, 4.669}  -- must be power-of-two long
 local FH1 = WT.fht(f)   -- Hadamard transform
 local fh1 = WT.fhtinv(FH1)
 -- or
 local FW2 = WT.fwt(f)   -- Walsh transform
 local fw2 = WT.fwtinv(FW2)
 local FH2 = WT.walsh2hadamard(FW2)

 local PS  = WT.power_spectrum(f)

 local whats_going_on = WT.biggest(9,WT.fwt(WT.sublist(time_series,-16)))
 local EVENT1 = WT.fwt(WT.sublist(time_series,478,16))
 local EVENT2 = WT.fwt(WT.sublist(time_series,2316,16))
 local EVENT3 = WT.fwt(WT.sublist(time_series,3261,16))
 EVENT1[1]=0.0; EVENT2[1]=0.0; EVENT3[1]=0.0 -- ignore constant
 EVENT1 = WT.normalise(EVENT1) -- ignore scale
 EVENT2 = WT.normalise(EVENT2)
 EVENT3 = WT.normalise(EVENT3)
 local TYPICAL_EVENT = WT.average(EVENT1, EVENT2, EVENT3)
 ...
 local NOW = WT.fwt(WT.sublist(time_series,-16))
 NOW[1] = 0.0
 NOW = WT.normalise(NOW)
 if WT.distance(NOW, TYPICAL_EVENT) < .28 then get_worried() end

=head1 DESCRIPTION

These routines implement fast Hadamard and Walsh Transforms
and their inverse transforms.

Also included are routines
for converting a Hadamard to a Walsh transform and vice versa,
for calculating the Logical Convolution of two lists,
or the Logical Autocorrelation of a list,
and for calculating the Walsh Power Spectrum -
in short, almost everything you ever wanted to do with a Walsh Transform.

Intelligible speech can be reconstructed by transforming
blocks of, say, 64 samples, deleting all but the several largest
transform components, and inverse-transforming;
in other words, these transforms extract from a time-series
the most significant things that are going on.
They should be useful for B<noticing important things>,
for example in software that monitors time-series data such as
system or network administration data, share-price, currency,
ecological, opinion poll, process management data, and so on.

As from version 1.10, Math::WalshTransform uses C routines to perform
the transforms. This runs 25 to 30 times faster than previous versions.

Not yet included are multi-dimensional Hadamard and Walsh Transforms,
conversion between Logical and Arithmetic Autocorrelation Functions,
or conversion between the Walsh Power Spectrum and the Fourier Power Spectrum.

Version 1.17

=head1 SUBROUTINES

Routines which take just one array as argument expect the array itself;
those which take more than one array expect a list of references.