{-# OPTIONS -fno-implicit-prelude #-} {- | Copyright : (c) Henning Thielemann 2006 License : GPL Maintainer : synthesizer@henning-thielemann.de Stability : provisional Portability : requires multi-parameter type classes -} module Synthesizer.Plain.Filter.NonRecursive where import qualified Synthesizer.Plain.Control as Ctrl import qualified Synthesizer.Plain.Signal as Sig import qualified Algebra.Transcendental as Trans import qualified Algebra.Module as Module import qualified Algebra.RealField as RealField import qualified Algebra.Field as Field import qualified Algebra.Ring as Ring import qualified Algebra.Additive as Additive import Algebra.Module(linearComb, (*>)) import Synthesizer.Utility (nest, ) import Data.List (tails) -- import Control.Monad.State (StateT) -- import Control.Monad.Writer (WriterT) import PreludeBase import NumericPrelude {- * Envelope application -} amplify :: (Ring.C a) => a -> Sig.T a -> Sig.T a amplify v = map (v*) amplifyVector :: (Module.C a v) => a -> Sig.T v -> Sig.T v amplifyVector = (*>) envelope :: (Ring.C a) => Sig.T a {-^ the envelope -} -> Sig.T a {-^ the signal to be enveloped -} -> Sig.T a envelope = zipWith (*) envelopeVector :: (Module.C a v) => Sig.T a {-^ the envelope -} -> Sig.T v {-^ the signal to be enveloped -} -> Sig.T v envelopeVector = zipWith (*>) fadeInOut :: (Field.C a) => Int -> Int -> Int -> Sig.T a -> Sig.T a fadeInOut tIn tHold tOut xs = let leadIn = take tIn $ Ctrl.linearMultiscale ( recip (fromIntegral tIn)) 0 leadOut = take tOut $ Ctrl.linearMultiscale (- recip (fromIntegral tOut)) 1 (partIn, partHoldOut) = splitAt tIn xs (partHold, partOut) = splitAt tHold partHoldOut in envelope leadIn partIn ++ partHold ++ envelope leadOut partOut fadeInOutAlt :: (Field.C a) => Int -> Int -> Int -> Sig.T a -> Sig.T a fadeInOutAlt tIn tHold tOut = zipWith id ((map (\x y -> y * fromIntegral x / fromIntegral tIn) [0..tIn-1]) ++ (replicate tHold id) ++ (map (\x y -> y * fromIntegral x / fromIntegral tOut) [tOut-1,tOut-2..0])) {- * Shift -} delay :: Additive.C y => Int -> Sig.T y -> Sig.T y delay = delayPad zero delayPad :: y -> Int -> Sig.T y -> Sig.T y delayPad z n = if n<0 then drop (negate n) else (replicate n z ++) {- * Smoothing -} {-| Unmodulated non-recursive filter -} generic :: Module.C a v => Sig.T a -> Sig.T v -> Sig.T v generic m x = let mr = reverse m xp = delay (pred (length m)) x in map (linearComb mr) (init (tails xp)) {-| Unmodulated non-recursive filter Output has same length as the input. It is elegant but leaks memory. -} genericAlt :: Module.C a v => Sig.T a -> Sig.T v -> Sig.T v genericAlt m x = map (linearComb m) (tail (scanl (flip (:)) [] x)) propGeneric :: (Module.C a v, Eq v) => Sig.T a -> Sig.T v -> Bool propGeneric m x = -- generic m x == genericAlt m x and $ zipWith (==) (generic m x) (genericAlt m x) {- | @eps@ is the threshold relatively to the maximum. That is, if the gaussian falls below @eps * gaussian 0@, then the function truncated. -} gaussian :: (Trans.C a, RealField.C a, Module.C a v) => a -> a -> a -> Sig.T v -> Sig.T v gaussian eps ratio freq = let var = ratioFreqToVariance ratio freq area = var * sqrt (2*pi) gau t = exp (-(t/var)^2/2) / area width = ceiling (var * sqrt (-2 * log eps)) -- inverse gau gauSmp = map (gau . fromIntegral) [-width .. width] in drop width . generic gauSmp {- GNUPlot.plotList [] (take 1000 $ gaussian 0.001 0.5 0.04 (Filter.Test.chirp 5000) :: [Double]) The filtered chirp must have amplitude 0.5 at 400 (0.04*10000). -} {- We want to approximate a Gaussian by a binomial filter. The latter one can be implemented by a convolutional power. However we still require a number of operations per sample which is proportional to the variance. -} binomial :: (Trans.C a, RealField.C a, Module.C a v) => a -> a -> Sig.T v -> Sig.T v binomial ratio freq = let width = ceiling (2 * ratioFreqToVariance ratio freq ^ 2) in drop width . nest (2*width) ((asTypeOf 0.5 freq *>) . binomial1) {- exp (-(t/var)^2/2) / area *> cis (2*pi*f*t) == exp (-(t/var)^2/2 +: 2*pi*f*t) / area == exp ((-t^2 +: 2*var^2*2*pi*f*t) / (2*var^2)) / area == exp ((t^2 - i*2*var^2*2*pi*f*t) / (-2*var^2)) / area == exp (((t^2 - i*var^2*2*pi*f)^2 + (var^2*2*pi*f)^2) / (-2*var^2)) / area == exp (((t^2 - i*var^2*2*pi*f)^2 / (-2*var^2) - (var*2*pi*f)^2/2)) / area sumMap (\t -> exp (-(t/var)^2/2) / area *> cis (2*pi*f*t)) [-infinity..infinity] ~ sumMap (\t -> exp (-(t/var)^2/2)) [-infinity..infinity] * exp (-(var*2*pi*f)^2/2) / area = exp (-(var*2*pi*f)^2/2) -} {- | Compute the variance of the Gaussian such that its Fourier transform has value @ratio@ at frequency @freq@. -} ratioFreqToVariance :: (Trans.C a) => a -> a -> a ratioFreqToVariance ratio freq = sqrt (-2 * log ratio) / (2*pi*freq) -- inverse of the fourier transformed gaussian binomial1 :: (Additive.C v) => Sig.T v -> Sig.T v binomial1 xt@(x:xs) = x : (xs + xt) binomial1 [] = [] {- | Moving (uniformly weighted) average in the most trivial form. This is very slow and needs about @n * length x@ operations. -} sums :: (Additive.C v) => Int -> Sig.T v -> Sig.T v sums n = map (sum . take n) . init . tails sumsDownsample2 :: (Additive.C v) => Sig.T v -> Sig.T v sumsDownsample2 (x0:x1:xs) = (x0+x1) : sumsDownsample2 xs sumsDownsample2 xs = xs downsample2 :: Sig.T a -> Sig.T a downsample2 (x0:_:xs) = x0 : downsample2 xs downsample2 xs = xs {- | Given a list of numbers and a list of sums of (2*k) of successive summands, compute a list of the sums of (2*k+1) or (2*k+2) summands. Eample for 2*k+1 @ [0+1+2+3, 2+3+4+5, 4+5+6+7, ...] -> [0+1+2+3+4, 1+2+3+4+5, 2+3+4+5+6, 3+4+5+6+7, 4+5+6+7+8, ...] @ Example for 2*k+2 @ [0+1+2+3, 2+3+4+5, 4+5+6+7, ...] -> [0+1+2+3+4+5, 1+2+3+4+5+6, 2+3+4+5+6+7, 3+4+5+6+7+8, 4+5+6+7+8+9, ...] @ -} sumsUpsampleOdd :: (Additive.C v) => Int -> Sig.T v -> Sig.T v -> Sig.T v sumsUpsampleOdd n {- 2*k -} xs ss = let xs2k = drop n xs in (head ss + head xs2k) : concat (zipWith3 (\s x0 x2k -> [x0+s, s+x2k]) (tail ss) (downsample2 (tail xs)) (tail (downsample2 xs2k))) sumsUpsampleEven :: (Additive.C v) => Int -> Sig.T v -> Sig.T v -> Sig.T v sumsUpsampleEven n {- 2*k -} xs ss = sumsUpsampleOdd (n+1) xs (zipWith (+) ss (downsample2 (drop n xs))) sumsPyramid :: (Additive.C v) => Int -> Sig.T v -> Sig.T v sumsPyramid = let aux 1 ys = ys aux 2 ys = ys + tail ys -- binomial1 aux m ys = let ysd = sumsDownsample2 ys in if even m then sumsUpsampleEven (m-2) ys (aux (div (m-2) 2) ysd) else sumsUpsampleOdd (m-1) ys (aux (div (m-1) 2) ysd) in aux {- *Synthesizer.Plain.Filter.NonRecursive> movingAverageModulated 10 (replicate 10 (3::Double) ++ [1.1,2.2,2.6,0.7,0.1,0.1]) (repeat (1::Double)) [0.5,0.6666666666666666,0.8333333333333333,1.0,1.0,1.0,1.0,1.0,1.0,1.0,0.9999999999999999,1.0,0.9999999999999998,0.999999999999999,0.9999999999999942,0.9999999999999942] -} {- * Filter operators from calculus -} {- | Forward difference quotient. Shortens the signal by one. Inverts 'Synthesizer.Plain.Filter.Recursive.Integration.run' in the sense that @differentiate (zero : integrate x) == x@. The signal is shifted by a half time unit. -} differentiate :: Additive.C v => Sig.T v -> Sig.T v differentiate x = zipWith subtract x (tail x) {- | Central difference quotient. Shortens the signal by two elements, and shifts the signal by one element. (Which can be fixed by prepending an appropriate value.) For linear functions this will yield essentially the same result as 'differentiate'. You obtain the result of 'differentiateCenter' if you smooth the one of 'differentiate' by averaging pairs of adjacent values. ToDo: Vector variant -} differentiateCenter :: Field.C v => Sig.T v -> Sig.T v differentiateCenter x = map ((1/2)*) $ zipWith subtract x (tail (tail x)) {- | Second derivative. It is @differentiate2 == differentiate . differentiate@ but 'differentiate2' should be faster. -} differentiate2 :: Additive.C v => Sig.T v -> Sig.T v differentiate2 xs0 = let xs1 = tail xs0 xs2 = tail xs1 in zipWith3 (\x0 x1 x2 -> x0+x2-(x1+x1)) xs0 xs1 xs2