-----------------------------------------------------------------------------
-- |
-- Module      :  Numeric.Random.Spectrum.Pink
-- Copyright   :  (c) Matthew Donadio 2003
-- License     :  GPL
--
-- Maintainer  :  m.p.donadio@ieee.org
-- Stability   :  experimental
-- Portability :  portable
--
-- Functions for pinking noise
--
-- <http://www.firstpr.com.au/dsp/pink-noise/>
--
-----------------------------------------------------------------------------

module Numeric.Random.Spectrum.Pink (kellet,
				     voss) where

import DSP.Basic (downsample, upsampleAndHold)
import Data.List (tails)

-------------------------------------------------------------------------------

-- rb-j filter

-- pole            zero
-- ----            ----
-- 0.99572754      0.98443604
-- 0.94790649      0.83392334
-- 0.53567505      0.07568359

-------------------------------------------------------------------------------

-- | Kellet's filter

-- b0 = 0.99886 * b0 + white * 0.0555179;
-- b1 = 0.99332 * b1 + white * 0.0750759;
-- b2 = 0.96900 * b2 + white * 0.1538520;
-- b3 = 0.86650 * b3 + white * 0.3104856;
-- b4 = 0.55000 * b4 + white * 0.5329522;
-- b5 = -0.7616 * b5 - white * 0.0168980;
-- pink = b0 + b1 + b2 + b3 + b4 + b5 + b6 + white * 0.5362;
-- b6 = white * 0.115926;

kellet :: [Double] -- ^ noise
       -> [Double] -- ^ pinked noise

kellet :: [Double] -> [Double]
kellet [Double]
w = forall {t}.
Fractional t =>
[t] -> t -> t -> t -> t -> t -> t -> t -> [t]
kellet' [Double]
w Double
0 Double
0 Double
0 Double
0 Double
0 Double
0 Double
0
    where kellet' :: [t] -> t -> t -> t -> t -> t -> t -> t -> [t]
kellet' []         t
_  t
_  t
_  t
_  t
_  t
_  t
_  = []
          kellet' (t
white:[t]
ws) t
b0 t
b1 t
b2 t
b3 t
b4 t
b5 t
b6 = t
pink forall a. a -> [a] -> [a]
: [t] -> t -> t -> t -> t -> t -> t -> t -> [t]
kellet' [t]
ws t
b0' t
b1' t
b2' t
b3' t
b4' t
b5' t
b6'
	      where b0' :: t
b0' = t
0.99886 forall a. Num a => a -> a -> a
* t
b0 forall a. Num a => a -> a -> a
+ t
white forall a. Num a => a -> a -> a
* t
0.0555179
		    b1' :: t
b1' = t
0.99332 forall a. Num a => a -> a -> a
* t
b1 forall a. Num a => a -> a -> a
+ t
white forall a. Num a => a -> a -> a
* t
0.0750759
		    b2' :: t
b2' = t
0.96900 forall a. Num a => a -> a -> a
* t
b2 forall a. Num a => a -> a -> a
+ t
white forall a. Num a => a -> a -> a
* t
0.1538520
		    b3' :: t
b3' = t
0.86650 forall a. Num a => a -> a -> a
* t
b3 forall a. Num a => a -> a -> a
+ t
white forall a. Num a => a -> a -> a
* t
0.3104856
		    b4' :: t
b4' = t
0.55000 forall a. Num a => a -> a -> a
* t
b4 forall a. Num a => a -> a -> a
+ t
white forall a. Num a => a -> a -> a
* t
0.5329522
		    b5' :: t
b5' = -t
0.7616 forall a. Num a => a -> a -> a
* t
b5 forall a. Num a => a -> a -> a
- t
white forall a. Num a => a -> a -> a
* t
0.0168980
		    pink :: t
pink = t
b0 forall a. Num a => a -> a -> a
+ t
b1 forall a. Num a => a -> a -> a
+ t
b2 forall a. Num a => a -> a -> a
+ t
b3 forall a. Num a => a -> a -> a
+ t
b4 forall a. Num a => a -> a -> a
+ t
b5 forall a. Num a => a -> a -> a
+ t
b6 forall a. Num a => a -> a -> a
+ t
white forall a. Num a => a -> a -> a
* t
0.5362
		    b6' :: t
b6' = t
white forall a. Num a => a -> a -> a
* t
0.115926

-------------------------------------------------------------------------------

-- voss algorithm

add :: Num a => [[a]] -> [a]
add :: forall a. Num a => [[a]] -> [a]
add = forall (t :: * -> *) a. Foldable t => (a -> a -> a) -> t a -> a
foldl1 (forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith forall a. Num a => a -> a -> a
(+))

split :: Int -> [a] -> [[a]]
split :: forall a. Int -> [a] -> [[a]]
split Int
n = forall a. Int -> [a] -> [a]
take Int
n forall b c a. (b -> c) -> (a -> b) -> a -> c
. forall a b. (a -> b) -> [a] -> [b]
map (forall a. Int -> [a] -> [a]
downsample Int
n) forall b c a. (b -> c) -> (a -> b) -> a -> c
. (forall a. [a] -> [a] -> [a]
++forall a. a -> [a]
repeat []) forall b c a. (b -> c) -> (a -> b) -> a -> c
. forall a. [a] -> [[a]]
tails


mkOctaves :: [[a]] -> [[a]]
mkOctaves :: forall a. [[a]] -> [[a]]
mkOctaves = forall {a}. Int -> [[a]] -> [[a]]
mkOctaves' Int
1
    where mkOctaves' :: Int -> [[a]] -> [[a]]
mkOctaves' Int
_ []       = []
	  mkOctaves' Int
n ([a]
xs:[[a]]
xss) = forall a. Int -> [a] -> [a]
upsampleAndHold Int
n [a]
xs forall a. a -> [a] -> [a]
: Int -> [[a]] -> [[a]]
mkOctaves' (Int
2forall a. Num a => a -> a -> a
*Int
n) [[a]]
xss

-- | Voss's algorithm
--
-- UNTESTED, but the algorithm looks like it is working based on my hand
-- tests.

-- TODO: Since the input noise is consumed in different speed for the octaves
-- the function will certainly leak memory.

voss :: Int      -- ^ number of octaves to sum
     -> [Double] -- ^ noise
     -> [Double] -- ^ pinked noise

voss :: Int -> [Double] -> [Double]
voss Int
n [Double]
w = forall a. Num a => [[a]] -> [a]
add forall a b. (a -> b) -> a -> b
$ forall a. [[a]] -> [[a]]
mkOctaves forall a b. (a -> b) -> a -> b
$ forall a. Int -> [a] -> [[a]]
split Int
n [Double]
w

-------------------------------------------------------------------------------

-- voss-mccartney algorithm

-------------------------------------------------------------------------------

-- vm w =