packages feed

synthesizer-core 0.4.0.4 → 0.4.1

raw patch · 97 files changed

+3964/−444 lines, 97 filesdep +explicit-exceptiondep +timeitdep ~basedep ~containersdep ~numeric-preludenew-component:exe:fouriertest

Dependencies added: explicit-exception, timeit

Dependency ranges changed: base, containers, numeric-prelude, storablevector

Files

+ speedtest/Fourier.hs view
@@ -0,0 +1,94 @@+module Main where++import qualified Synthesizer.Storable.Signal as SigSt+import qualified Synthesizer.Generic.Fourier as Fourier+import qualified Synthesizer.Generic.Noise as NoiseG+import qualified Synthesizer.Generic.Signal as SigG+import qualified Synthesizer.Generic.Cut as CutG+import qualified Synthesizer.State.Noise as NoiseS+import qualified Synthesizer.State.Signal as SigS++import qualified Data.StorableVector as SV++import qualified Algebra.Ring as Ring+import qualified Number.Complex as NPComplex++import System.TimeIt (timeIt, )++import NumericPrelude.Numeric+import NumericPrelude.Base+import Prelude ()++++test0 :: IO ()+test0 =+   SigSt.writeFile "fouriertest.f64" $+   SigG.take 65536 $+   (NoiseG.white SigG.defaultLazySize :: SigSt.T Double)++test1 :: IO ()+test1 =+   SigSt.writeFile "fouriertest.f64" $+   SigG.fromState SigG.defaultLazySize $+   SigS.take 65536 $+   SigS.map (NPComplex.+: 0) $+   (NoiseS.white :: SigS.T Double)++test2 :: Int -> IO ()+test2 n =+   writeFile "fouriertest.cache" $+   show $ Fourier.cacheBackward $+   (\sig ->+      SigG.fromState SigG.defaultLazySize sig ::+         SigSt.T (NPComplex.T Double)) $+   SigS.take n $+   SigS.map (NPComplex.+: 0) $+   NoiseS.white++test3 :: Int -> IO ()+test3 n =+   let sig :: SigSt.T (NPComplex.T Double)+       sig =+          SigG.fromState SigG.defaultLazySize $+          SigS.take n $+          SigS.map (NPComplex.+: 0) $+          NoiseS.white+       cache =+          Fourier.cacheBackward sig+   in  do timeIt $ writeFile "fouriertest.cache" $ show cache+          timeIt $ SigSt.writeFile "fouriertest.f64" $+             Fourier.transformWithCache cache sig++test4 :: Int -> IO ()+test4 n =+   let sig :: SV.Vector (NPComplex.T Double)+       sig =+          SigS.toStrictStorableSignal n $+          SigS.take n $+          SigS.map (NPComplex.+: 0) $+          NoiseS.white+       cache =+          Fourier.cacheBackward sig+   in  do -- timeIt $ writeFile "fouriertest.cache" $ show cache+          timeIt $ SV.writeFile "fouriertest.f64" $+             Fourier.transformWithCache cache sig+++main :: IO ()+main =+--   timeIt $ test2 (4096*3+1)+--   test4 (4096*3+1)+   sequence_ $+   timeIt test0 : timeIt test1 :+   map test4+      (16384 : (4096*3) : (4096*3+1) : 11025 :+       (3^9) : (5^6) : (7^5) :+       (6^6) : (3*5*7*11*13) :+       [])+{-+      (65536 : 65537 : 44100 :+       (3^10) : (5^7) : (7^5) :+       (6^6) : (2*3*5*7*11*13) :+       [])+-}
src/Synthesizer/Basic/ComplexModule.hs view
@@ -6,11 +6,10 @@ import qualified Number.Complex as Complex import qualified Algebra.Module as Module import Number.Complex ((+:), )-import Algebra.Module ((*>), ) -import qualified Prelude as P--- import NumericPrelude.Base import NumericPrelude.Numeric+-- import NumericPrelude.Base+import Prelude ()   {-# INLINE scale #-}
src/Synthesizer/Basic/Distortion.hs view
@@ -16,9 +16,8 @@ import qualified Algebra.Transcendental        as Trans import qualified Algebra.RealField             as RealField import qualified Algebra.Field                 as Field-import qualified Algebra.RealRing                  as RealRing+import qualified Algebra.RealRing              as RealRing import qualified Algebra.Ring                  as Ring-import qualified Algebra.Additive              as Additive  import Data.Ord.HT (limit, )
src/Synthesizer/Basic/DistortionControlled.hs view
@@ -10,9 +10,7 @@ import qualified Algebra.Transcendental        as Trans import qualified Algebra.RealField             as RealField import qualified Algebra.Field                 as Field-import qualified Algebra.RealRing                  as RealRing-import qualified Algebra.Ring                  as Ring-import qualified Algebra.Additive              as Additive+import qualified Algebra.RealRing              as RealRing  import Data.Ord.HT (limit, )
+ src/Synthesizer/Basic/NumberTheory.hs view
@@ -0,0 +1,896 @@+{-# LANGUAGE NoImplicitPrelude #-}+{-+Some of these functions might be moved to NumericPrelude.++Wikipedia: (primitive) roots of unity modulo n+   (primitive) roots must be units and all units are (primitive) roots+   maximum possible order for primitive roots - Carmichael+   all possible orders: divisor of Carmichael (proof? statement already in Carmichael-function-article)+   sum of primitive roots that vanishes+   order of primitive root is a divisor of each possible exponent+      proof with GCD and diophantine in exponent+   check for primitive root: fast exponentiation,+      primitivity: check exponents that are prime divisors+   how to find a primitive root: just try+   sum of powers of a primitive root is zero+   number of primitive roots of given order+      g(n,k) > 0 if k|lambda(n)+      g(n,k) = 0 else+      g(n,1) = 1+      g(4,2) = 1+      g(2^n,2) = 3 for n>=3  ((-1) is always a square root of 1)+      g(2^n,2^k) = 2^k for k>=2 && k<n-1+      g(n,2) = 1 for n>=3 and n in OEIS:A033948+      sum(g(n,k), k\in\N) = phi(n)+      There are only a few patterns that occur as rows of g,+      but a row of g (i.e. g(n)) does functionally depend on+      either lambda(n) nor phi(n)+      lambda(14) = 6   nozeros(g(14)) = [1,1,2,2]   (f ~ [1,2,3,6])+      lambda(21) = 6   nozeros(g(21)) = [1,3,2,6]   (f ~ [1,4,3,12])+      phi(13) = 12   nozeros(g(13)) = [1,1,2,2,2,4]   (f ~ [1,2,3,4,6,12])+      phi(21) = 12   nozeros(g(21)) = [1,3,2,6]       (f ~ [1,4,3,12])+      However length(nozeros(f(n))) = numberofdivisors(lambda(n))+      numberofdivisors=A000005+   number of roots of given order+      easier to compute+      k|m => f(n,k) | f(n,m)+      g(n,k) = f(n,k) - sum(f(n,d), d|k and k/d prime) + ...+         inclusion-exclusion-principle+      better to write the other round:+      f(n,k) = sum(g(n,d), d|k) - this is Dirichlet convolution+      RUNM says f(n,k) is multiplicative+         list it in multiplicative function+      f(n,1) = 1 for n>=2+      f(n,lambda(n)) = phi(n)+      f(n,a·b) = f(n,a)·f(n,b) if a and b are coprime (conjecture)+      f(n,lcm(a,b)) = lcm(f(n,a),f(n,b)) (conjecture)+      If this conjecture is true, then we only need to know f(n,p^i).+      The following conjecture is wrong:+         for prime p it is   f(n,p^i) = gcd(lambda(n),p^i)+      counterexamples+         f(8,2) = 4, lambda(8)=2+         f(63,3) = 9, lambda(63)=6+         f(275,5) = 25, lambda(275)=20+         f(1247,7) = 49, lambda(1247)=84+      It seems to be:+         for prime p it is   f(n,p^i) = p^j for some j+   How to find a modulus where there is a primitive root of order o?+      just try numbers from the sequence o+1, 2*o+1, 3*o+1+      Because of [[Dirichlet's theorem on arithmetic progressions]]+      you will somewhen find a prime p,+      and its Carmichael value is p-1, which is a multiple of o.+      In this ring even 'o' is a unit.+   How to find a modulus where there are primitive roots of orders o1,..,ok?+      Just search for a modulus with roots of order lcm(o1,...,ok).+      The smallest such modulus should also be the smallest one+      that has primitive roots of orders o1,..,ok?+      Proof: If a ring has primitive roots of orders o1,..,ok+      then all orders divide the carmichael value of that ring,+      thus lcm(o1,...,ok) divides the carmichael value of that ring,+      thus there is a primitive root of order lcm(o1,...,ok).+-}+module Synthesizer.Basic.NumberTheory where++import qualified Synthesizer.State.Signal as SigS++import qualified Data.Set as Set+import qualified Data.Map as Map++import qualified Algebra.Ring as Ring+import qualified Algebra.Units as Units+import qualified Algebra.PrincipalIdealDomain as PID+import qualified Algebra.IntegralDomain as Integral+import qualified Algebra.ZeroTestable as ZeroTestable++import qualified Number.ResidueClass.Check as RC+import Number.ResidueClass.Check ((/:), )++import qualified Number.FixedPoint as FP+import Data.Bits (Bits, (.&.), (.|.), shiftR, )++import qualified Data.List.HT as ListHT+import Data.List (unfoldr, mapAccumL, genericDrop, genericSplitAt, )+import Data.Tuple.HT (mapFst, mapSnd, mapPair, swap, )+import Data.Maybe.HT (toMaybe, )++import Test.QuickCheck (Arbitrary(arbitrary), )++import NumericPrelude.Numeric+import NumericPrelude.Base+++{- |+The first pair member in @powerOfTwoFactors n@+is the maximum factor of @n@, that is a power of two.+-}+powerOfTwoFactors ::+   (Bits a, Integral.C a) => a -> (a, a)+powerOfTwoFactors n =+   let powerOfTwo = n .&. (-n)+   in  (powerOfTwo, div n powerOfTwo)+++{- |+List all factorizations of an odd number+where the first factor is at most the second factor+and the first factors are in descending order.+-}+fermatFactors :: Integer -> [(Integer,Integer)]+fermatFactors n =+   let root = FP.sqrt 1 n+   in  map (\(a,b) -> (b-a,b+a)) $+       mergeAndFilter+          (zip (scanl (+) n [1,3..]) [0 .. div (n-1) 2])+          (zip (scanl (+) (root*root) $ iterate (2+) (2*root+1)) [root..])++mergeAndFilter :: (Ord a) => [(a,b)] -> [(a,c)] -> [(b,c)]+mergeAndFilter ((a0,b):a0s) ((a1,c):a1s) =+   case compare a0 a1 of+      LT -> mergeAndFilter a0s ((a1,c):a1s)+      GT -> mergeAndFilter ((a0,b):a0s) a1s+      EQ -> (b,c) : mergeAndFilter a0s a1s+mergeAndFilter _ _ = []++++{- |+Argument must be a prime.+Usage of Set for efficient filtering of candidates seems to be overkill,+since the multiplicative generator seems to be small in most cases,+i.e. 2 or 3.+-}+multiplicativeGenerator :: Integer -> Integer+multiplicativeGenerator p =+   let search candidates =+          case Set.minView candidates of+             Nothing -> error $ show p ++ " is not an odd prime"+             Just (x,rest) ->+                case orbitSet $ orbit p x of+                   new ->+                      -- fromIntegral (Set.size new) == p-2+                      if new == Set.fromList [1..p-1]+                        then x+                        else search (Set.difference rest new)+   in  search (Set.fromList [2..p-1])+++newtype Order = Order {getOrder :: Integer}+   deriving (Show, Eq, Ord)++instance Arbitrary Order where+   arbitrary = fmap (Order . (1+) . abs) arbitrary++instance Enum Order where+   succ (Order n) = Order (n+1)+   pred (Order n) = Order (n-1)+   fromEnum (Order n) = fromEnum n+   toEnum n = Order (toEnum n)+   enumFrom (Order from) =+      map Order $ enumFrom from+   enumFromThen (Order from) (Order thn) =+      map Order $ enumFromThen from thn+   enumFromTo (Order from) (Order to) =+      map Order $ enumFromTo from to+   enumFromThenTo (Order from) (Order thn) (Order to) =+      map Order $ enumFromThenTo from thn to++countOrder :: [a] -> Order+countOrder = foldl (\o _ -> succ o) (Order 0)++dividesOrder :: Order -> Order -> Bool+dividesOrder (Order k) (Order n) =+   divides k n+++-- class Integral.C a => PrimitiveRoot a where+class PID.C a => PrimitiveRoot a where+   primitiveRootCandidates :: a -> [a]+   maximumOrderOfPrimitiveRootsOfUnity :: a -> Order++instance PrimitiveRoot Integer where+   primitiveRootCandidates modu = [1 .. modu-1]+   maximumOrderOfPrimitiveRootsOfUnity =+      maximumOrderOfPrimitiveRootsOfUnityInteger++{-+For 'ordersOfPrimitiveRootsOfUnityInteger'+and the connection to Euler's totient function+we have chosen to have++> primitiveRootsOfUnity m 1 == [1].+-}+primitiveRootsOfUnity ::+   (PrimitiveRoot a, Eq a) => a -> Order -> [a]+primitiveRootsOfUnity =+   primitiveRootsOfUnityPower++{-+Verifying that a ring has no primitive root of the wanted order+takes a long time if we do it by exhaustive search.+In the case of a=Integer we could first check,+whether the considered residue ring has a primitive root of wanted order+using the Carmichael function.+We could certainly count the number of primitive roots+and stop searching if we reach that count.+-}+primitiveRootsOfUnityPower ::+   (PrimitiveRoot a, Eq a) => a -> Order -> [a]+primitiveRootsOfUnityPower modu (Order order) =+   let greatDivisors = map (div order) $ uniquePrimeFactors order+   in  filter+          (\n ->+             let pow y = RC.representative $ (n /: modu) ^ y+             in  PID.coprime n modu+                 &&+                 pow order == one+                 &&+                 all (\y -> pow y /= one) greatDivisors) $+       primitiveRootCandidates modu++primitiveRootsOfUnityNaive ::+   (PrimitiveRoot a, Eq a) => a -> Order -> [a]+primitiveRootsOfUnityNaive _ (Order 0) = []+primitiveRootsOfUnityNaive modu (Order expo) =+   filter+      (\n ->+         let (prefix,end:_) =+                genericSplitAt (expo-1) $ SigS.toList $ orbit modu n+         in  all (1/=) prefix && end==1) $+   primitiveRootCandidates modu++orbitSet :: Ord a => SigS.T a -> Set.Set a+orbitSet list =+   SigS.foldR+      (\new cont seen ->+         if Set.member new seen+           then seen+           else cont (Set.insert new seen))+      id list Set.empty++orbit :: (Integral.C a) => a -> a -> SigS.T a+orbit p x = SigS.iterate (\y -> mod (x*y) p) x+++{- |+Does not emit values in ascending order+and may return duplicates (e.g. primitiveRootsOfUnityFullOrbit 70000 10).+I hoped it would be faster than the other implementations+since it eliminates non-roots in large blocks.+However it seems that managing the root candidates in a Set+reduces performance significantly.++The idea:+Start with a seed that is a unit.+Compute its orbit until a one is reached.+Since it is a unit, we always encounter a one.+We do not need to check for non-unit seeds,+since (gcd modu seed) will be a divisor of all seed powers.+In the orbit all numbers are powers of each other.+Now finding the roots is a matter of solving+a Diophantine equation of the exponents.+In one such step all powers in an orbit are classified as roots or non-roots+and thus we can remove them all from the set of root candidates+and continue with the remaining candidates.+Duplicates can occur if a seed+in a later iteration is found again as power of another seed.+-}+primitiveRootsOfUnityFullOrbit ::+   (PrimitiveRoot a, Ord a) => a -> Order -> [a]+primitiveRootsOfUnityFullOrbit modu expo =+   let search candidates =+          flip fmap (Set.minView candidates) $ \(x,rest) ->+          mapSnd (Set.difference rest . Set.fromList) $+          primitiveRootsOfOrbit modu expo x+   in  concat $ unfoldr search $ Set.fromList $+       -- needed for modules with many divisors+       filter (PID.coprime modu) $+       primitiveRootCandidates modu++primitiveRootsOfUnityFullOrbitTest ::+   (PrimitiveRoot a, Ord a) => a -> Order -> [(a,[a])]+primitiveRootsOfUnityFullOrbitTest modu expo =+   let search candidates =+          flip fmap (Set.minView candidates) $ \(x,rest) ->+          mapPair ((,) x,+                   Set.difference rest . Set.fromList) $+          primitiveRootsOfOrbit modu expo x+   in  unfoldr search $ Set.fromList $+       filter (PID.coprime modu) $+       primitiveRootCandidates modu++primitiveRootsOfOrbit ::+   (PrimitiveRoot a, Ord a) => a -> Order -> a -> ([a], [a])+primitiveRootsOfOrbit modu (Order expo) x =+   let orb = (1:) $ takeWhile (1/=) $ iterate (\y -> mod (x*y) modu) x+       (Order orbitSize) = countOrder orb+   in  (if expo==0+          then []+          else+            {-+            size = length orb+            Search for m and k with 0<k and 0<m and m<size+            and expo*m = size*k+            such that for all l with 0<l and l<k+            m does not divide size*l.+            To this end we ask for every m+            for the smallest r such that size divides r*m.+            If r=expo then x^m is a primitive root of order expo.+            If r<expo then x^m has order smaller than expo.+            The searched r is div size (gcd size m).+            However expo = div size (gcd size m) implies,+            that expo is a divisor of size.+                expo = div size (gcd size m)+             => div size expo = gcd size m+                s = gcd size m+            We have to consider for m+            only the multiples of s.+            Then divide both sides of the equation by s, yielding+                1 = gcd expo m'+            -}+            case divMod orbitSize expo of+               (s,0) ->+                  map snd $ filter (PID.coprime expo . fst) $+                  zip+                     [0 .. expo-1]+                     -- (ListHT.sieve s $ orb)+                     (map head $ iterate (genericDrop s) orb)+               _ -> [],+        orb)+++hasPrimitiveRootOfUnityNaive ::+   (PrimitiveRoot a, Ord a) => a -> Order -> Bool+hasPrimitiveRootOfUnityNaive modu expo =+   any (dividesOrder expo . snd) $+   ordersOfPrimitiveRootsOfUnityTest modu++{-+This should be a maximum both with respect to magnitude and to divisibility.+-}+maximumOrderOfPrimitiveRootsOfUnityNaive ::+   (PrimitiveRoot a, Ord a) => a -> Order+maximumOrderOfPrimitiveRootsOfUnityNaive =+   foldl max (Order 1) . map snd . ordersOfPrimitiveRootsOfUnityTest++{- |+Computes a list of seeds and according maximum orders of roots of unity.+All divisors of those maximum orders are possible orders of roots of unity, too.+-}+ordersOfPrimitiveRootsOfUnityTest ::+   (PrimitiveRoot a, Ord a) => a -> [(a, Order)]+ordersOfPrimitiveRootsOfUnityTest modu =+   let search candidates =+          flip fmap (Set.minView candidates) $ \(x,rest) ->+          mapPair ((,) x,+                   Set.difference rest . Set.fromList) $+          orderOfOrbit modu x+   in  unfoldr search $ Set.fromList $+       filter (PID.coprime modu) $+       primitiveRootCandidates modu++{- |+modu and x must be coprime.+If they are not,+then none of the numbers in the orbit is a root of unity+and the function enters an infinite loop.+-}+orderOfOrbit ::+   (PrimitiveRoot a, Ord a) => a -> a -> (Order, [a])+orderOfOrbit modu x =+   let cyc = takeWhile (one/=) $ SigS.toList $ orbit modu x+   in  (succ $ countOrder cyc, cyc)+++{-+This test speeds up 'hasPrimitiveRootOfUnityNaive' considerably+by considering the prime factors of modu.+If modu is a prime, then the ring has a multiplicative generator,+i.e. a primitive root of unity of order modu-1.+That is, all primitive roots of unity are of an order that divides modu-1.+It seems that if modu is a power of a prime,+then the according ring has also multiplicative generator.+I think this is the reason for generalising the Rader transform+to signals of prime power length.+-}+hasPrimitiveRootOfUnityInteger ::+   Integer -> Order -> Bool+hasPrimitiveRootOfUnityInteger modu expo =+   dividesOrder expo $+   maximumOrderOfPrimitiveRootsOfUnityInteger modu++{-+Carmichael theorem:+If the integer residue rings with coprime moduli m0 and m1+have primitive roots of maximum order o0 and o1, respectively,+then the integer ring with modulus m0*m1+has maximum order (lcm o0 o1).+-}++{-+This is the Carmichael function.+OEIS-A002322+-}+maximumOrderOfPrimitiveRootsOfUnityInteger ::+   Integer -> Order+maximumOrderOfPrimitiveRootsOfUnityInteger =+   Order .+   lcmMulti .+   map+      (\(e,p) ->+         if p == 2 && e > 2+           then p^(e-2)+           else p^(e-1) * (p-1)) .+   map (mapFst fromIntegral) .+   primeFactors+++{-+The sum of the sub-lists should equal the Euler totient function values+(OEIS-A000010).+-}+ordersOfPrimitiveRootsOfUnityInteger :: [[Int]]+ordersOfPrimitiveRootsOfUnityInteger =+   flip map [1..] $ \modu ->+   let maxOrder = maximumOrderOfPrimitiveRootsOfUnity (modu::Integer)+   in  map (length . primitiveRootsOfUnityPower modu) $+--       filter (flip divides maxOrder) $+       [Order 1 .. maxOrder]++ordersOfRootsOfUnityInteger :: [[Int]]+ordersOfRootsOfUnityInteger =+   flip map [1..] $ \modu ->+   map (length . rootsOfUnityPower (modu::Integer)) $+   [Order 1 ..]+{-+mapM_ print $ map (\n -> (n, maximumOrderOfPrimitiveRootsOfUnityInteger (fromIntegral n), take 30 $ ordersOfRootsOfUnityInteger !! (n-1))) [2..30]++mapM_ print $ map (\n -> (n, maximumOrderOfPrimitiveRootsOfUnityInteger (fromIntegral n), let row = ordersOfRootsOfUnityInteger !! (n-1) in map (row!!) $ map pred $ take 10 $ iterate (2*) 1)) [2..30]+-}++ordersOfRootsOfUnityIntegerCondensed :: [[Int]]+ordersOfRootsOfUnityIntegerCondensed =+   flip map [1..] $ \modu ->+   let maxOrder = maximumOrderOfPrimitiveRootsOfUnity (modu::Integer)+   in  map (length . rootsOfUnityPower modu) $+--       filter (flip divides maxOrder) $+       [Order 1 .. maxOrder]++rootsOfUnityPower ::+   (PrimitiveRoot a, Eq a) => a -> Order -> [a]+rootsOfUnityPower modu (Order expo) =+   filter+      (\n ->+         PID.coprime n modu+         &&+         RC.representative ((n /: modu) ^ expo) == one) $+   primitiveRootCandidates modu++{-+Corollary from the Carmichael function properties:+If in Z_n there is a primitive root r of unity of order o,+then for every Z_{m \cdot n} there is also a primitive root of unity+with the same order.++Collary:+If in Z_n1 you have a primitive root of order o1,+and so on for Z_{n_k} and order ok,+then Z_{lcm(n1,...,nk)} has primitive roots for every of the order o1,...,on.++Conjecture:+If Z_n has a total number of m primitive roots of unity of order o,+then Z_{k*n} has at least m primitive roots of unity of order o.++Conjecture:+Precondition: In Z_n there is a primitive root of prime order o.+Claims:+a) There are natural numbers m and k with n = m*(k*o+1) or n = m*o.+b) The smallest such n is of the form k*o+1 with k>1.++Counterexample for a) and non-prime o: o=15, n=77+Counterexample for b) and non-prime o:+   o=20, n=25; o=27, n=81; o=54, n=81; o=55, n=121++Corollary from definition of Carmichael function:+For n>1, Z_{2^{n+2}} has primitive root of unity of order 2^n.+-}++{- |+Given an order find integer residue rings+where roots of unity of this order exist.+The way they are constructed also warrants,+that 'order' is a unit (i.e. invertible) in those rings.++The list is not exhaustive+but computes suggestions quickly.+The first found modulus seems to be smallest one that exist.+However, the first modulus is not the smallest one+among the ones that only have the wanted primitive root,+but where 'order' is not necessarily a unit.+E.g.++> ringsWithPrimitiveRootOfUnityAndUnit 840 == 2521 : 3361 : ...++but the smallest modulus is 1763.++Most of the numbers are primes.+For these the recursion property of the Carmichael function+immediately yields that there are primitive roots of unity of order 'order'.+-}+ringsWithPrimitiveRootOfUnityAndUnit :: Order -> [Integer]+ringsWithPrimitiveRootOfUnityAndUnit order@(Order k) =+   filter (flip hasPrimitiveRootOfUnityInteger order) $+   iterate (k+) 1+++ringsWithPrimitiveRootsOfUnityAndUnitsNaive :: [Order] -> [Integer] -> [Integer]+ringsWithPrimitiveRootsOfUnityAndUnitsNaive rootOrders units =+   filter+      (\n ->+         all (hasPrimitiveRootOfUnityInteger n) rootOrders &&+         all (PID.coprime n) units)+      [1..]+++{-+It would be nice to have the Omega monad here+in order to enumerate all rings.+-}+ringWithPrimitiveRootsOfUnityAndUnits :: [Order] -> [Integer] -> Integer+ringWithPrimitiveRootsOfUnityAndUnits rootOrders units =+   let p = lcmMulti units+   in  lcmMulti $+       map (head . filter (PID.coprime p) .+            ringsWithPrimitiveRootOfUnityAndUnit) $+       rootOrders++{-+Search for an appriopriate modulus+using the recursive definition of the Carmichael function.+We chose the prime factors of the Carmichael function argument+such that we get at least the prime factors in the function value that we need.++The modulus constructed this way is usually not the smallest possible+and it also does not provide that 'n' is a unit in the residue ring.+The simpler function 'ringsWithPrimitiveRootOfUnityAndUnit'+will usually produce a smaller modulus.+-}+ringWithPrimitiveRootsOfUnity :: Order -> Integer+ringWithPrimitiveRootsOfUnity (Order n) =+   case n of+      0 -> 2+      _ ->+         product $ map (uncurry ringPower) $ snd $+         mapAccumL+            (\factors (e,p) ->+               if Map.findWithDefault 0 p factors >= e+                 then (factors, (0,p))+                 else+                   if p==2+                     then+                       (factors,+                        case e of+                           0 -> (0,2)+                           1 -> (1,3)+                           2 -> (1,5)+                           _ -> (e+2, 2))+                     else+                       (Map.unionWith max factors $+                           Map.fromList $ map swap $ primeFactors $ p-1,+                        (e+1, p)))+            Map.empty $+         reverse $ primeFactors $ lcmMulti $+         n : map (subtract 1) (partialPrimes n)++lcmMulti :: (PID.C a) => [a] -> a+lcmMulti = foldl lcm one+++{- |+List all numbers that only contain prime factors 2 and 3 in ascending order.+OEIS:A003586+-}+numbers3Smooth :: [Integer]+numbers3Smooth =+   foldr+      (\(x0:x1:xs) ys -> x0 : x1 : ListHT.mergeBy (<=) xs ys)+      (error "numbers3Smooth: infinite list should not have an end") $+   iterate (map (3*)) $+   iterate (2*) 1++numbers3SmoothAlt :: [Integer]+numbers3SmoothAlt =+   unfoldr+      (fmap (\(m,rest) -> (m, Set.union rest $ Set.fromAscList [2*m,3*m])) .+       Set.minView) $+   Set.singleton 1++{-+OEIS:A051037+-}+numbers5Smooth :: [Integer]+numbers5Smooth =+   foldr+      (\(x0:x1:x2:xs) ys -> x0 : x1 : x2 : ListHT.mergeBy (<=) xs ys)+      (error "numbers5Smooth: infinite list should not have an end") $+   iterate (map (5*)) $+   numbers3Smooth++numbers5SmoothAlt :: [Integer]+numbers5SmoothAlt =+   unfoldr+      (fmap (\(m,rest) -> (m, Set.union rest $ Set.fromAscList [2*m,3*m,5*m])) .+       Set.minView) $+   Set.singleton 1++ceilingPowerOfTwo :: (Ring.C a, Bits a) => a -> a+ceilingPowerOfTwo 0 = 1+ceilingPowerOfTwo n =+   (1+) $ fst $ head $+   dropWhile (uncurry (/=)) $+   ListHT.mapAdjacent (,) $+   scanl (\m d -> shiftR m d .|. m) (n-1) $+   iterate (2*) 1++divideByMaximumPower ::+   (Integral.C a, ZeroTestable.C a) => a -> a -> a+divideByMaximumPower b n =+   last $+   n : unfoldr (\m -> case divMod m b of (q,r) -> toMaybe (isZero r) (q,q)) n++divideByMaximumPowerRecursive ::+   (Integral.C a, Eq a, ZeroTestable.C a) => a -> a -> a+divideByMaximumPowerRecursive b =+   let recourse n =+          case divMod b n of+             (q,0) -> recourse q+             _ -> n+   in  recourse++getMaximumExponent ::+   (Integral.C a, ZeroTestable.C a) =>+   a -> a -> (Int,a)+getMaximumExponent b n =+   last $ (0,n) :+   unfoldr+      (\(e,m) ->+         let (q,r) = divMod m b+             eq = (e+1,q)+         in  toMaybe (isZero r) (eq,eq))+      (0,n)++is3Smooth :: Integer -> Bool+is3Smooth =+   (1==) .+   divideByMaximumPower 3 .+   divideByMaximumPower 2++is5Smooth :: Integer -> Bool+is5Smooth =+   (1==) .+   divideByMaximumPower 5 .+   divideByMaximumPower 3 .+   divideByMaximumPower 2++{- |+Compute the smallest composite of 2 and 3 that is as least as large as the input.+This can be interpreted as solving an integer linear programming problem with+min (\(a,b) -> a * log 2 + b * log 3)+over the domain {(a,b) : a>=0, b>=0, a * log 2 + b * log 3 >= log n}+-}+{-+Problem: We cannot just start with the ceilingPowerOfTwo+and then multiply with 3/4 until we fall below n,+since the 3/4 decreases too fast.+27/32 is closer to one,+and higher powers of 3 and 2 in the ratio make the ratio even closer to one.+-}+ceiling3Smooth :: Integer -> Integer+ceiling3Smooth n =+   head $ dropWhile (<n) numbers3Smooth++ceiling5Smooth :: Integer -> Integer+ceiling5Smooth n =+   head $ dropWhile (<n) numbers5Smooth++ceiling3SmoothNaive :: Integer -> Integer+ceiling3SmoothNaive =+   head .+   dropWhile (not . is3Smooth) .+   iterate (1+)++ceiling5SmoothNaive :: Integer -> Integer+ceiling5SmoothNaive =+   head .+   dropWhile (not . is5Smooth) .+   iterate (1+)+++{- |+Compute all primes that occur in the course of dividing+a Fourier transform into sub-transforms.+-}+partialPrimes :: Integer -> [Integer]+partialPrimes =+   let primeFactorSet =+          Set.fromAscList . uniquePrimeFactors+   in  unfoldr+         (fmap+             (\(p,set) ->+                (p, Set.union (primeFactorSet (p-1)) set)) .+          Set.maxView)+       .+       primeFactorSet++-- cf. htam:NumberTheory+uniquePrimeFactors ::+   (Integral.C a, Bits a, ZeroTestable.C a, Ord a) =>+   a -> [a]+--   map snd . primeFactors+uniquePrimeFactors n =+   let oddFactors =+          foldr+             (\p go m ->+                let (q,r) = divMod m p+                in  if r==0+                      then p : go (divideByMaximumPower p q)+                      else+                        if q >= p+                          then go m+                          else if m==1 then [] else m : [])+             (error "uniquePrimeFactors: end of infinite list")+             (iterate (2+) 3)+   in  case powerOfTwoFactors n of+          (1,m) -> oddFactors m+          (_,m) -> 2 : oddFactors m++{- |+Prime factors and their exponents in ascending order.+-}+primeFactors ::+   (PrimitiveRoot a, Ord a) => a -> [(Int, a)]+primeFactors n =+   let oddFactors =+          foldr+             (\p go m ->+                let (q0,r) = divMod m p+                in  if r==0+                      then+                        case getMaximumExponent p q0 of+                          (e,q1) -> (e+1,p) : go q1+                      else+                        if q0 >= p+                          then go m+                          else if m==1 then [] else (1,m) : [])+             (const [])+             (filter (not . Units.isUnit) $+              primitiveRootCandidates n)+   in  case getMaximumExponent 2 n of+          (0,m) -> oddFactors m+          (e,m) -> (e,2) : oddFactors m++{-+cf. htam:NumberTheory++Shall this be moved to NumericPrelude?++It should be replaced by a more sophisticated prime test.+-}+isPrime :: Integer -> Bool+isPrime n =+   case primeFactors n of+      [] -> False+      (e,m):_ -> e==1 && m==n++{- |+Find lengths of signals that require many interim Rader transforms+and end with the given length.++raderWorstCases 2  =  OEIS-A061092+raderWorstCases 5  =  tail OEIS-A059411++Smallest raderWorstCase numbers are 2,5,13,17,19,31,37,41,43,61,...+This matches the definition of OEIS-A061303.+-}+raderWorstCases :: Integer -> [Integer]+raderWorstCases =+   iterate+      (\n ->+         head $ dropWhile (not . isPrime) $+         tail $ iterate (n+) 1)++{- |+This is usually faster than 'fastFourierRing'+since it does not need to factor large numbers.+However, the generated modulus is usually much greater.+-}+{-+I see the following opportunities for optimization:++1. Speedup 'fastFourierRing' by+   faster primality test (Miller-Rabin) and+   faster prime factorization (Pollard-Rho-method).+   These are also needed for+   maximumOrderOfPrimitiveRootsOfUnityInteger+   that is used by Fourier.Element.primitiveRoot+   in order to compute a root with maximum order.++2. Reduce the moduli produced by 'fastFourierRingAlt'+   by merging some orders that are passed to+   ringWithPrimitiveRootsOfUnityAndUnits,+   such that an LCM of a group of orders can still be handled.+   This is a kind of knapsack problem.+   Maybe we could collect the factors in a way+   such that (lcm orderGroup + 1) is prime.++3. Avoid to compute factorizations of numbers+   where we already know the factors+   or where we do not need the factors at all.+   Use the factors returned by partialPrimes+   in order to compute a prime factorization+   of lcmMulti (map pred (partialPrimes n)).+   Call this (product ps).+   Now search for rings of moduli (1 + k * product ps),+   where there are (small) primitive roots of order (product ps).+   We only need to check whether there are small numbers+   such as 2, 3, 5, 6, 7 that have a (product ps)-th power that is 1,+   using fast exponentiation.+   If there is a power being 1 then check for primitivity+   by computing (k * product ps / p)-th powers+   for all prime factors p of (k * product ps).+   If there is no primitive root <= 7,+   there may still be a primitive root of wanted order,+   but it is then cheaper to seek for larger moduli.++   If we finally have a nice modulus+   it is still stupid to factorize (modulus-1)+   and search for a primitive root+   in each invocation of Fourier.Element.primitiveRoot.+   We could define a special datatype analogously to ResidueClass,+   that stores the primitive root and its order+   additional to the ResidueClass modulus.+-}+fastFourierRingAlt :: Int -> Integer+fastFourierRingAlt n =+   case n of+      0 -> 2+      1 -> 2+      _ ->+         let ni = fromIntegral n+             ps = filter (>1) (map (subtract 1) (partialPrimes ni))+         in  ringWithPrimitiveRootsOfUnityAndUnits (map Order $ ni : ps) ps++{- |+Determine an integer residue ring+in which a Fast Fourier transform of size n can be performed.+It must contain certain primitive roots.+If we choose a non-primitive root,+then some off-diagonal values in F^-1·F are non-zero.+-}+{-+When we need roots of orders o1,...,ok and according inverse elements+we can also ask for a ring, where there is a root of order lcm(o1,...,ok).+The answer to both questions is the same set of rings.+This can be proven using the statement,+that the order of any primitive root+divides the carmichael value of the modulus.++Since ringWithPrimitiveRootsOfUnityAndUnits+multiplies the moduli of rings for o1,...,ok,+it will produce large moduli.+-}+fastFourierRing :: Int -> Integer+fastFourierRing n =+   case n of+      0 -> 2+      1 -> 2+      _ ->+         let ni = fromIntegral n+         in  {-+             We cannot use ringsWithPrimitiveRootOfUnityAndUnit+             since for 359 we already get an Int overflow.+             For 719, 1439, 2879 we also get a very large value.+             -}+             head $ filter isPrime $+             (\order -> iterate (order +) 1) $+             lcmMulti $+             ni : map (subtract 1) (partialPrimes ni)
src/Synthesizer/Basic/ToneModulation.hs view
@@ -5,16 +5,11 @@  import Synthesizer.Interpolation (Margin, marginOffset, marginNumber, ) --- import qualified Algebra.Transcendental        as Trans import qualified Algebra.RealField             as RealField import qualified Algebra.Field                 as Field--- import qualified Algebra.RealRing                  as RealRing import qualified Algebra.Ring                  as Ring-import qualified Algebra.Additive              as Additive  import NumericPrelude.Numeric---- import qualified Prelude as P import NumericPrelude.Base
+ src/Synthesizer/Causal/Analysis.hs view
@@ -0,0 +1,34 @@+{-# LANGUAGE NoImplicitPrelude #-}+module Synthesizer.Causal.Analysis where++import qualified Synthesizer.Causal.Filter.Recursive.Integration as Integration++import qualified Synthesizer.Causal.Process as Causal+import qualified Synthesizer.Plain.Analysis as Ana++import qualified Algebra.RealRing              as RealRing++import Control.Arrow (second, (^<<), (<<^), )++-- import qualified Prelude as P+import NumericPrelude.Base+import NumericPrelude.Numeric+++deltaSigmaModulation ::+   RealRing.C y => Causal.T y Ana.BinaryLevel+deltaSigmaModulation =+   Causal.feedback+      ((Ana.binaryLevelFromBool . (zero <=)) ^<<+       Integration.run <<^+       uncurry (-))+      (Causal.consInit zero <<^ Ana.binaryLevelToNumber)++deltaSigmaModulationPositive ::+   RealRing.C y => Causal.T (y, y) y+deltaSigmaModulationPositive =+   Causal.feedback+      ((\(threshold,xi) -> if threshold<=xi then threshold else zero) ^<<+       second Integration.run <<^+       (\((threshold,xi),cum) -> (threshold,xi-cum)))+      (Causal.consInit zero)
src/Synthesizer/Causal/Displacement.hs view
@@ -7,7 +7,6 @@ import qualified Algebra.Ring                  as Ring import qualified Algebra.Additive              as Additive -import qualified Prelude as P import NumericPrelude.Base import NumericPrelude.Numeric
src/Synthesizer/Causal/Interpolation.hs view
@@ -24,8 +24,6 @@ import qualified Algebra.RealRing  as RealRing import qualified Algebra.Additive  as Additive -import Algebra.Additive(zero)-  import NumericPrelude.Base import NumericPrelude.Numeric
src/Synthesizer/Causal/Oscillator.hs view
@@ -29,13 +29,10 @@ import qualified Algebra.Transcendental        as Trans import qualified Algebra.RealField             as RealField import qualified Algebra.RealRing              as RealRing-import qualified Algebra.Ring                  as Ring -import Control.Arrow ((^<<), (<<^), (<<<), (***), second, )+import Control.Arrow ((^<<), (<<^), (<<<), (***), )  import NumericPrelude.Numeric--import qualified Prelude as P import NumericPrelude.Base
src/Synthesizer/Causal/Oscillator/Core.hs view
@@ -21,13 +21,10 @@ import qualified Synthesizer.State.Signal as Sig  import qualified Algebra.RealRing             as RealRing-import qualified Algebra.Additive              as Additive  import Control.Arrow ((^<<), (&&&), second, returnA, )  import NumericPrelude.Numeric--import qualified Prelude as P import NumericPrelude.Base
+ src/Synthesizer/Causal/Spatial.hs view
@@ -0,0 +1,24 @@+{-# LANGUAGE NoImplicitPrelude #-}+module Synthesizer.Causal.Spatial where++import qualified Algebra.NormedSpace.Euclidean as Euc+import qualified Algebra.Field                 as Field++import Control.Arrow (Arrow, arr, )++import NumericPrelude.Base+import NumericPrelude.Numeric+++{-|+simulate an moving sounding object++convert the way of the object through 3D space+into a delay and attenuation information,+sonicDelay is the reciprocal of the sonic velocity+-}+receive3Dsound ::+   (Field.C a, Euc.C a v, Arrow arrow) =>+   a -> a -> v -> arrow v (a,a)+receive3Dsound att sonicDelay ear =+   arr ((\dist -> (sonicDelay*dist, 1/(att+dist)^2)) . Euc.norm . subtract ear)
src/Synthesizer/ChunkySize/Cut.hs view
@@ -9,45 +9,24 @@ import qualified Synthesizer.Generic.Cut as Cut import qualified Synthesizer.Generic.Signal as SigG --- import qualified Synthesizer.Plain.Signal as Sig import qualified Synthesizer.State.Signal as SigS--- import qualified Synthesizer.Storable.Signal as SigSt import qualified Data.StorableVector.Lazy.Pattern as SigStV import qualified Data.StorableVector.Lazy as Vector  -import qualified Algebra.Ring as Ring--- import qualified Algebra.ToInteger as ToInteger---- import qualified Number.NonNegative as NonNegW--- import qualified Algebra.NonNegative as NonNeg import qualified Number.NonNegativeChunky as Chunky -{---- import qualified Numeric.NonNegative.Wrapper as NonNegW98-import qualified Numeric.NonNegative.Class as NonNeg98-import qualified Numeric.NonNegative.Chunky as Chunky98--}- import Foreign.Storable (Storable)  import qualified Data.List as List import qualified Data.List.Match as Match import Data.Tuple.HT (mapPair, ) -import qualified Data.Monoid as Monoid import Data.Monoid (Monoid, ) -import qualified Prelude as P+import Prelude () import NumericPrelude.Numeric import NumericPrelude.Base hiding (splitAt, Read, )-{--import Prelude-   (Bool, Int, String, (++), error, const,-    pred, (<=), (>=), (<), (>), ($),-    (.), not, (||), (&&),-    Maybe(Just, Nothing), )--}   class Cut.Read sig => Read sig where
src/Synthesizer/ChunkySize/Signal.hs view
@@ -14,10 +14,6 @@ import qualified Data.StorableVector.Lazy.Pattern as SigStV import qualified Data.StorableVector.Lazy as Vector --- import qualified Algebra.NonNegative as NonNeg--- import qualified Algebra.Module   as Module--- import qualified Algebra.Additive as Additive- import Foreign.Storable (Storable)  import qualified Data.List.Match as Match@@ -28,7 +24,7 @@  -- import NumericPrelude.Numeric import Prelude-   (Bool, Int, Maybe(Just), fst, (.), id, )+   (Maybe(Just), fst, (.), id, )   class (SigG.Write sig y, Cut.Transform (sig y)) => Write sig y where
src/Synthesizer/Generic/Analysis.hs view
@@ -8,31 +8,15 @@ import qualified Synthesizer.Generic.Signal as SigG import qualified Synthesizer.Generic.Signal2 as SigG2 --- import qualified Synthesizer.Plain.Control as Ctrl---- import qualified Algebra.Module                as Module--- import qualified Algebra.Transcendental        as Trans import qualified Algebra.Algebraic             as Algebraic--- import qualified Algebra.RealField             as RealField import qualified Algebra.Field                 as Field import qualified Algebra.RealRing              as RealRing-import qualified Algebra.Absolute              as Absolute import qualified Algebra.Ring                  as Ring-import qualified Algebra.Additive              as Additive  import qualified Algebra.NormedSpace.Maximum   as NormedMax import qualified Algebra.NormedSpace.Euclidean as NormedEuc import qualified Algebra.NormedSpace.Sum       as NormedSum --- import qualified Data.Array as Array---- import qualified Data.IntMap as IntMap---- import Algebra.Module((*>))---- import Data.Array (accumArray)--import qualified Prelude as P import NumericPrelude.Base import NumericPrelude.Numeric @@ -310,18 +294,16 @@             then not(x<lower)             else x>upper) -{- {- | Almost naive implementation of the chirp transform, a generalization of the Fourier transform.  More sophisticated algorithms like Rader, Cooley-Tukey, Winograd, Prime-Factor may follow. -}-chirpTransform :: Ring.C y =>-   y -> sig y -> sig y-chirpTransform z xs =-   let powers = Ctrl.curveMultiscaleNeutral (*) z one-       powerPowers =-          SigG.map (\zn -> Ctrl.curveMultiscaleNeutral (*) zn one) powers-   in  SigG.map (scalarProduct xs) powerPowers--}+chirpTransform ::+   (SigG.Write sig y, Ring.C y) =>+   SigG.LazySize -> y -> sig y -> sig y+chirpTransform size z =+   SigG.fromState size .+   Ana.chirpTransform z .+   SigG.toState
src/Synthesizer/Generic/Control.hs view
@@ -12,17 +12,11 @@ import qualified Algebra.Transcendental        as Trans import qualified Algebra.RealRing              as RealRing import qualified Algebra.Field                 as Field-import qualified Algebra.Ring                  as Ring import qualified Algebra.Additive              as Additive -import Algebra.Module((*>))- import Number.Complex (cis,real) import qualified Number.Complex as Complex --- import Control.Applicative ((<$>), )--import qualified Prelude as P import NumericPrelude.Base import NumericPrelude.Numeric
src/Synthesizer/Generic/Cut.hs view
@@ -11,21 +11,18 @@ import qualified Synthesizer.State.Signal as SigS -- import qualified Synthesizer.Storable.Signal as SigSt import qualified Data.StorableVector as SV-import qualified Data.StorableVector.Lazy as Vector+import qualified Data.StorableVector.Lazy as SVL  import qualified Algebra.ToInteger as ToInteger import qualified Algebra.Ring as Ring  import qualified Data.EventList.Relative.BodyTime as EventList --- import qualified Number.NonNegative as NonNegW import qualified Algebra.NonNegative as NonNeg import qualified Number.NonNegativeChunky as Chunky --- import qualified Numeric.NonNegative.Wrapper as NonNegW98 import qualified Numeric.NonNegative.Class as NonNeg98 import qualified Numeric.NonNegative.Chunky as Chunky98-import Numeric.NonNegative.Class ((-|), )  import Foreign.Storable (Storable, ) import Control.DeepSeq (NFData, rnf, )@@ -41,10 +38,10 @@ import qualified Prelude as P import NumericPrelude.Numeric import Prelude-   (Bool, Int, String, (++), error,+   (Bool, String, (++), error,     pred, (<=), (>=), (<),     (.), ($), const, snd,-    not, (||), (&&), min, )+    not, (||), (&&), min, max, )   class Read sig where@@ -81,48 +78,74 @@    reverse :: sig -> sig  +instance Storable y => Read (SV.Vector y) where+   {-# INLINE null #-}+   null = SV.null+   {-# INLINE length #-}+   length = SV.length++instance (Storable y) => NormalForm (SV.Vector y) where+   {-# INLINE evaluateHead #-}+   evaluateHead x =+      if SV.null x then () else ()++instance Storable y => Transform (SV.Vector y) where+   {-# INLINE take #-}+   take = SV.take+   {-# INLINE drop #-}+   drop = SV.drop+   {-# INLINE splitAt #-}+   splitAt = SV.splitAt+   {-# INLINE dropMarginRem #-}+   dropMarginRem n m xs =+      let d = min m $ max 0 $ SV.length xs - n+      in  (m-d, SV.drop d xs)+   {-# INLINE reverse #-}+   reverse = SV.reverse++ -- instance Storable y => Read SigSt.T y where-instance Storable y => Read (Vector.Vector y) where+instance Storable y => Read (SVL.Vector y) where    {-# INLINE null #-}-   null = Vector.null+   null = SVL.null    {-# INLINE length #-}-   length = Vector.length+   length = SVL.length -instance (Storable y) => NormalForm (Vector.Vector y) where+instance (Storable y) => NormalForm (SVL.Vector y) where    {-# INLINE evaluateHead #-}    evaluateHead =-      ListHT.switchL () (\x _ -> if SV.null x then () else ()) . Vector.chunks---      ListHT.switchL () (\x _ -> rnf x) . Vector.chunks+      ListHT.switchL () (\x _ -> evaluateHead x) . SVL.chunks+--      ListHT.switchL () (\x _ -> rnf x) . SVL.chunks --   evaluateHead x =---      if Vector.null x then () else ()+--      if SVL.null x then () else ()  {--instance (Storable y, NFData y) => NormalForm (Vector.Vector y) where+instance (Storable y, NFData y) => NormalForm (SVL.Vector y) where    {-# INLINE evaluateHead #-}-   evaluateHead x = Vector.switchL () (\x _ -> rnf x)+   evaluateHead x = SVL.switchL () (\x _ -> rnf x) -} -instance Storable y => Transform (Vector.Vector y) where+instance Storable y => Transform (SVL.Vector y) where    {-    {-# INLINE empty #-}-   empty = Vector.empty+   empty = SVL.empty    {-# INLINE cycle #-}-   cycle = Vector.cycle+   cycle = SVL.cycle    {-# INLINE append #-}-   append = Vector.append+   append = SVL.append    {-# INLINE concat #-}-   concat = Vector.concat+   concat = SVL.concat    -}    {-# INLINE take #-}-   take = Vector.take+   take = SVL.take    {-# INLINE drop #-}-   drop = Vector.drop+   drop = SVL.drop    {-# INLINE splitAt #-}-   splitAt = Vector.splitAt+   splitAt = SVL.splitAt    {-# INLINE dropMarginRem #-}-   dropMarginRem = Vector.dropMarginRem+   dropMarginRem = SVL.dropMarginRem    {-# INLINE reverse #-}-   reverse = Vector.reverse+   reverse = SVL.reverse   instance Read ([] y) where@@ -301,8 +324,8 @@    dropMarginRem n m x =       let (z,~(b,d)) =              Chunky.minMaxDiff-                (intToChunky "dropMargin/n" n)-                (x NonNeg.-| intToChunky "dropMargin/m" m)+                (intToChunky "dropMargin/n" m)+                (x NonNeg.-| intToChunky "dropMargin/m" n)       in  (if b then 0 else fromIntegral (Chunky.toNumber d),            x NonNeg.-| z)    {-# INLINE splitAt #-}@@ -346,8 +369,8 @@    dropMarginRem n m x =       let (z,~(b,d)) =              NonNeg98.split-                (intToChunky98 "dropMargin/n" n)-                (x NonNeg98.-| intToChunky98 "dropMargin/m" m)+                (intToChunky98 "dropMargin/n" m)+                (x NonNeg98.-| intToChunky98 "dropMargin/m" n)       in  (if b then 0 else P.fromIntegral (Chunky98.toNumber d),            x NonNeg98.-| z)    {-# INLINE splitAt #-}
+ src/Synthesizer/Generic/Cyclic.hs view
@@ -0,0 +1,192 @@+module Synthesizer.Generic.Cyclic where++import qualified Synthesizer.Generic.Filter.NonRecursive as FiltNRG+import qualified Synthesizer.Generic.Analysis as AnaG+import qualified Synthesizer.Generic.Signal2 as SigG2+import qualified Synthesizer.Generic.Signal as SigG+import qualified Synthesizer.Generic.Cut as CutG+import qualified Synthesizer.State.Signal as Sig++import qualified Algebra.Ring as Ring+import qualified Algebra.Additive as Additive++import NumericPrelude.Numeric+import NumericPrelude.Base+import Prelude ()+++fromSignal ::+   (SigG.Write sig yv, Additive.C yv) =>+   SigG.LazySize -> Int -> sig yv -> sig yv+fromSignal chunkSize n =+   {- almost Sig.sum -}+   Sig.foldL SigG.mix (SigG.replicate chunkSize n zero) .+   CutG.sliceVertical n++reverse ::+   (SigG.Transform sig y) =>+   sig y -> sig y+reverse sig =+   SigG.switchL sig+      (\y ys -> SigG.cons y (SigG.reverse ys)) sig+++{- |+It must hold @n <= CutG.length x@.+-}+reperiodize ::+   (SigG.Transform sig yv, Additive.C yv) =>+   Int -> sig yv -> sig yv+reperiodize n =+   {- Sig.sum -}+   Sig.foldL SigG.mix CutG.empty .+   CutG.sliceVertical n++{- |+length of the input signals must be equal+-}+convolve ::+   (SigG2.Transform sig y y, Ring.C y) =>+   sig y -> sig y -> sig y+convolve x y =+   reperiodize (CutG.length x) $+   FiltNRG.karatsubaFinite (*) x y++++{- |+The size of both input signals must be equal.++Could be optimized by computing only first (length x) elements.+-}+filterNaive ::+   (SigG.Transform sig y, Ring.C y) =>+   sig y -> sig y -> sig y+filterNaive x y =+   SigG.takeStateMatch y $+   SigG.toState $+   SigG.mapTails+      (AnaG.scalarProduct x)+      (SigG.append y y)++convolveNaive ::+   (SigG.Transform sig y, Ring.C y) =>+   sig y -> sig y -> sig y+convolveNaive x y =+   SigG.takeStateMatch y $+   SigG.toState $+   SigG.mapTails+      (AnaG.scalarProduct (SigG.reverse x))+      (SigG.laxTail $ SigG.append y y)+++{-+Some small size convolutions using the Karatsuba trick.+We do not use Toom-3 multiplication,+because this requires division by 2 and 6.++In principle we could implement them+by calling the corresponding functions in Filter.NonRecursive+and periodize them afterwards.+However the custom implementations below+allow a litte bit more optimization,+namely sharing of some sums.+-}++type Pair y = (y,y)++{-# INLINE convolvePair #-}+convolvePair ::+   (Ring.C y) =>+   Pair y -> Pair y -> Pair y+convolvePair a b =+   snd $ sumAndConvolvePair a b++{-# INLINE sumAndConvolvePair #-}+sumAndConvolvePair ::+   (Ring.C y) =>+   Pair y -> Pair y -> ((y,y), Pair y)+sumAndConvolvePair (a0,a1) (b0,b1) =+   let sa01 = a0+a1+       sb01 = b0+b1+       ab0ab1 = a0*b0+a1*b1+   in  ((sa01, sb01), (ab0ab1, sa01*sb01-ab0ab1))+++type Triple y = (y,y,y)++{-# INLINE convolveTriple #-}+convolveTriple ::+   (Ring.C y) =>+   Triple y -> Triple y -> Triple y+convolveTriple a b =+   snd $ sumAndConvolveTriple a b++{-# INLINE sumAndConvolveTriple #-}+sumAndConvolveTriple ::+   (Ring.C y) =>+   Triple y -> Triple y -> ((y,y), Triple y)+sumAndConvolveTriple (a0,a1,a2) (b0,b1,b2) =+   let ab0 = a0*b0+       dab12 = a1*b1 - a2*b2+       sa01 = a0+a1; sb01 = b0+b1; tab01 = sa01*sb01 - ab0+       sa02 = a0+a2; sb02 = b0+b2; tab02 = sa02*sb02 - ab0+       sa012 = sa01+a2+       sb012 = sb01+b2++       d0 = sa012*sb012 - tab01 - tab02+       d1 = tab01 - dab12+       d2 = tab02 + dab12+   in  ((sa012, sb012), (d0, d1, d2))++{-# INLINE sumAndConvolveTripleAlt #-}+sumAndConvolveTripleAlt ::+   (Ring.C y) =>+   Triple y -> Triple y -> ((y,y), Triple y)+sumAndConvolveTripleAlt (a0,a1,a2) (b0,b1,b2) =+   let ab0 = a0*b0+       ab1 = a1*b1+       ab2 = a2*b2+       sa01 = a0+a1; sb01 = b0+b1+       ab01 = sa01*sb01 - (ab0+ab1)+       sa02 = a0+a2; sb02 = b0+b2+       ab02 = sa02*sb02 - (ab0+ab2)+       sa12 = a1+a2; sb12 = b1+b2+       ab12 = sa12*sb12 - (ab1+ab2)+   in  ((sa01+a2, sb01+b2), (ab0+ab12, ab2+ab01, ab1+ab02))+++type Quadruple y = (y,y,y,y)++{-# INLINE convolveQuadruple #-}+convolveQuadruple ::+   (Ring.C y) =>+   Quadruple y -> Quadruple y -> Quadruple y+convolveQuadruple a b =+   snd $ sumAndConvolveQuadruple a b++{-# INLINE sumAndConvolveQuadruple #-}+sumAndConvolveQuadruple ::+   (Ring.C y) =>+   Quadruple y -> Quadruple y -> ((y,y), Quadruple y)+sumAndConvolveQuadruple (a0,a1,a2,a3) (b0,b1,b2,b3) =+   let ab0 = a0*b0+       ab1 = a1*b1+       sa01 = a0+a1; sb01 = b0+b1+       ab01 = sa01*sb01 - (ab0+ab1)+       ab2 = a2*b2+       ab3 = a3*b3+       sa23 = a2+a3; sb23 = b2+b3+       ab23 = sa23*sb23 - (ab2+ab3)+       c0 = ab0  + ab2 - (ab1 + ab3)+       c1 = ab01 + ab23+       ab02 = (a0+a2)*(b0+b2)+       ab13 = (a1+a3)*(b1+b3)+       sa0123 = sa01+sa23+       sb0123 = sb01+sb23+       ab0123 = sa0123*sb0123 - (ab02+ab13)+       d0 = ab13   + c0+       d1 =          c1+       d2 = ab02   - c0+       d3 = ab0123 - c1+   in  ((sa0123, sb0123), (d0, d1, d2, d3))
src/Synthesizer/Generic/Filter/NonRecursive.hs view
@@ -1,7 +1,7 @@ {-# LANGUAGE NoImplicitPrelude #-} {-# LANGUAGE FlexibleContexts #-} {- |-Copyright   :  (c) Henning Thielemann 2008-2009+Copyright   :  (c) Henning Thielemann 2008-2011 License     :  GPL  Maintainer  :  synthesizer@henning-thielemann.de@@ -12,8 +12,10 @@  import qualified Synthesizer.Generic.Signal as SigG import qualified Synthesizer.Generic.Signal2 as SigG2-+import qualified Synthesizer.Generic.Cut as CutG import qualified Synthesizer.Generic.Control as Ctrl+import qualified Synthesizer.Generic.LengthSignal as SigL+ import qualified Synthesizer.State.Signal as SigS import qualified Synthesizer.Plain.Filter.NonRecursive as Filt import qualified Synthesizer.State.Filter.NonRecursive as FiltS@@ -25,8 +27,6 @@ import qualified Algebra.Ring           as Ring import qualified Algebra.Additive       as Additive -import Algebra.Module( {- linearComb, -} (*>), )- import Control.Monad (mplus, ) import Data.Function.HT (nest, ) import Data.Tuple.HT (mapSnd, mapPair, )@@ -36,7 +36,7 @@ import NumericPrelude.Numeric as NP  -{- * Envelope application -}+-- * Envelope application  {-# INLINE negate #-} negate ::@@ -99,7 +99,7 @@        envelope leadOut partOut  -{- * Smoothing -}+-- * Delay  {-# INLINE delay #-} delay :: (Additive.C y, SigG.Write sig y) =>@@ -153,6 +153,7 @@    SigG.append (SigG.replicate size n zero)  +-- * smoothing  binomialMask ::    (Field.C a, SigG.Write sig a) =>@@ -169,16 +170,6 @@ property: must sum up to 1 -} -{-| Unmodulated non-recursive filter -}-{-# INLINE generic #-}-generic ::-   (Module.C a v, SigG.Transform sig a, SigG.Write sig v) =>-   sig a -> sig v -> sig v-generic m x =-   let mr = SigG.reverse m-       xp = delayPos (pred (SigG.length m)) x-   in  SigG.mapTails (SigG.linearComb mr) xp- {- {- | @eps@ is the threshold relatively to the maximum.@@ -588,3 +579,335 @@    (Additive.C v, SigG.Transform sig v) =>    sig v -> sig v differentiate2 = differentiate . differentiate+++-- * general non-recursive filters++{-|+Unmodulated non-recursive filter (convolution)++Brute force implementation.+-}+{-# INLINE generic #-}+generic ::+   (Module.C a v, SigG.Transform sig a, SigG.Write sig v) =>+   sig a -> sig v -> sig v+generic m x =+   if SigG.null m || SigG.null x+     then CutG.empty+     else+       let mr = SigG.reverse m+           xp = delayPos (pred (SigG.length m)) x+       in  SigG.mapTails (SigG.linearComb mr) xp+++{- |+Both should signals should have similar length.+If they have considerably different length,+then better use 'karatsubaFiniteInfinite'.++Implementation using Karatsuba trick and split-and-overlap-add.+This way we stay in a ring, are faster than quadratic runtime+but do not reach log-linear runtime.+-}+karatsubaFinite ::+   (Additive.C a, Additive.C b, Additive.C c,+    SigG2.Transform sig a c, SigG2.Transform sig b c) =>+   (a -> b -> c) ->+   sig a -> sig b -> sig c+karatsubaFinite mul a b =+   SigL.toSignal $+   karatsubaBounded mul+      (SigL.fromSignal a) (SigL.fromSignal b)++{-# INLINE karatsubaBounded #-}+karatsubaBounded ::+   (Additive.C a, Additive.C b, Additive.C c,+    SigG2.Transform sig a c, SigG2.Transform sig b c) =>+   (a -> b -> c) ->+   SigL.T (sig a) -> SigL.T (sig b) -> SigL.T (sig c)+karatsubaBounded mul a b =+   case (SigL.length a, SigL.length b) of+      (0,_) -> CutG.empty+      (_,0) -> CutG.empty+      (1,_) ->+         SigG.switchL+            (error "karatsubaBounded: empty signal")+            (\y _ -> fmap (SigG2.map (mul y)) b) $+         SigL.body a+      (_,1) ->+         SigG.switchL+            (error "karatsubaBounded: empty signal")+            (\y _ -> fmap (SigG2.map (flip mul y)) a) $+         SigL.body b+      (2,2) ->+         let [a0,a1] = SigG.toList (SigL.toSignal a)+             [b0,b1] = SigG.toList (SigL.toSignal b)+             (c0,c1,c2) = convolvePair mul (a0,a1) (b0,b1)+         in  SigL.Cons 3 $ rechunk a b $+             c0 : c1 : c2 : []+      (2,3) ->+         let [a0,a1]    = SigG.toList (SigL.toSignal a)+             [b0,b1,b2] = SigG.toList (SigL.toSignal b)+             (c0,c1,c2,c3) =+                convolvePairTriple mul (a0,a1) (b0,b1,b2)+         in  SigL.Cons 4 $ rechunk a b $+             c0 : c1 : c2 : c3 : []+      (3,2) ->+         let [a0,a1,a2] = SigG.toList (SigL.toSignal a)+             [b0,b1]    = SigG.toList (SigL.toSignal b)+             (c0,c1,c2,c3) =+                convolvePairTriple (flip mul) (b0,b1) (a0,a1,a2)+         in  SigL.Cons 4 $ rechunk a b $+             c0 : c1 : c2 : c3 : []+      (3,3) ->+         let [a0,a1,a2] = SigG.toList (SigL.toSignal a)+             [b0,b1,b2] = SigG.toList (SigL.toSignal b)+             (c0,c1,c2,c3,c4) =+                convolveTriple mul (a0,a1,a2) (b0,b1,b2)+         in  SigL.Cons 5 $ rechunk a b $+             c0 : c1 : c2 : c3 : c4 : []+      (4,4) ->+         let [a0,a1,a2,a3] = SigG.toList (SigL.toSignal a)+             [b0,b1,b2,b3] = SigG.toList (SigL.toSignal b)+             (c0,c1,c2,c3,c4,c5,c6) =+                convolveQuadruple mul (a0,a1,a2,a3) (b0,b1,b2,b3)+         in  SigL.Cons 7 $ rechunk a b $+             c0 : c1 : c2 : c3 : c4 : c5 : c6 : []+      (lenA,lenB) ->+         let n2 = div (max lenA lenB) 2+             (a0,a1) = SigL.splitAt n2 a+             (b0,b1) = SigL.splitAt n2 b+             (c0,c1,c2) =+                convolvePair+                   (karatsubaBounded mul)+                   (a0,a1) (b0,b1)+         in  fmap (rechunk a b) $+             SigL.addShiftedSimple n2 c0 $+             SigL.addShiftedSimple n2 c1 c2++{-# INLINE rechunk #-}+rechunk ::+   (SigG2.Transform sig1 a c, SigG2.Transform sig1 b c,+    SigG.Transform sig0 c) =>+   SigL.T (sig1 a) -> SigL.T (sig1 b) -> sig0 c -> sig1 c+rechunk a b c =+   let (ac,bc) = CutG.splitAt (SigL.length a) c+   in  SigG2.takeStateMatch (SigL.body a) (SigG.toState ac)+       `SigG.append`+       SigG2.takeStateMatch (SigL.body b) (SigG.toState bc)+++{- |+The first operand must be finite and+the second one can be infinite.+For efficient operation we expect that the second signal+is longer than the first one.+-}+{-+Implemented by overlap-add of pieces that are convolved by Karatsuba trick.+Is it more efficient to round the chunk size up to the next power of two?+Can we make use of the fact,+that the first operand is always split in the same way?+-}+karatsubaFiniteInfinite ::+   (Additive.C a, Additive.C b, Additive.C c,+    SigG2.Transform sig a c, SigG2.Transform sig b c) =>+   (a -> b -> c) ->+   sig a -> sig b -> sig c+karatsubaFiniteInfinite mul a b =+   let al = SigL.fromSignal a+   in  case SigL.length al of+          0 -> CutG.empty+          alen ->+             SigS.foldR (addShiftedSimple alen) CutG.empty $+             SigS.map SigL.toSignal $+             SigS.map (karatsubaBounded mul al . SigL.fromSignal) $+             SigG.sliceVertical alen b+++karatsubaInfinite ::+   (Additive.C a, Additive.C b, Additive.C c,+    SigG2.Transform sig a c, SigG2.Transform sig b c) =>+   (a -> b -> c) ->+   sig a -> sig b -> sig c+karatsubaInfinite mul =+   let recourse n a b =+          let (a0,a1) = SigG.splitAt n a+              (b0,b1) = SigG.splitAt n b+              {-+              We could also apply Karatsuba's trick to these pairs.+              But this requires Additive (sig a) constraint+              and I do not know whether this is actually an optimization.+              -}+              ab00 =+                 SigL.toSignal $+                 karatsubaBounded mul+                    (SigL.fromSignal a0) (SigL.fromSignal b0)+              ab01 = karatsubaFiniteInfinite mul a0 b1+              ab10 = karatsubaFiniteInfinite (flip mul) b0 a1+              ab11 = recourse (2*n) a1 b1+          in  if SigG.null a || SigG.null b+                then CutG.empty+                else+                  addShiftedSimple n ab00 $+                  addShiftedSimple n (SigG.mix ab01 ab10) ab11+   in  recourse 1+++{- |+It must hold @delay <= length a@.+-}+{-+It is crucial that 'mix' uses the chunk size structure of the second operand.+This way we avoid unnecessary and even infinite look-ahead.+-}+{-# INLINE addShiftedSimple #-}+addShiftedSimple ::+   (Additive.C a, SigG.Transform sig a) =>+   Int -> sig a -> sig a -> sig a+addShiftedSimple del a b =+   uncurry CutG.append $+   mapSnd (flip SigG.mix b) $+   CutG.splitAt del a+++-- ** hard-wired convolutions for small sizes++{-+Some small size convolutions using the Karatsuba trick.+We do not use Toom-3 multiplication,+because this requires division by 2 and 6.+With Karatsuba we can stay in a ring.+-}++type Pair a = (a,a)++{- |+Reasonable choices for the multiplication operation are '(*)', '(*>)', 'convolve'.+-}+{-# INLINE convolvePair #-}+convolvePair ::+   (Additive.C a, Additive.C b, Additive.C c) =>+   (a -> b -> c) ->+   Pair a -> Pair b -> Triple c+convolvePair mul a b =+   snd $ sumAndConvolvePair mul a b++{-# INLINE sumAndConvolvePair #-}+sumAndConvolvePair ::+   (Additive.C a, Additive.C b, Additive.C c) =>+   (a -> b -> c) ->+   Pair a -> Pair b -> ((a,b), Triple c)+sumAndConvolvePair (!*!) (a0,a1) (b0,b1) =+   let sa01 = a0+a1+       sb01 = b0+b1+       ab0 = a0!*!b0+       ab1 = a1!*!b1+   in  ((sa01, sb01), (ab0, sa01!*!sb01-(ab0+ab1), ab1))++type Triple a = (a,a,a)++{-# INLINE convolvePairTriple #-}+convolvePairTriple ::+   (Additive.C a, Additive.C b, Additive.C c) =>+   (a -> b -> c) ->+   Pair a -> Triple b -> (c,c,c,c)+convolvePairTriple (!*!) (a0,a1) (b0,b1,b2) =+   let ab0 = a0!*!b0+       ab1 = a1!*!b1+       sa01 = a0+a1; sb01 = b0+b1; ab01 = sa01!*!sb01+   in  (ab0, ab01 - (ab0+ab1),+        a0!*!b2 + ab1, a1!*!b2)+++{-# INLINE convolveTriple #-}+convolveTriple ::+   (Additive.C a, Additive.C b, Additive.C c) =>+   (a -> b -> c) ->+   Triple a -> Triple b -> (c,c,c,c,c)+convolveTriple mul a b =+   snd $ sumAndConvolveTriple mul a b++{-# INLINE sumAndConvolveTriple #-}+sumAndConvolveTriple ::+   (Additive.C a, Additive.C b, Additive.C c) =>+   (a -> b -> c) ->+   Triple a -> Triple b -> ((a,b), (c,c,c,c,c))+sumAndConvolveTriple (!*!) (a0,a1,a2) (b0,b1,b2) =+   let ab0 = a0!*!b0+       ab1 = a1!*!b1+       ab2 = a2!*!b2+       sa01 = a0+a1; sb01 = b0+b1; ab01 = sa01!*!sb01+       sa02 = a0+a2; sb02 = b0+b2; ab02 = sa02!*!sb02+       sa012 = sa01+a2+       sb012 = sb01+b2+   in  ((sa012, sb012),+        (ab0, ab01 - (ab0+ab1),+         ab02 + ab1 - (ab0+ab2),+         sa012!*!sb012 - ab02 - ab01 + ab0, ab2))++{-# INLINE sumAndConvolveTripleAlt #-}+sumAndConvolveTripleAlt ::+   (Additive.C a, Additive.C b, Additive.C c) =>+   (a -> b -> c) ->+   Triple a -> Triple b -> ((a,b), (c,c,c,c,c))+sumAndConvolveTripleAlt (!*!) (a0,a1,a2) (b0,b1,b2) =+   let ab0 = a0!*!b0+       ab1 = a1!*!b1+       ab2 = a2!*!b2+       sa01 = a0+a1; sb01 = b0+b1+       ab01 = sa01!*!sb01 - (ab0+ab1)+       sa02 = a0+a2; sb02 = b0+b2+       ab02 = sa02!*!sb02 - (ab0+ab2)+       sa12 = a1+a2; sb12 = b1+b2+       ab12 = sa12!*!sb12 - (ab1+ab2)+   in  ((sa01+a2, sb01+b2),+        (ab0, ab01, ab1+ab02, ab12, ab2))++type Quadruple a = (a,a,a,a)++{-# INLINE convolveQuadruple #-}+convolveQuadruple ::+   (Additive.C a, Additive.C b, Additive.C c) =>+   (a -> b -> c) ->+   Quadruple a -> Quadruple b -> (c,c,c,c,c,c,c)+convolveQuadruple mul a b =+   snd $ sumAndConvolveQuadruple mul a b++{-# INLINE sumAndConvolveQuadruple #-}+sumAndConvolveQuadruple ::+   (Additive.C a, Additive.C b, Additive.C c) =>+   (a -> b -> c) ->+   Quadruple a -> Quadruple b -> ((a,b), (c,c,c,c,c,c,c))+sumAndConvolveQuadruple (!*!) (a0,a1,a2,a3) (b0,b1,b2,b3) =+   let ab0 = a0!*!b0+       ab1 = a1!*!b1+       sa01 = a0+a1; sb01 = b0+b1+       ab01 = sa01!*!sb01 - (ab0+ab1)+       ab2 = a2!*!b2+       ab3 = a3!*!b3+       sa23 = a2+a3; sb23 = b2+b3+       ab23 = sa23!*!sb23 - (ab2+ab3)+       ab02 = (a0+a2)!*!(b0+b2)+       ab13 = (a1+a3)!*!(b1+b3)+       sa0123 = sa01+sa23+       sb0123 = sb01+sb23+       ab0123 = sa0123!*!sb0123 - (ab02+ab13)+   in  ((sa0123, sb0123),+        (ab0, ab01, ab1+ab02-(ab0+ab2),+         ab0123 - (ab01+ab23),+         ab2+ab13-(ab1+ab3), ab23, ab3))++{-# INLINE sumAndConvolveQuadrupleAlt #-}+sumAndConvolveQuadrupleAlt ::+   (Additive.C a, Additive.C b, Additive.C c) =>+   (a -> b -> c) ->+   Quadruple a -> Quadruple b -> ((a,b), (c,c,c,c,c,c,c))+sumAndConvolveQuadrupleAlt mul (a0,a1,a2,a3) (b0,b1,b2,b3) =+   let (((sa02,sa13), (sb02,sb13)),+        ((c00,c01,c02), (c10,c11,c12), (c20,c21,c22))) =+          sumAndConvolvePair (convolvePair mul)+             ((a0,a1),(a2,a3)) ((b0,b1),(b2,b3))+   in  ((sa02+sa13, sb02+sb13),+        (c00,c01,c02+c10,c11,c12+c20,c21,c22))
src/Synthesizer/Generic/Filter/Recursive/Comb.hs view
@@ -19,11 +19,9 @@ import qualified Synthesizer.Generic.Signal as SigG  import qualified Algebra.Module                as Module--- import qualified Algebra.Field                 as Field import qualified Algebra.Ring                  as Ring import qualified Algebra.Additive              as Additive -import qualified Prelude as P import NumericPrelude.Base import NumericPrelude.Numeric
+ src/Synthesizer/Generic/Fourier.hs view
@@ -0,0 +1,997 @@+{- |+Complete implementation for Fast Fourier Transform for any signal length.+Although defined for all kinds of signal storage,+we need fast access to arbitrary indices.+-}+{-+further thoughts+ - Test the algorithms using remainder polynomials+     with respect to [-1,0,...,0,1]+     Problem: We would need large polynomial degrees,+     namely LCM of the size of all sub-transforms.+     Those numbers are in the same magnitude+     as the integers we use for our integer residue class arithmetic.+ - a z-transform by convolving with a chirp would be nice,+     however we need a square of the primitive root of unity+     in order to compute cis((i/n)^2/2)+ - Can we write the Fourier transforms for lengths larger than the input signal length+     with implicit zero padding?+     This would be useful for Fourier based convolution.+     Our frequent use of 'rechunk' would be a problem, though.+     transformCoprime also needs explicit zero padding.+ - a type class could unify all Level generators+     and thus they would allow for a generic way to call a certain sub-transform+-}+{-# LANGUAGE NoImplicitPrelude #-}+module Synthesizer.Generic.Fourier (+   Element(..),+   -- * conversion between time and frequency domain (spectrum)+   transformForward,+   transformBackward,+   cacheForward,+   cacheBackward,+   cacheDuplex,+   transformWithCache,+   -- * convolution based on Fourier transform+   convolveCyclic,+   Window,+   window,+   convolveWithWindow,+   ) where++import qualified Synthesizer.Generic.Signal as SigG+import qualified Synthesizer.Generic.Cut as CutG+import qualified Synthesizer.Generic.Cyclic as Cyclic+import qualified Synthesizer.Generic.Filter.NonRecursive as FiltNRG++import qualified Synthesizer.Generic.Permutation as Permutation+import qualified Synthesizer.Basic.NumberTheory as NumberTheory++import qualified Synthesizer.State.Analysis as Ana+import qualified Synthesizer.State.Signal as SigS++import qualified Algebra.Transcendental as Trans+-- import qualified Algebra.Field as Field+import qualified Algebra.Ring as Ring+import qualified Algebra.PrincipalIdealDomain as PID+import qualified Algebra.IntegralDomain as Integral++import qualified Number.ResidueClass.Check as RC+import Number.ResidueClass.Check ((/:), )++import qualified Number.Complex as Complex+import Number.Complex ((+:))++import qualified Data.Map as Map+import qualified Control.Monad.Trans.State as State+import Control.Monad (liftM2, )+import Control.Applicative ((<$>), )++import Data.Tuple.HT (mapPair, )++import NumericPrelude.Numeric+import NumericPrelude.Base hiding (head, )++++class Ring.C y => Element y where+   recipInteger :: (SigG.Read sig y) => sig y -> y+   addId :: (SigG.Read sig y) => sig y -> y+   multId :: (SigG.Read sig y) => sig y -> y+   {- |+   It must hold:++   > uncurry (*) (conjugatePrimitiveRootsOfUnity n) = 1++   > mapPair ((^m), (^m)) (conjugatePrimitiveRootsOfUnity (n*m) y)+   >    == conjugatePrimitiveRootsOfUnity n y@++   since we need for caching that the cache is uniquely determined+   by singal length and transform direction.+   -}+   conjugatePrimitiveRootsOfUnity :: (SigG.Read sig y) => sig y -> (y,y)++instance Trans.C a => Element (Complex.T a) where+   recipInteger sig = recip (fromIntegral (SigG.length sig)) +: zero+   addId _sig = zero+   multId _sig = one+   conjugatePrimitiveRootsOfUnity sig =+      (\x -> (x, Complex.conjugate x)) $+      case SigG.length sig of+         1 -> one+         2 -> negate one+         3 -> (negate one +: sqrt 3) / 2+         4 -> zero +: one+         5 ->+            let sqrt5 = sqrt 5+            in  ((sqrt5 - 1) +: sqrt 2 * sqrt(5 + sqrt5)) / 4+         6 -> (one +: sqrt 3) / 2+         8 -> Complex.scale (sqrt 2 / 2) (one +: one)+         12 -> (sqrt 3 +: one) / 2+         n -> Complex.cis (2*pi / fromIntegral n)++instance (NumberTheory.PrimitiveRoot a, PID.C a, Eq a) => Element (RC.T a) where+   recipInteger sig =+      recip (fromIntegral (SigG.length sig) /: RC.modulus (head sig))+   addId sig = zero /: RC.modulus (head sig)+   multId sig = one /: RC.modulus (head sig)+   {-+   We cannot simply compute+     NumberTheory.primitiveRootsOfUnity modu (SigG.length sig)+   since we have to fulfill the laws.+   In order to fulfill them,+   we choose a root with maximum order,+   this will always be the same,+   and it is a root of all primitive roots+   of any possible order in that ring.+   -}+   conjugatePrimitiveRootsOfUnity sig =+      let modu = RC.modulus (head sig)+          order@(NumberTheory.Order expo) =+             NumberTheory.maximumOrderOfPrimitiveRootsOfUnity modu+          r:_ = NumberTheory.primitiveRootsOfUnity modu order+          n = Integral.divChecked expo (fromIntegral (SigG.length sig))+          z = (r /: modu) ^ n+      in  (z, recip z)+++head :: (SigG.Read sig y) => sig y -> y+head =+   SigG.switchL (error "Generic.Signal.head: empty signal") const .+   SigG.toState+++directionPrimitiveRootsOfUnity ::+   (Element y, SigG.Read sig y) =>+   sig y -> ((Direction,y), (Direction,y))+directionPrimitiveRootsOfUnity x =+   let (z,zInv) =+          conjugatePrimitiveRootsOfUnity x+   in  ((Forward,z), (Backward,zInv))++transformForward ::+   (Element y, SigG.Transform sig y) =>+   sig y -> sig y+transformForward xs =+   transformWithCache (cacheForward xs) xs++{- |+Shall we divide the result values by the length of the signal?+Our dimensional wrapper around the Fourier transform does not expect this.+-}+transformBackward ::+   (Element y, SigG.Transform sig y) =>+   sig y -> sig y+transformBackward xs =+   transformWithCache (cacheBackward xs) xs++{- |+The size of the signal must match the size, that the plan was generated for.+-}+_transformPlan ::+   (Element y, SigG.Transform sig y) =>+   Plan -> (Direction,y) -> sig y -> sig y+_transformPlan p z xs =+   transformWithCache (cacheFromPlan p z xs) xs++{- |+The size and type of the signal must match the parameters,+that the cache was generated for.+-}+transformWithCache ::+   (Element y, SigG.Transform sig y) =>+   Cache sig y -> sig y -> sig y+transformWithCache cache xs =+   case cache of+      CacheIdentity -> xs+      CacheSmall size ->+         case size of+            LevelCache2 zs -> transform2 zs xs+            LevelCache3 zs -> transform3 zs xs+            LevelCache4 zs -> transform4 zs xs+            LevelCache5 zs -> transform5 zs xs+      CacheNaive level ->+         transformNaive level xs+      CacheRadix2 level subCache ->+         transformRadix2InterleavedFrequency level subCache xs+      CachePrime level subCaches ->+         transformPrime level subCaches xs+      CacheCoprime level subCaches ->+         transformCoprime level subCaches xs+      CacheComposite level subCaches ->+         transformComposite level subCaches xs+++{- |+Memorize factorizations of the data size and permutation vectors.+-}+data Plan =+     PlanIdentity+   | PlanSmall LevelSmall+   | PlanNaive  -- mainly for debugging+   | PlanRadix2 LevelRadix2 Plan+   | PlanPrime LevelPrime Plan+   | PlanCoprime LevelCoprime (Plan, Plan)+   | PlanComposite LevelComposite (Plan, Plan)+   deriving (Show)++{-+efficient swallow comparison+only correct for Plans generated by 'plan'.+-}+instance Eq Plan where+   p0 == p1  =  compare p0 p1 == EQ++{-+Needed for keys in CacheMap+-}+instance Ord Plan where+   compare p0 p1  =+      case (p0,p1) of+         (PlanIdentity, PlanIdentity) -> EQ+         (PlanIdentity, _) -> LT+         (_, PlanIdentity) -> GT+         (PlanSmall l0, PlanSmall l1) -> compare l0 l1+         (PlanSmall _, _) -> LT+         (_, PlanSmall _) -> GT+         (PlanNaive, PlanNaive) -> EQ+         (PlanNaive, _) -> LT+         (_, PlanNaive) -> GT+         (PlanRadix2 l0 _, PlanRadix2 l1 _) -> compare l0 l1+         (PlanRadix2 _ _, _) -> LT+         (_, PlanRadix2 _ _) -> GT+         (PlanPrime l0 _, PlanPrime l1 _) -> compare l0 l1+         (PlanPrime _ _, _) -> LT+         (_, PlanPrime _ _) -> GT+         (PlanCoprime l0 _, PlanCoprime l1 _) -> compare l0 l1+         (PlanCoprime _ _, _) -> LT+         (_, PlanCoprime _ _) -> GT+         (PlanComposite l0 _, PlanComposite l1 _) -> compare l0 l1+++plan :: Integer -> Plan+plan n =+   State.evalState (planWithMapUpdate n) smallPlanMap++type PlanMap = Map.Map Integer Plan++smallPlanMap :: PlanMap+smallPlanMap =+   Map.fromAscList $ zip [0..] $+   PlanIdentity :+   PlanIdentity :+   PlanSmall Level2 :+   PlanSmall Level3 :+   PlanSmall Level4 :+   PlanSmall Level5 :+   []++{- |+Detect and re-use common sub-plans.+-}+planWithMap :: Integer -> State.State PlanMap Plan+planWithMap n =+   case divMod n 2 of+      (n2,0) -> PlanRadix2 (levelRadix2 n2) <$> planWithMapUpdate n2+      _ ->+         let facs = NumberTheory.fermatFactors n+         in  -- find unitary divisors+             case filter (\(a,b) -> a>1 && gcd a b == 1) facs of+                q2 : _ ->+                   PlanCoprime (levelCoprime q2) <$>+                   planWithMapUpdate2 q2+                _ ->+                   let (q2 : _) = facs+                   in  if fst q2 == 1+                         then PlanPrime (levelPrime $ snd q2) <$>+                              planWithMapUpdate (n-1)+                         else PlanComposite (levelComposite q2) <$>+                              planWithMapUpdate2 q2++planWithMapUpdate :: Integer -> State.State PlanMap Plan+planWithMapUpdate n = do+   item <- State.gets (Map.lookup n)+   case item of+      Just p -> return p+      Nothing ->+         planWithMap n >>= \m -> State.modify (Map.insert n m) >> return m++planWithMapUpdate2 :: (Integer, Integer) -> State.State PlanMap (Plan, Plan)+planWithMapUpdate2 =+   uncurry (liftM2 (,)) .+   mapPair (planWithMapUpdate,planWithMapUpdate)+++{- |+Cache powers of the primitive root of unity+in a storage compatible to the processed signal.+-}+data Cache sig y =+     CacheIdentity+   | CacheSmall (LevelCacheSmall y)+   | CacheNaive (LevelCacheNaive y)+   | CacheRadix2 (LevelCacheRadix2 sig y) (Cache sig y)+   | CachePrime (LevelCachePrime sig y) (Cache sig y, Cache sig y)+   | CacheCoprime LevelCoprime (Cache sig y, Cache sig y)+   | CacheComposite (LevelCacheComposite sig y) (Cache sig y, Cache sig y)+   deriving (Show)++{- |+The expression @cacheForward prototype@+precomputes all data that is needed for forward Fourier transforms+for signals of the type and length @prototype@.+You can use this cache in 'transformWithCache'.+-}+cacheForward ::+   (Element y, SigG.Transform sig y) =>+   sig y -> Cache sig y+cacheForward xs =+   cacheFromPlan+      (plan $ fromIntegral $ SigG.length xs)+      (fst $ directionPrimitiveRootsOfUnity xs)+      xs++{- |+See 'cacheForward'.+-}+cacheBackward ::+   (Element y, SigG.Transform sig y) =>+   sig y -> Cache sig y+cacheBackward xs =+   cacheFromPlan+      (plan $ fromIntegral $ SigG.length xs)+      (snd $ directionPrimitiveRootsOfUnity xs)+      xs++{- |+It is @(cacheForward x, cacheBackward x) = cacheDuplex x@+but 'cacheDuplex' shared common data of both caches.+-}+cacheDuplex ::+   (Element y, SigG.Transform sig y) =>+   sig y -> (Cache sig y, Cache sig y)+cacheDuplex xs =+   let p = plan $ fromIntegral $ SigG.length xs+       (z,zInv) = directionPrimitiveRootsOfUnity xs+   in  State.evalState+          (cacheFromPlanWithMapUpdate2 (p,p) (z,zInv) (xs,xs)) $+       Map.empty+++data Direction = Forward | Backward+   deriving (Show, Eq, Ord)++type CacheMap sig y = Map.Map (Plan,Direction) (Cache sig y)++cacheFromPlan ::+   (Element y, SigG.Transform sig y) =>+   Plan -> (Direction, y) -> sig y -> Cache sig y+cacheFromPlan p z xs =+   State.evalState (cacheFromPlanWithMapUpdate p z xs) $+   Map.empty++{- |+Detect and re-use common sub-caches.+-}+cacheFromPlanWithMap ::+   (Element y, SigG.Transform sig y) =>+   Plan -> (Direction,y) -> sig y ->+   State.State (CacheMap sig y) (Cache sig y)+cacheFromPlanWithMap p (d,z) xs =+   case p of+      PlanIdentity -> return $ CacheIdentity+      PlanSmall size -> return $ CacheSmall $+         case size of+            Level2 -> LevelCache2 $ cache2 z+            Level3 -> LevelCache3 $ cache3 z+            Level4 -> LevelCache4 $ cache4 z+            Level5 -> LevelCache5 $ cache5 z+      PlanNaive ->+         return $ CacheNaive $ LevelCacheNaive z+      PlanRadix2 level@(LevelRadix2 n2) subPlan ->+         let subxs = CutG.take n2 xs+         in  CacheRadix2 (levelCacheRadix2 level z subxs) <$>+             cacheFromPlanWithMapUpdate subPlan (d,z*z) subxs+      PlanPrime level@(LevelPrime (perm,_,_)) subPlan ->+         (\subCaches ->+            CachePrime+               (levelCachePrime level (fst subCaches) z xs) subCaches)+         <$>+         let subxs = CutG.take (Permutation.size perm) xs+         in  cacheFromPlanWithMapUpdate2 (subPlan,subPlan)+                (directionPrimitiveRootsOfUnity subxs)+                (subxs,subxs)+      PlanCoprime level@(LevelCoprime (n,m) _) subPlans ->+         CacheCoprime level <$>+         cacheFromPlanWithMapUpdate2 subPlans ((d,z^m), (d,z^n))+            (CutG.take (fromInteger n) xs, CutG.take (fromInteger m) xs)+      PlanComposite level@(LevelComposite (n,m) _) subPlans ->+         CacheComposite (levelCacheComposite level z xs) <$>+         cacheFromPlanWithMapUpdate2 subPlans ((d,z^m), (d,z^n))+            (CutG.take (fromInteger n) xs, CutG.take (fromInteger m) xs)++cacheFromPlanWithMapUpdate ::+   (Element y, SigG.Transform sig y) =>+   Plan -> (Direction,y) -> sig y ->+   State.State (CacheMap sig y) (Cache sig y)+cacheFromPlanWithMapUpdate p z xs = do+   let key = (p, fst z)+   item <- State.gets (Map.lookup key)+   case item of+      Just c -> return c+      Nothing ->+         cacheFromPlanWithMap p z xs >>= \m ->+         State.modify (Map.insert key m) >>+         return m++cacheFromPlanWithMapUpdate2 ::+   (Element y, SigG.Transform sig y) =>+   (Plan, Plan) -> ((Direction,y),(Direction,y)) -> (sig y, sig y) ->+   State.State (CacheMap sig y) (Cache sig y, Cache sig y)+cacheFromPlanWithMapUpdate2 (p0,p1) (z0,z1) (xs0,xs1) =+   liftM2 (,)+      (cacheFromPlanWithMapUpdate p0 z0 xs0)+      (cacheFromPlanWithMapUpdate p1 z1 xs1)+++newtype LevelCacheNaive y =+      LevelCacheNaive y+   deriving (Show)++transformNaive ::+   (Element y, SigG.Transform sig y) =>+   LevelCacheNaive y -> sig y -> sig y+transformNaive (LevelCacheNaive z) sig =+   SigG.takeStateMatch sig $+   SigS.map+      (scalarProduct1 (SigG.toState sig) . powers sig)+      (powers sig z)++scalarProduct1 ::+   (Ring.C a) =>+   SigS.T a -> SigS.T a -> a+scalarProduct1 xs ys =+   SigS.foldL1 (+) $ SigS.zipWith (*) xs ys++_transformRing ::+   (Ring.C y, SigG.Transform sig y) =>+   y -> sig y -> sig y+_transformRing z sig =+   SigG.takeStateMatch sig $+   Ana.chirpTransform z $ SigG.toState sig++powers ::+   (Element y, SigG.Read sig y) =>+   sig y -> y -> SigS.T y+powers sig c = SigS.iterate (c*) $ multId sig+++data LevelSmall = Level2 | Level3 | Level4 | Level5+   deriving (Show, Eq, Ord, Enum)++data LevelCacheSmall y =+     LevelCache2 y+   | LevelCache3 (y,y)+   | LevelCache4 (y,y,y)+   | LevelCache5 (y,y,y,y)+   deriving (Show)++cache2 :: (Ring.C y) => y -> y+cache3 :: (Ring.C y) => y -> (y,y)+cache4 :: (Ring.C y) => y -> (y,y,y)+cache5 :: (Ring.C y) => y -> (y,y,y,y)++cache2 z = z+cache3 z = (z, z*z)+cache4 z = let z2=z*z in (z,z2,z*z2)+cache5 z = let z2=z*z in (z,z2,z*z2,z2*z2)+++transform2 ::+   (Ring.C y, SigG.Transform sig y) =>+   y -> sig y -> sig y+transform2 z sig =+   let x0:x1:_ = SigG.toList sig+   in  SigG.takeStateMatch sig $+       SigS.fromList [x0+x1, x0+z*x1]++transform3 ::+   (Ring.C y, SigG.Transform sig y) =>+   (y,y) -> sig y -> sig y+transform3 (z,z2) sig =+   let x0:x1:x2:_ = SigG.toList sig+{- Rader's algorithm with convolution by 2-size-Fourier-transform+       xf1 = x1+x2+       xf2 = x1-x2+       zf1 = z+z2+       zf2 = z-z2+       xzf1 = xf1*zf1+       xzf2 = xf2*zf2+       xz1 = (xzf1+xzf2)/2+       xz2 = (xzf1-xzf2)/2+-}+{- naive+       [x0+x1+x2, x0+z*x1+z2*x2, x0+z2*x1+z*x2]+-}+       ((s,_), (zx1,zx2)) = Cyclic.sumAndConvolvePair (x1,x2) (z,z2)+   in  SigG.takeStateMatch sig $+       SigS.fromList [x0+s, x0+zx1, x0+zx2]++transform4 ::+   (Ring.C y, SigG.Transform sig y) =>+   (y,y,y) -> sig y -> sig y+transform4 (z,z2,z3) sig =+   let x0:x1:x2:x3:_ = SigG.toList sig+       x02a = x0+x2; x02b = x0+z2*x2+       x13a = x1+x3; x13b = x1+z2*x3+   in  SigG.takeStateMatch sig $+       SigS.fromList [x02a+   x13a, x02b+z *x13b,+                      x02a+z2*x13a, x02b+z3*x13b]+{-+This needs also five multiplications,+but in complex numbers it is z=i, and thus multiplications are cheap+and we should better make use of distributive law in order to save additions.++       x02a = x0+x2; x02b = x0+z2*x2+       x1_2 = z2*x1; x3_2 = z2*x3+   in  SigG.takeStateMatch sig $+       SigS.fromList [x02a + x1   + x3  , x02b+z*(x1   + x3_2),+                      x02a + x1_2 + x3_2, x02b+z*(x1_2 + x3  )]+-}++{-+Use Rader's trick for mapping the transform to a convolution+and apply Karatsuba's trick at two levels (i.e. total three times)+to that convolution.++0 0 0 0 0+0 1 2 3 4+0 2 4 1 3+0 3 1 4 2+0 4 3 2 1++Permutation.T: 0 1 2 4 3++0 0 0 0 0+0 1 2 4 3+0 2 4 3 1+0 4 3 1 2+0 3 1 2 4+-}+transform5 ::+   (Ring.C y, SigG.Transform sig y) =>+   (y,y,y,y) -> sig y -> sig y+transform5 (z1,z2,z3,z4) sig =+   let x0:x1:x2:x3:x4:_ = SigG.toList sig+       ((s,_), (d1,d2,d4,d3)) =+          Cyclic.sumAndConvolveQuadruple (x1,x3,x4,x2) (z1,z2,z4,z3)+   in  SigG.takeStateMatch sig $+       SigS.fromList [x0+s, x0+d1, x0+d2, x0+d3, x0+d4]++{-+transform7++Toom-3-multiplication at the highest level and Karatsuba below?+Toom-2.5-multiplication with manual addition of the missing parts?++Toom-3-multiplication with complex interpolation nodes?+Still requires division by 4 and then complex multiplication in the frequency domain.+A:=matrix(5,5,[1,0,0,0,0,1,1,1,1,1,1,-1,1,-1,1,1,I,-1,-I,1,0,0,0,0,1]);+A:=matrix(5,5,[1,0,0,0,0,1,1,1,1,1,1,-1,1,-1,1,1,I,-1,-I,1,1,-I,-1,I,1]);++Karatsuba at three levels for convolution of signal of size 8 with zero padding?++Modify the 3x3 Fourier matrix by multiplying a regular matrix+to make it more convenient to work with?+We will hardly get rid of the irrational numbers.+-}++newtype LevelRadix2 = LevelRadix2 Int+   deriving (Show, Eq, Ord)++levelRadix2 :: Integer -> LevelRadix2+levelRadix2 =+   LevelRadix2 . fromIntegral+++data LevelCacheRadix2 sig y =+   LevelCacheRadix2 Int (sig y)+   deriving (Show)++levelCacheRadix2 ::+   (Element y, SigG.Transform sig y) =>+   LevelRadix2 -> y -> sig y -> LevelCacheRadix2 sig y+levelCacheRadix2 (LevelRadix2 n2) z sig =+   LevelCacheRadix2 n2+      (SigG.takeStateMatch sig $ powers sig z)+++{- |+Cooley-Tukey specialised to one factor of the size being 2.++Size of the input signal must be even.+-}+transformRadix2InterleavedFrequency ::+   (Element y, SigG.Transform sig y) =>+   LevelCacheRadix2 sig y -> Cache sig y -> sig y -> sig y+transformRadix2InterleavedFrequency+      (LevelCacheRadix2 n2 twiddle) subCache sig =+   let (xs0,xs1) = SigG.splitAt n2 sig+       fs0 = transformWithCache subCache $ SigG.zipWith (+) xs0 xs1+       fs1 = transformWithCache subCache $+                SigG.zipWith3+                   (\w x0 x1 -> w*(x0-x1))+                   twiddle xs0 xs1+   in  SigG.takeStateMatch sig $+       SigS.interleave (SigG.toState fs0) (SigG.toState fs1)+++data LevelComposite =+   LevelComposite+      (Integer, Integer)+      (Permutation.T, Permutation.T)+   deriving (Show)++instance Eq LevelComposite where+   a == b  =  compare a b == EQ++instance Ord LevelComposite where+   compare (LevelComposite a _) (LevelComposite b _)  =+      compare a b++levelComposite :: (Integer, Integer) -> LevelComposite+levelComposite (n,m) =+   let ni = fromInteger n+       mi = fromInteger m+   in  LevelComposite (n,m)+          (Permutation.transposition ni mi,+           Permutation.transposition mi ni)+++data LevelCacheComposite sig y =+   LevelCacheComposite+      (Integer, Integer)+      (Permutation.T, Permutation.T)+      (sig y)+   deriving (Show)++levelCacheComposite ::+   (Element y, SigG.Transform sig y) =>+   LevelComposite -> y -> sig y -> LevelCacheComposite sig y+levelCacheComposite (LevelComposite (n,m) transpose) z sig =+   LevelCacheComposite (n,m) transpose $+   SigG.takeStateMatch sig $+   flip SigS.generateInfinite (n, multId sig, multId sig) $ \(i,zi,zij) ->+   (zij,+    case pred i of+      0 -> (n, zi*z, multId sig)+      i1 -> (i1, zi, zij*zi))+{-+   {-# SCC "levelCacheComposite:rechunk" #-}+   concatRechunk sig $+   {-# SCC "levelCacheComposite:subpowers" #-}+   SigS.map+      (SigG.takeStateMatch (SigG.take (fromIntegral n) sig) . powers sig)+      ({-# SCC "levelCacheComposite:powers" #-}+       powers sig z)+-}+{-+   SigS.map+      (SigG.takeStateMatch sig . SigS.take (fromIntegral n) . powers sig)+      ({-# SCC "levelCacheComposite:powers" #-}+       powers sig z)+-}+{- suffers from big inefficiency of repeated 'append'+   SigG.takeStateMatch sig $+   SigS.monoidConcat $+   SigS.map (SigS.take (fromIntegral n) . powers sig) $+   SigS.take (fromIntegral m) $ -- necessary for strict storable vectors+   powers sig z+-}++{- |+For @transformComposite z (n,m) sig@,+the parameters @n@ and @m@ must be relatively prime+and @n*m == length sig@ and @z ^ length sig == 1@.++Cooley-Tukey-algorithm+-}+transformComposite ::+   (Element y, SigG.Transform sig y) =>+   LevelCacheComposite sig y -> (Cache sig y, Cache sig y) -> sig y -> sig y+transformComposite+      (LevelCacheComposite (n,m) (transposeNM, transposeMN) twiddle)+      (subCacheN,subCacheM) sig =+   Permutation.apply transposeMN .+       concatRechunk sig .+       SigS.map (transformWithCache subCacheM) .+       SigG.sliceVertical (fromInteger m) .+       Permutation.apply transposeNM .+--       concatRechunk sig .+       SigG.zipWith (*) twiddle .+       SigS.monoidConcat .+       SigS.map (transformWithCache subCacheN) .+       SigG.sliceVertical (fromInteger n) .+       Permutation.apply transposeMN $+       sig+++data LevelCoprime =+   LevelCoprime+      (Integer, Integer)+      (Permutation.T, Permutation.T, Permutation.T)+   deriving (Show)++instance Eq LevelCoprime where+   a == b  =  compare a b == EQ++instance Ord LevelCoprime where+   compare (LevelCoprime a _) (LevelCoprime b _)  =+      compare a b++{-+Fourier exponent matrix of a signal of size 6.++0 0 0 0 0 0     0               0   0   0       0     0+0 1 2 3 4 5               0       2   0   4       3     0+0 2 4 0 2 4  =          0    *  0   2   4    *      0     0+0 3 0 3 0 3           0           0   0   0     0     3+0 4 2 0 4 2         0           0   4   2         0     0+0 5 4 3 2 1       0               4   0   2         0     3+-}+levelCoprime :: (Integer, Integer) -> LevelCoprime+levelCoprime (n,m) =+   let ni = fromInteger n+       mi = fromInteger m+   in  LevelCoprime (n,m)+          (Permutation.skewGrid mi ni,+           Permutation.transposition ni mi,+           Permutation.skewGridCRTInv ni mi)+++{- |+For @transformCoprime z (n,m) sig@,+the parameters @n@ and @m@ must be relatively prime+and @n*m == length sig@ and @z ^ length sig == 1@.++Good-Thomas algorithm+-}+{-+A very elegant way would be to divide the signal into chunks of size n,+define ring operations on these chunks+and perform one (length/n)-size-sub-transform in this chunk-ring.+This way we would also only have to plan the sub-transform once.+On StorableVectors the chunking could be performed in-place+in terms of a virtual reshape operation.+In the general case the performance can become very bad+if the chunks are very small, say 2 or 3 elements.+-}+transformCoprime ::+   (Element y, SigG.Transform sig y) =>+   LevelCoprime -> (Cache sig y, Cache sig y) -> sig y -> sig y+transformCoprime+      (LevelCoprime (n,m) (grid, transpose, gridInv)) (subCacheN,subCacheM) =+   let subTransform cache j sig =+          concatRechunk sig .+          SigS.map (transformWithCache cache) .+          SigG.sliceVertical (fromIntegral j) $ sig+   in  Permutation.apply gridInv .+       subTransform subCacheM m .+       Permutation.apply transpose .+       subTransform subCacheN n .+       Permutation.apply grid+++-- concatenate and reorganize for faster indexing+concatRechunk ::+   (SigG.Transform sig y) =>+   sig y -> SigS.T (sig y) -> sig y+concatRechunk pattern =+   SigG.takeStateMatch pattern .+   SigG.toState .+   SigS.monoidConcat+++data LevelPrime =+   LevelPrime (Permutation.T, Permutation.T, Permutation.T)+      deriving (Show)++instance Eq LevelPrime where+   a == b  =  compare a b == EQ++instance Ord LevelPrime where+   compare (LevelPrime (a,_,_)) (LevelPrime (b,_,_))  =+      compare (Permutation.size a) (Permutation.size b)++{-+Fourier exponent matrix of a signal of size 7.++0 0 0 0 0 0 0+0 1 2 3 4 5 6+0 2 4 6 1 3 5+0 3 6 2 5 1 4+0 4 1 5 2 6 3+0 5 3 1 6 4 2+0 6 5 4 3 2 1++multiplicative generator in Z7: 3+permutation of rows and columns by powers of 3: 1 3 2 6 4 5++0 0 0 0 0 0 0+0 1 3 2 6 4 5+0 3 2 6 4 5 1+0 2 6 4 5 1 3+0 6 4 5 1 3 2+0 4 5 1 3 2 6+0 5 1 3 2 6 4++Inverse permutation: 1 3 2 5 6 4+The inverse permutations seems not to be generated by a multiplication.+-}+levelPrime :: Integer -> LevelPrime+levelPrime n =+   let perm = Permutation.multiplicative $ fromIntegral n+   in  LevelPrime+          (perm, Permutation.reverse perm, Permutation.inverse perm)+++data LevelCachePrime sig y =+   LevelCachePrime (Permutation.T, Permutation.T) (sig y)+      deriving (Show)++levelCachePrime ::+   (Element y, SigG.Transform sig y) =>+   LevelPrime -> Cache sig y -> y -> sig y -> LevelCachePrime sig y+levelCachePrime (LevelPrime (perm, rev, inv)) subCache z sig =+   LevelCachePrime (rev, inv)+      ((\zs -> FiltNRG.amplify (recipInteger zs) zs) $+       transformWithCache subCache $+       Permutation.apply perm $+       SigG.takeStateMatch sig $+       SigS.iterate (z*) z)++{- |+Rader's algorithm for prime length signals.+-}+transformPrime ::+   (Element y, SigG.Transform sig y) =>+   LevelCachePrime sig y -> (Cache sig y, Cache sig y) -> sig y -> sig y+transformPrime (LevelCachePrime (rev, inv) zs) subCaches =+   SigG.switchL (error "transformPrime: empty signal") $+   \x0 rest ->+      SigG.cons (SigG.foldL (+) x0 rest) $+      SigG.map (x0+) $+      Permutation.apply inv $+      convolveSpectrumCyclicCache subCaches zs $+      Permutation.apply rev rest++{-+Cyclic.reverse xs = shiftR 1 (reverse xs)+Cyclic.reverse (xs <*> ys) = Cyclic.reverse xs <*> Cyclic.reverse ys+Cyclic.reverse (Cyclic.reverse xs) = xs++We could move the 'Cyclic.reverse' over to the z-vector,+but then we would have to reverse again after convolution.++zs <*> Cyclic.reverse rest+ = Cyclic.reverse (Cyclic.reverse zs <*> rest)+-}++{-+This uses Cyclic.filter instead of Cyclic.convolve.+This is simpler, but Fourier.convolveCyclic is a bit simpler than Fourier.filterCyclic,+since it does not need to reverse an operand.+-}+_transformPrimeAlt ::+   (Ring.C y, SigG.Transform sig y) =>+   LevelPrime -> y -> sig y -> sig y+_transformPrimeAlt (LevelPrime (perm, _, inv)) z =+   SigG.switchL (error "transformPrime: empty signal") $+   \x0 rest ->+      SigG.cons (SigG.foldL (+) x0 rest) $+      SigG.map (x0+) $+      Permutation.apply inv $+      Cyclic.filterNaive+         (Permutation.apply perm rest)+         (Permutation.apply perm (SigG.takeStateMatch rest (SigS.iterate (z*) z)))++++{- |+Filter window stored as spectrum+such that it can be applied efficiently to long signals.+-}+data Window sig y =+   Window Int (Cache sig y, Cache sig y) (sig y)+   deriving (Show)+++window ::+   (Element y, SigG.Transform sig y) =>+   sig y -> Window sig y+window x =+   if CutG.null x+     then Window 0 (CacheIdentity, CacheIdentity) CutG.empty+     else+       let size  = CutG.length x+           size2 = 2 * NumberTheory.ceilingPowerOfTwo size+           padded =+              SigG.take size2 $+              CutG.append x $+                 let pad = SigG.takeStateMatch x $ SigS.repeat $ addId x+                 in  CutG.append pad (SigG.append pad pad)+           caches@(cache, _cacheInv) =+              cacheDuplex padded+       in  Window+              (size2-size+1)+              caches+              (transformWithCache cache $+               FiltNRG.amplify (recipInteger padded) padded)++{- |+Efficient convolution of a large filter window+with a probably infinite signal.+-}+convolveWithWindow ::+   (Element y, SigG.Transform sig y) =>+   Window sig y -> sig y -> sig y+convolveWithWindow (Window blockSize caches spectrum) b =+   if blockSize==zero+     then CutG.empty+     else+       let windowSize = SigG.length spectrum - blockSize+       in  SigS.foldR (FiltNRG.addShiftedSimple blockSize) CutG.empty $+           SigS.map+              (\block ->+                 SigG.take (windowSize + SigG.length block) $+                 convolveSpectrumCyclicCache caches spectrum $+                 flip CutG.append+                    {-+                    The last block may be shorter than blockSize+                    and thus needs more padding.+                    -}+                    (SigG.takeStateMatch spectrum $ SigS.repeat $ addId b) $+                 block) $+           SigG.sliceVertical blockSize b+++{- |+Signal must have equal size and must not be empty.+-}+convolveCyclic ::+   (Element y, SigG.Transform sig y) =>+   sig y -> sig y -> sig y+convolveCyclic x =+   let len = fromIntegral $ SigG.length x+       (z,zInv) =+          directionPrimitiveRootsOfUnity x+   in  convolveCyclicCache+          (cacheFromPlan (plan len) z x,+           cacheFromPlan (plan len) zInv x)+          x++convolveCyclicCache ::+   (Element y, SigG.Transform sig y) =>+   (Cache sig y, Cache sig y) -> sig y -> sig y -> sig y+convolveCyclicCache caches x =+   convolveSpectrumCyclicCache caches $+   FiltNRG.amplify (recipInteger x) $ transformWithCache (fst caches) x++{- |+This function does not apply scaling.+That is you have to scale the spectrum by @recip (length x)@+if you want a plain convolution.+-}+convolveSpectrumCyclicCache ::+   (Element y, SigG.Transform sig y) =>+   (Cache sig y, Cache sig y) -> sig y -> sig y -> sig y+convolveSpectrumCyclicCache (cache,cacheInv) x y =+   transformWithCache cacheInv $+   SigG.zipWith (*) x $+   transformWithCache cache y++{-+Test:++let xs = [0,1,0,0,0,0 :: Complex.T Double]; z = fst $ conjugatePrimitiveRootsOfUnity xs in print (transformNaive z xs) >> print (transformCoprime z (2,3) xs)+-}
src/Synthesizer/Generic/Interpolation.hs view
@@ -21,11 +21,8 @@ import qualified Algebra.Module    as Module import qualified Algebra.RealField as RealField import qualified Algebra.RealRing  as RealRing--- import qualified Algebra.Field     as Field--- import qualified Algebra.Ring      as Ring import qualified Algebra.Additive  as Additive -import Algebra.Additive(zero, ) import Data.Maybe (fromMaybe, )  import NumericPrelude.Base
+ src/Synthesizer/Generic/LengthSignal.hs view
@@ -0,0 +1,62 @@+{-# LANGUAGE NoImplicitPrelude #-}+{-# LANGUAGE FlexibleInstances #-}+module Synthesizer.Generic.LengthSignal where++import qualified Synthesizer.Generic.Signal as SigG+import qualified Synthesizer.Generic.Cut as CutG++import qualified Algebra.Additive       as Additive++import Data.Monoid (Monoid, mempty, mappend, )+import Data.Tuple.HT (mapSnd, )++import NumericPrelude.Numeric as NP+import NumericPrelude.Base hiding (length, splitAt, )+++data T sig = Cons {length :: Int, body :: sig}+   deriving (Show)++fromSignal :: (CutG.Read sig) => sig -> T sig+fromSignal xs  =  Cons (CutG.length xs) xs++toSignal :: T sig -> sig+toSignal  =  body++{- |+Each fmap must preserve the signal length.+-}+instance Functor T where+   fmap f (Cons xl xs) = Cons xl (f xs)++instance (Additive.C a, SigG.Transform sig a) => Additive.C (T (sig a)) where+   zero = mempty+   negate xs = xs{body = SigG.map negate (body xs)}+   (Cons xl xs) + (Cons yl ys) =+      Cons (max xl yl) (SigG.mix xs ys)++instance (Monoid sig) => Monoid (T sig) where+   mempty = Cons zero mempty+   mappend (Cons xl xs) (Cons yl ys) =+      Cons (xl+yl) (mappend xs ys)++splitAt :: (CutG.Transform sig) => Int -> T sig -> (T sig, T sig)+splitAt n (Cons xl xs) =+   let (ys,zs) = SigG.splitAt n xs+   in  (Cons (min n xl) ys, Cons (max n xl - n) zs)++{- |+It must hold @delay <= length a@.+-}+{-+It is crucial that 'mix' uses the chunk size structure of the second operand.+This way we avoid unnecessary and even infinite look-ahead.+-}+{-# INLINE addShiftedSimple #-}+addShiftedSimple ::+   (Additive.C a, SigG.Transform sig a) =>+   Int -> T (sig a) -> T (sig a) -> T (sig a)+addShiftedSimple del a b =+   uncurry mappend $+   mapSnd (flip (+) b) $+   splitAt del a
src/Synthesizer/Generic/Loop.hs view
@@ -30,9 +30,7 @@ import qualified Algebra.Transcendental as Trans import qualified Algebra.Module         as Module import qualified Algebra.RealField      as RealField-import qualified Algebra.RealRing           as RealRing-import qualified Algebra.Ring           as Ring-import qualified Algebra.Additive       as Additive+import qualified Algebra.RealRing       as RealRing  import NumericPrelude.Numeric import NumericPrelude.Base
src/Synthesizer/Generic/Noise.hs view
@@ -13,9 +13,7 @@ import qualified Algebra.Ring                  as Ring  import System.Random (Random, RandomGen, randomR, mkStdGen, )-import qualified System.Random as Rnd -import qualified Prelude as P import NumericPrelude.Base import NumericPrelude.Numeric
src/Synthesizer/Generic/Oscillator.hs view
@@ -29,24 +29,10 @@  import Control.Arrow ((>>>), ) -{--import qualified Algebra.RealTranscendental    as RealTrans-import qualified Algebra.Module                as Module-import qualified Algebra.VectorSpace           as VectorSpace--import Algebra.Module((*>))--} import qualified Algebra.Transcendental        as Trans 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 qualified Number.NonNegative       as NonNeg- import NumericPrelude.Numeric---- import qualified Prelude as P import NumericPrelude.Base
+ src/Synthesizer/Generic/Permutation.hs view
@@ -0,0 +1,151 @@+{- |+Permutations of signals as needed for Fast Fourier transforms.+Most functions are independent of the Signal framework.+We could move them as well to Synthesizer.Basic.+-}+module Synthesizer.Generic.Permutation where++import qualified Synthesizer.Basic.NumberTheory as NumberTheory++import qualified Synthesizer.Generic.Signal as SigG+import qualified Synthesizer.State.Signal as SigS++import qualified Data.StorableVector.ST.Strict as SVST+import qualified Data.StorableVector as SV++import qualified Algebra.PrincipalIdealDomain as PID++++type T = SV.Vector Int++apply ::+   (SigG.Transform sig y) =>+   T -> sig y -> sig y+apply p xs =+   SigG.takeStateMatch xs $+   SigS.map (SigG.index xs) $+   SigS.fromStrictStorableSignal p+++size :: T -> Int+size = SV.length+++{- |+> inverse (transposition n m) = transposition m n+-}+transposition ::+   Int -> Int -> T+transposition n m =+   fst $ SV.unfoldrN (n*m)+      (\(i,j,k0) -> Just (i,+         case pred k0 of+            0  -> let j1 = j+1 in (j1, j1, m)+            k1 -> (i+n, j, k1)))+      (0,0,m)+++{-+In general the inverse of a skewGrid+does not look like even a generalized skewGrid.+E.g. @inverse $ skewGrid 3 4@.+-}+skewGrid ::+   Int -> Int -> T+skewGrid n m =+   let len = n*m+   in  fst $ SV.unfoldrN len+          (\(i0,k0) -> Just (i0,+             let k1 = pred k0+                 i1 = i0+n+             in  if k1==0+                   then (mod (i1+m) len, m)+                   else (mod i1 len, k1)))+          (0,m)++{- |+> inverse (skewGrid n m) == skewGridInv n m++In general the inverse of a skewGrid+cannot be expressed like skewGrid or skewGridCRT.+E.g. @inverse $ skewGrid 3 4@.+-}+skewGridInv ::+   Int -> Int -> T+skewGridInv n m =+   SV.pack $+   map (\k ->+      let Just (i,j) = PID.diophantine k n m+      in  mod i m + mod j n * m) $+   take (n*m) $ iterate (1+) 0++skewGridCRT ::+   Int -> Int -> T+skewGridCRT n m =+   let len = n*m+       (ni,mi) = snd $ PID.extendedGCD n m+   in  fst $ SV.unfoldrN len+          (\(i0,k0) -> Just (i0,+             let k1 = pred k0+                 i1 = i0+ni*n+             in  if k1==0+                   then (mod (i1+mi*m) len, m)+                   else (mod i1 len, k1)))+          (0,m)++skewGridCRTInv ::+   Int -> Int -> T+skewGridCRTInv n m =+   fst $ SV.packN (n*m) $+   map (\k -> mod k m + mod k n * m) $+   iterate (1+) 0+++{- |+Beware of 0-based indices stored in the result vector.+-}+multiplicative :: Int -> T+multiplicative ni =+   let n = fromIntegral ni+       gen = NumberTheory.multiplicativeGenerator n+   in  {-+       Since 'gen' is usually 2 or 3,+       the error should occur really only for huge signals.+       -}+       if gen * n > fromIntegral (maxBound :: Int)+         then error "signal too long for Int indexing"+         else fst $ SV.unfoldrN (ni-1)+                 (\x -> Just (x-1, mod (fromInteger gen * x) ni)) 1++{- |+We only need to compute the inverse permutation explicitly,+because not all signal structures support write to arbitrary indices,+thus Generic.Write does not support it.+For strict StorableVector it would be more efficient+to build the vector directly.++It holds:++> inverse . inverse == id+-}+inverse :: T -> T+inverse perm =+   SVST.runSTVector+      (do inv <- SVST.new_ (SV.length perm)+          SigS.sequence_ $+             SigS.zipWith (SVST.write inv)+                (SigS.fromStrictStorableSignal perm)+                (SigS.iterate (1+) 0)+          return inv)++reverse :: T -> T+reverse perm =+   fst $ SV.unfoldrN (SV.length perm)+      (\mn -> Just $+         case mn of+            Nothing -> (SV.head perm, Just $ SV.length perm)+            Just n ->+               let n1 = n-1+               in  (SV.index perm n1, Just n1))+      Nothing
src/Synthesizer/Generic/Piece.hs view
@@ -24,10 +24,7 @@ import qualified Algebra.Transcendental        as Trans 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 qualified Prelude as P import NumericPrelude.Base import NumericPrelude.Numeric
src/Synthesizer/Generic/Signal.hs view
@@ -35,13 +35,21 @@ import qualified Synthesizer.Plain.Signal as Sig import qualified Synthesizer.State.Signal as SigS import qualified Synthesizer.Storable.Signal as SigSt-import qualified Data.StorableVector.Lazy as Vector+import qualified Data.StorableVector.Lazy as SVL+import qualified Data.StorableVector as SV  import qualified Synthesizer.Plain.Modifier as Modifier -import qualified Algebra.NonNegative as NonNeg+import qualified Algebra.ToInteger    as ToInteger+import qualified Algebra.ToRational   as ToRational+import qualified Algebra.Absolute     as Absolute+import qualified Algebra.RealIntegral as RealIntegral+import qualified Algebra.IntegralDomain as Integral+import qualified Algebra.NonNegative  as NonNeg+import qualified Algebra.ZeroTestable as ZeroTestable  import qualified Algebra.Module   as Module+import qualified Algebra.Ring     as Ring import qualified Algebra.Additive as Additive import qualified Algebra.Monoid   as Monoid @@ -59,14 +67,16 @@ import qualified Data.List.Stream as List import Data.Tuple.HT (mapPair, mapFst, ) +import qualified Test.QuickCheck as QC+ -- import NumericPrelude.Numeric import qualified Prelude as P import Prelude-   (Bool, Int, Maybe(Just), maybe, snd,-    (==), (<), (>), (<=), (>=),+   (Bool, Int, Maybe(Just), maybe, fst, snd,+    (==), (<), (>), (<=), (>=), compare, Ordering(..),     flip, uncurry, const, (.), ($), (&&), id, (++),     fmap, return, error, show,-    Eq, Ord, Show, max, min, )+    Eq, Ord, Show, min, )   class Cut.Read (sig y) => Read sig y where@@ -108,7 +118,10 @@ but we need control over packet size in applications with feedback. -} newtype LazySize = LazySize Int-   deriving (Eq, Ord, Show, Additive.C)+   deriving (Eq, Ord, Show,+             Additive.C, Ring.C, ZeroTestable.C,+             ToInteger.C, ToRational.C, Absolute.C,+             RealIntegral.C, Integral.C)  instance Monoid.C LazySize where    idt = LazySize 0@@ -117,6 +130,12 @@ instance NonNeg.C LazySize where    split = NonNeg.splitDefault (\(LazySize n) -> n) LazySize +instance QC.Arbitrary LazySize where+   arbitrary =+      case defaultLazySize of+         LazySize n -> fmap LazySize (QC.choose (1, 2 P.* n))++ {- | This can be used for internal signals that have no observable effect on laziness.@@ -126,7 +145,7 @@ -} defaultLazySize :: LazySize defaultLazySize =-   let (Vector.ChunkSize size) = Vector.defaultChunkSize+   let (SVL.ChunkSize size) = SVL.defaultChunkSize    in  LazySize size  {- |@@ -148,65 +167,108 @@   -- instance Storable y => Read SigSt.T y where-instance Storable y => Read Vector.Vector y where+instance Storable y => Read SVL.Vector y where    {-# INLINE toList #-}-   toList = Vector.unpack+   toList = SVL.unpack    {-# INLINE toState #-}    toState = SigS.fromStorableSignal    {-# INLINE foldL #-}-   foldL = Vector.foldl+   foldL = SVL.foldl    {-# INLINE foldR #-}-   foldR = Vector.foldr+   foldR = SVL.foldr    {-# INLINE index #-}-   index = Vector.index+   index = SVL.index -instance Storable y => Transform Vector.Vector y where+instance Storable y => Transform SVL.Vector y where    {-# INLINE cons #-}-   cons = Vector.cons+   cons = SVL.cons    {-# INLINE takeWhile #-}-   takeWhile = Vector.takeWhile+   takeWhile = SVL.takeWhile    {-# INLINE dropWhile #-}-   dropWhile = Vector.dropWhile+   dropWhile = SVL.dropWhile    {-# INLINE span #-}-   span = Vector.span+   span = SVL.span     {-# INLINE viewL #-}-   viewL = Vector.viewL+   viewL = SVL.viewL    {-# INLINE viewR #-}-   viewR = Vector.viewR+   viewR = SVL.viewR     {-# INLINE map #-}-   map = Vector.map+   map = SVL.map    {-# INLINE scanL #-}-   scanL = Vector.scanl+   scanL = SVL.scanl    {-# INLINE crochetL #-}-   crochetL = Vector.crochetL+   crochetL = SVL.crochetL    {-# INLINE zipWithAppend #-}    zipWithAppend = SigSt.zipWithAppend    withStorableContext ::-   (Vector.ChunkSize -> a) -> (LazySize -> a)+   (SVL.ChunkSize -> a) -> (LazySize -> a) withStorableContext f =-   \(LazySize size) -> f (Vector.ChunkSize size)+   \(LazySize size) -> f (SVL.ChunkSize size) -instance Storable y => Write Vector.Vector y where+instance Storable y => Write SVL.Vector y where    {-# INLINE fromList #-}-   fromList = withStorableContext $ \size -> Vector.pack size+   fromList = withStorableContext $ \size -> SVL.pack size    {-# INLINE repeat #-}-   repeat = withStorableContext $ \size -> Vector.repeat size+   repeat = withStorableContext $ \size -> SVL.repeat size    {-# INLINE replicate #-}-   replicate = withStorableContext $ \size -> Vector.replicate size+   replicate = withStorableContext $ \size -> SVL.replicate size    {-# INLINE iterate #-}-   iterate = withStorableContext $ \size -> Vector.iterate size+   iterate = withStorableContext $ \size -> SVL.iterate size    {-# INLINE unfoldR #-}-   unfoldR = withStorableContext $ \size -> Vector.unfoldr size+   unfoldR = withStorableContext $ \size -> SVL.unfoldr size    {-# INLINE iterateAssociative #-}-   iterateAssociative = withStorableContext $ \size op x -> Vector.iterate size (op x) x -- should be optimized+   iterateAssociative = withStorableContext $ \size op x -> SVL.iterate size (op x) x -- should be optimized   +instance Storable y => Read SV.Vector y where+   {-# INLINE toList #-}+   toList = SV.unpack+   {-# INLINE toState #-}+   toState = SigS.fromStrictStorableSignal+   {-# INLINE foldL #-}+   foldL = SV.foldl+   {-# INLINE foldR #-}+   foldR = SV.foldr+   {-# INLINE index #-}+   index = SV.index++instance Storable y => Transform SV.Vector y where+   {-# INLINE cons #-}+   cons = SV.cons+   {-# INLINE takeWhile #-}+   takeWhile = SV.takeWhile+   {-# INLINE dropWhile #-}+   dropWhile = SV.dropWhile+   {-# INLINE span #-}+   span = SV.span++   {-# INLINE viewL #-}+   viewL = SV.viewL+   {-# INLINE viewR #-}+   viewR = SV.viewR++   {-# INLINE map #-}+   map = SV.map+   {-# INLINE scanL #-}+   scanL = SV.scanl+   {-# INLINE crochetL #-}+   crochetL f acc =+      fst . SVL.crochetLChunk f acc+      -- fst . SV.crochetContL f acc+   {-# INLINE zipWithAppend #-}+   zipWithAppend f xs ys =+      case compare (SV.length xs) (SV.length ys) of+         EQ -> SV.zipWith f xs ys+         LT -> SV.append (SV.zipWith f xs ys) (SV.drop (SV.length xs) ys)+         GT -> SV.append (SV.zipWith f xs ys) (SV.drop (SV.length ys) xs)++ instance Read [] y where    {-# INLINE toList #-}    toList = id@@ -478,6 +540,10 @@    SigS.runSwitchL . toState  +{-# INLINE singleton #-}+singleton :: (Transform sig y) => y -> sig y+singleton x = cons x mempty+ {-# INLINE mix #-} mix :: (Additive.C y, Transform sig y) =>    sig y -> sig y -> sig y@@ -486,13 +552,13 @@ {-# INLINE zipWith #-} zipWith :: (Read sig a, Transform sig b) =>    (a -> b -> b) -> (sig a -> sig b -> sig b)-zipWith h =-   flip runViewL (\next ->-      crochetL-         (\x0 a0 ->-             do (y0,a1) <- next a0-                Just (h y0 x0, a1)))+zipWith h = zipWithState h . toState +{-# INLINE zipWith3 #-}+zipWith3 :: (Read sig a, Read sig b, Transform sig c) =>+   (a -> b -> c -> c) -> (sig a -> sig b -> sig c -> sig c)+zipWith3 h as bs = zipWithState3 h (toState as) (toState bs)+ {-# INLINE zipWithState #-} zipWithState :: (Transform sig b) =>    (a -> b -> b) -> SigS.T a -> sig b -> sig b@@ -508,6 +574,18 @@ zipWithState3 h a b =    zipWithState ($) (SigS.zipWith h a b) +{- |+@takeStateMatch len xs@+keeps a prefix of @xs@ of the same length and block structure as @len@+and stores it in the same type of container as @len@.+-}+{-# INLINE takeStateMatch #-}+takeStateMatch :: (Transform sig a) =>+   sig a -> SigS.T a -> sig a+takeStateMatch x y =+   zipWithState const y x++ {-# INLINE delay #-} delay :: (Write sig y) =>    LazySize -> y -> Int -> sig y -> sig y@@ -534,7 +612,7 @@                  It's absolutely necessary that this function preserves the chunk structure                  and that it does not look a chunk ahead.                  That's guaranteed for processes that do not look ahead at all,-                 like 'Vector.map', 'Vector.crochetL' and+                 like 'SVL.map', 'SVL.crochetL' and                  all of type @Causal.Process@. -}    -> sig y -- ^ input    -> sig y -- ^ output has the same length as the input@@ -547,6 +625,17 @@ {-# INLINE sum #-} sum :: (Additive.C a, Read sig a) => sig a -> a sum = foldL (Additive.+) Additive.zero++{-# INLINE sum1 #-}+sum1 :: (Additive.C a, Read sig a) => sig a -> a+sum1 = SigS.foldL1 (Additive.+) . toState+{-+sum1 :: (Additive.C a, Transform sig a) => sig a -> a+sum1 =+   switchL+      (error "Generic.Signal.sum1: signal must be non-empty in order to avoid to use a non-existing zero")+      (foldL (Additive.+))+-}  {-# INLINE monoidConcatMap #-} monoidConcatMap :: (Read sig a, Monoid m) => (a -> m) -> sig a -> m
src/Synthesizer/Generic/Signal2.hs view
@@ -36,8 +36,8 @@ import Data.Tuple.HT (fst3, snd3, thd3, ) import Prelude    (Integral,-    Bool, Int, Maybe(Just), maybe, fst, snd,-    flip, ($), (.),+    Maybe(Just), maybe, fst, snd,+    flip, const, ($), (.),     return, )  @@ -205,3 +205,14 @@    (a -> b -> c -> d) -> (SigS.T a -> SigS.T b -> sig c -> sig d) zipWithState3 h a b =    zipWithState ($) (SigS.zipWith h a b)++{- |+@takeStateMatch len xs@+keeps a prefix of @xs@ of the same length and block structure as @len@+and stores it in the same type of container as @len@.+-}+{-# INLINE takeStateMatch #-}+takeStateMatch :: (Transform sig a b) =>+   sig a -> SigS.T b -> sig b+takeStateMatch x y =+   zipWithState const y x
src/Synthesizer/Interpolation/Class.hs view
@@ -10,7 +10,6 @@ import qualified Algebra.Module as Module import qualified Algebra.PrincipalIdealDomain as PID import qualified Algebra.Ring as Ring-import qualified Algebra.Additive as Additive  import qualified Sound.Frame.NumericPrelude.Stereo as Stereo import qualified Number.Ratio as Ratio
src/Synthesizer/Interpolation/Custom.hs view
@@ -24,8 +24,6 @@    )  import qualified Algebra.Field     as Field-import qualified Algebra.Ring      as Ring-import qualified Algebra.Additive  as Additive  import Synthesizer.Interpolation.Class ((+.*), )
src/Synthesizer/Interpolation/Module.hs view
@@ -24,12 +24,7 @@    )  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((*>))  import Control.Applicative (liftA2, ) import Synthesizer.ApplicativeUtility (liftA4, )
src/Synthesizer/Plain/Analysis.hs view
@@ -25,12 +25,9 @@  import qualified Data.IntMap as IntMap --- import Algebra.Module((*>))- import Data.Array (accumArray) import Data.List (foldl', ) -import qualified Prelude as P import NumericPrelude.Base import NumericPrelude.Numeric @@ -323,10 +320,9 @@ chirpTransform :: Ring.C y =>    y -> Sig.T y -> Sig.T y chirpTransform z xs =-   let powers = Ctrl.curveMultiscaleNeutral (*) z one-       powerPowers =-          map (\zn -> Ctrl.curveMultiscaleNeutral (*) zn one) powers-   in  map (scalarProduct xs) powerPowers+   map (scalarProduct xs) $+   map (\zn -> Ctrl.curveMultiscaleNeutral (*) zn one) $+   Ctrl.curveMultiscaleNeutral (*) z one   binarySign ::@@ -335,11 +331,33 @@    map (binaryLevelFromBool . (zero <=))  {- |-The output type could be different from the input type-but then we would need a conversion from output to input for feedback.+A kind of discretization for signals with sample values between -1 and 1.+If you smooth the resulting signal+(after you transformed with 'map binaryLevelToNumber'),+you should obtain an approximation to the input signal. -} deltaSigmaModulation ::    RealRing.C y => Sig.T y -> Sig.T BinaryLevel deltaSigmaModulation x =+   let y = binarySign (Integration.run (x - (zero : map binaryLevelToNumber y)))+   in  y+{-    let y = binarySign (Integration.runInit zero (x - map binaryLevelToNumber y))+   in  y+-}++{- |+A kind of discretization for signals with sample values between 0 and a threshold.+We accumulate input values and emit a threshold value+whenever the accumulator exceeds the threshold.+This is intended for generating clicks from input noise.++See also 'deltaSigmaModulation'.+-}+deltaSigmaModulationPositive ::+   RealRing.C y => y -> Sig.T y -> Sig.T y+deltaSigmaModulationPositive threshold x =+   let y =+          map (\xi -> if xi>=threshold then threshold else zero) $+          Integration.run (x - (zero:y))    in  y
src/Synthesizer/Plain/Builder.hs view
@@ -10,8 +10,7 @@ import qualified Algebra.ToInteger as ToInteger import qualified Algebra.RealField as RealField -import qualified Prelude as P98-+import Prelude () import NumericPrelude.Base import NumericPrelude.Numeric
src/Synthesizer/Plain/Control.hs view
@@ -14,14 +14,10 @@ import qualified Algebra.Ring                  as Ring import qualified Algebra.Additive              as Additive -import Algebra.Module((*>))--import Number.Complex (cis,real)--- import qualified Number.Complex as Complex+import Number.Complex (cis,real, ) import Data.List (zipWith4, tails, ) import Data.List.HT (iterateAssociative, ) -import qualified Prelude as P import NumericPrelude.Base import NumericPrelude.Numeric
src/Synthesizer/Plain/Effect.hs view
@@ -15,16 +15,18 @@ import Synthesizer.Plain.Effect.Glass(glass)  import qualified Synthesizer.Plain.File as File+import qualified Control.Monad.Exception.Synchronous as Exc import System.Exit(ExitCode) import System.Cmd(rawSystem)  main :: IO ExitCode main =    let rate = 44100-   in  do {- File.writeMono "test" rate+   in  Exc.toExitCodeT $+       do {- File.writeMono "test" rate                 (take (round (3*rate)) (soundD rate)) -}-          File.renderMonoToInt16 "test.aiff" rate soundE-          rawSystem "play" ["test.aiff"]+          Exc.fromExitCodeT $ File.renderMonoToInt16 "test.aiff" rate soundE+          Exc.fromExitCodeT $ rawSystem "play" ["test.aiff"]   soundE, soundB, soundA,
src/Synthesizer/Plain/Effect/Fly.hs view
@@ -1,15 +1,17 @@ {-# LANGUAGE NoImplicitPrelude #-} module Synthesizer.Plain.Effect.Fly where +import qualified Synthesizer.Causal.Spatial as Spatial+import qualified Synthesizer.Causal.Process as Causal+ import qualified Synthesizer.Plain.Oscillator as Osci import qualified Synthesizer.Plain.Filter.NonRecursive as FiltNR import qualified Synthesizer.Plain.Interpolation as Interpolation-import qualified Synthesizer.Plain.Miscellaneous as Syn  import qualified Synthesizer.Plain.File as File import System.Exit(ExitCode) -import System.Random+import System.Random (randomRs, mkStdGen, )  import qualified Algebra.NormedSpace.Euclidean as Euc @@ -18,7 +20,7 @@   {--  ghc -O -fvia-C -fglasgow-exts -fallow-undecidable-instances --make Fly.hs && echo start && time a.out+  ghc -O -fvia-C --make Fly.hs && echo start && time a.out -}  main :: IO ExitCode@@ -29,7 +31,7 @@ sampleRate :: Double sampleRate = 44100 -{-| stereo sound of a humming fly -}+{- | stereo sound of a humming fly -} fly :: [(Double,Double)] fly =    let pinkNoise seed freq range =@@ -49,11 +51,12 @@                (flyCoord 654891)         channel ear =-          let (phase,volumes) = Syn.receive3Dsound 1 0.1 ear trajectory+          let (phase,volumes) =+                 unzip $ Causal.apply (Spatial.receive3Dsound 1 0.1 ear) trajectory               -- (*sampleRate) in 'speed' and               -- (/sampleRate) in 'freqs' neutralizes-              speeds  = map (\v -> 250/sampleRate+2*(Euc.norm v))-                            (zipWith subtract (tail trajectory) trajectory)+              speeds  = map (\v -> 250/sampleRate + 2 * Euc.norm v)+                            (FiltNR.differentiate trajectory)               freqs   = zipWith (+) speeds (FiltNR.differentiate phase)               sound   = Osci.freqModSaw 0 freqs           in  zipWith (*) (map (10*) volumes) sound
src/Synthesizer/Plain/Effect/Glass.hs view
@@ -13,11 +13,12 @@  import qualified Algebra.Transcendental as Trans import qualified Algebra.RealField      as RealField-import qualified Algebra.Additive       as Additive import qualified Algebra.Module         as Module  import System.Random(randomRs, mkStdGen) +import qualified Data.List.HT as ListHT+ import NumericPrelude.Base import NumericPrelude.Numeric as NP @@ -63,8 +64,4 @@    in  diffs (NonNeg.fromNumber 0)  timeDiffs :: [Bool] -> [NonNeg.Int]-timeDiffs = map (NonNeg.fromNumber . length) . segmentBefore id--segmentBefore :: (a -> Bool) -> [a] -> [[a]]-segmentBefore p =-   foldr (\ x ~(y:ys) -> (if p x then ([]:) else id) ((x:y):ys)) [[]]+timeDiffs = map (NonNeg.fromNumber . length) . ListHT.segmentBefore id
src/Synthesizer/Plain/File.hs view
@@ -1,5 +1,22 @@ {-# LANGUAGE NoImplicitPrelude #-}-module Synthesizer.Plain.File where+module Synthesizer.Plain.File (+   render,+   renderToInt16,+   renderMonoToInt16,+   renderStereoToInt16,+   write,+   writeToInt16,+   writeMonoToInt16,+   writeStereoToInt16,+   writeRaw,+   writeRawCompressed,+   rawToAIFF,+   compress,+   readAIFFMono,+   readMonoFromInt16,+   -- will no longer be exported+   getInt16List,+   ) where  import qualified Sound.Sox.Convert as Convert import qualified Sound.Sox.Frame as Frame@@ -18,6 +35,8 @@ import Foreign.Storable (Storable, ) import Data.Int (Int16, ) +import qualified Control.Monad.Exception.Synchronous as Exc+import Control.Monad.Trans.Class (lift, ) import System.Cmd (rawSystem, ) import System.Exit (ExitCode, ) import Control.Monad (liftM2, )@@ -113,8 +132,9 @@ writeRawCompressed :: (RealRing.C a, Frame.C v, Storable v) =>    SoxOpt.T -> FilePath -> a -> [v] -> IO ExitCode writeRawCompressed opts fileName sampleRate signal =-   do writeRaw opts fileName sampleRate signal-      compress fileName+   Exc.toExitCodeT $+   do Exc.fromExitCodeT $ writeRaw opts fileName sampleRate signal+      Exc.fromExitCodeT $ compress fileName   {-# DEPRECATED rawToAIFF "If you want to generate AIFF, then just write to files with .aiff filename extension. If you want to convert files to AIFF, use Sound.Sox.Convert." #-}@@ -133,9 +153,9 @@           SoxOpt.none fileNameAIFF  compress :: FilePath -> IO ExitCode-compress fileName =-   do rawSystem "oggenc" ["--quality", "5", fileName]-      rawSystem "lame"+compress fileName = Exc.toExitCodeT $+   do Exc.fromExitCodeT $ rawSystem "oggenc" ["--quality", "5", fileName]+      Exc.fromExitCodeT $ rawSystem "lame"          ["-h", fileName, FilePath.replaceExtension fileName "mp3"]  @@ -155,10 +175,12 @@ -} readAIFFMono :: (Field.C a) => FilePath -> IO [a] readAIFFMono file =-   do --putStrLn ("sox "++file++" "++tmp)-      let tmp = FilePath.replaceExtension file "sw"-      Convert.simple SoxOpt.none file SoxOpt.none tmp-      fmap (map BinSmp.int16ToCanonical) (FileL.readInt16StreamStrict tmp)+   let tmp = FilePath.replaceExtension file "s16"+   in  Exc.resolveT (const $ return []) $ do+          -- lift $ putStrLn ("sox "++file++" "++tmp)+          Exc.fromExitCodeT $ Convert.simple SoxOpt.none file SoxOpt.none tmp+          fmap (map BinSmp.int16ToCanonical) $+             lift $ FileL.readInt16StreamStrict tmp   {- |@@ -167,12 +189,14 @@ readMonoFromInt16 :: (Field.C a) => FilePath -> IO [a] readMonoFromInt16 fileName =    Read.open SoxOpt.none fileName >>=-   Read.withHandle1 (fmap (Get.runGet getInt16List) . B.hGetContents) >>=+   Read.withHandle1 (fmap (Get.runGet getInt16ListPrivate) . B.hGetContents) >>=    return . map BinSmp.int16ToCanonical -getInt16List :: Get.Get [Int16]-getInt16List =+{-# DEPRECATED getInt16List "This function will no longer be exported" #-}+getInt16List, getInt16ListPrivate :: Get.Get [Int16]+getInt16List = getInt16ListPrivate+getInt16ListPrivate =    do b <- Get.isEmpty       if b         then return []-        else liftM2 (:) (fmap fromIntegral Get.getWord16host) getInt16List+        else liftM2 (:) (fmap fromIntegral Get.getWord16host) getInt16ListPrivate
src/Synthesizer/Plain/Filter/Delay.hs view
@@ -9,7 +9,6 @@  import qualified Algebra.Module    as Module import qualified Algebra.RealField as RealField-import qualified Algebra.Additive  as Additive  import qualified Synthesizer.Plain.Interpolation as Interpolation @@ -17,7 +16,6 @@ import qualified Synthesizer.Plain.Filter.Delay.List  as DelayList import qualified Synthesizer.Plain.Filter.Delay.Block as DelayBlock -import qualified Prelude as P import NumericPrelude.Base import NumericPrelude.Numeric
src/Synthesizer/Plain/Filter/Delay/Block.hs view
@@ -18,7 +18,6 @@  import Test.QuickCheck ((==>), Property) -import qualified Prelude as P import NumericPrelude.Base import NumericPrelude.Numeric
src/Synthesizer/Plain/Filter/Delay/List.hs view
@@ -8,7 +8,6 @@  import Data.List(tails) -import qualified Prelude as P import NumericPrelude.Base import NumericPrelude.Numeric
src/Synthesizer/Plain/Filter/Delay/ST.hs view
@@ -13,7 +13,6 @@ import Control.Monad.ST.Lazy(runST,strictToLazyST,ST) import Data.Array.ST -import qualified Prelude as P import NumericPrelude.Base import NumericPrelude.Numeric
src/Synthesizer/Plain/Filter/LinearPredictive.hs view
@@ -1,8 +1,6 @@ module Synthesizer.Plain.Filter.LinearPredictive where  import qualified Algebra.Field    as Field-import qualified Algebra.Ring     as Ring-import qualified Algebra.Additive as Additive import Synthesizer.Plain.Analysis (scalarProduct)  import qualified Data.List.Match as ListMatch
src/Synthesizer/Plain/Filter/NonRecursive.hs view
@@ -19,7 +19,7 @@ import qualified Algebra.Ring           as Ring import qualified Algebra.Additive       as Additive -import Algebra.Module(linearComb, (*>))+import Algebra.Module(linearComb, )  import Data.Function.HT (nest, ) import Data.Tuple.HT (mapPair, swap, )
src/Synthesizer/Plain/Filter/Recursive.hs view
@@ -12,15 +12,11 @@ module Synthesizer.Plain.Filter.Recursive where  import qualified Algebra.Module                as Module--- import qualified Algebra.Transcendental        as Trans--- import qualified Algebra.Field                 as Field--- import qualified Algebra.Ring                  as Ring import qualified Algebra.Additive              as Additive  import Algebra.Additive((+), (-), negate, ) import Algebra.Module((*>)) -import qualified Prelude as P import NumericPrelude.Base import NumericPrelude.Numeric
src/Synthesizer/Plain/Filter/Recursive/Allpass.hs view
@@ -28,9 +28,6 @@ import qualified Algebra.Transcendental        as Trans import qualified Algebra.Field                 as Field import qualified Algebra.Ring                  as Ring-import qualified Algebra.Additive              as Additive--import Algebra.Module((*>))  import qualified Number.Complex as Complex import Data.Tuple.HT (mapSnd, )
src/Synthesizer/Plain/Filter/Recursive/AllpassPoly.hs view
@@ -18,10 +18,6 @@ import qualified Algebra.RealTranscendental    as RealTrans import qualified Algebra.Transcendental        as Trans import qualified Algebra.Field                 as Field-import qualified Algebra.Ring                  as Ring-import qualified Algebra.Additive              as Additive--import Algebra.Module((*>))  import Number.Complex (cis,(+:),real,imag) import qualified Number.Complex as Complex
src/Synthesizer/Plain/Filter/Recursive/Butterworth.hs view
@@ -24,14 +24,11 @@  import qualified Algebra.Module                as Module import qualified Algebra.Transcendental        as Trans-import qualified Algebra.Field                 as Field import qualified Algebra.Ring                  as Ring-import qualified Algebra.Additive              as Additive  import qualified Data.StorableVector as SV import Foreign.Storable (Storable) -import qualified Prelude as P import NumericPrelude.Base import NumericPrelude.Numeric
src/Synthesizer/Plain/Filter/Recursive/Chebyshev.hs view
@@ -22,24 +22,16 @@ import qualified Synthesizer.Causal.Process as Causal import Control.Arrow ((>>>), (^>>), (&&&), ) --- import qualified Algebra.VectorSpace           as VectorSpace import qualified Algebra.Module                as Module import qualified Algebra.Transcendental        as Trans-import qualified Algebra.Field                 as Field import qualified Algebra.Ring                  as Ring-import qualified Algebra.Additive              as Additive -import Algebra.Module((*>))- import Number.Complex (real, imag, cis, ) import qualified Number.Complex as Complex --- import Control.Monad.Trans.State (State(..), evalState)- import qualified Data.StorableVector as SV import Foreign.Storable (Storable) -import qualified Prelude as P import NumericPrelude.Base import NumericPrelude.Numeric @@ -272,6 +264,14 @@               Filt2.adjustPassband kind                  (flip (partialParameterA order ratio) c) freq) $            circleVec)++{-# INLINE canonicalizeParameterA #-}+canonicalizeParameterA ::+   (Ring.C a, Storable a) =>+   ParameterA a -> Cascade.Parameter a+canonicalizeParameterA (amp, Cascade.Parameter p) =+   Cascade.Parameter+      (SV.switchL SV.empty (\h -> SV.cons (Filt2.amplify amp h)) p)   type ParameterB a = Cascade.Parameter a
src/Synthesizer/Plain/Filter/Recursive/Comb.hs view
@@ -21,13 +21,9 @@ import qualified Synthesizer.Plain.Control as Ctrl  import qualified Algebra.Module                as Module--- import qualified Algebra.Field                 as Field import qualified Algebra.Ring                  as Ring import qualified Algebra.Additive              as Additive -import Algebra.Module((*>))--import qualified Prelude as P import NumericPrelude.Base import NumericPrelude.Numeric
src/Synthesizer/Plain/Filter/Recursive/FirstOrder.hs view
@@ -2,7 +2,7 @@ {-# LANGUAGE MultiParamTypeClasses #-} {-# LANGUAGE FlexibleInstances #-} {- |-Copyright   :  (c) Henning Thielemann 2008+Copyright   :  (c) Henning Thielemann 2008-2011 License     :  GPL  Maintainer  :  synthesizer@henning-thielemann.de@@ -211,20 +211,26 @@   {-# INLINE step #-}-step :: (Ring.C a, Module.C a v) =>+step :: (Module.C a v) =>    Parameter a -> v -> State v (Result v) step c x =    fmap (\lp -> Result (x-lp) lp) (lowpassStep c x)  {-# INLINE modifierInit #-}-modifierInit :: (Ring.C a, Module.C a v) =>-   Modifier.Initialized v v (Parameter a) v v+modifierInit :: (Module.C a v) =>+   Modifier.Initialized v v (Parameter a) v (Result v) modifierInit =-   Modifier.Initialized id lowpassStep+   Modifier.Initialized id step  {-# INLINE modifier #-}-modifier :: (Ring.C a, Module.C a v) =>-   Modifier.Simple v (Parameter a) v v+modifier :: (Module.C a v) =>+   Modifier.Simple v (Parameter a) v (Result v) modifier =    Sig.modifierInitialize modifierInit zero +{-# INLINE causal #-}+causal ::+   (Module.C a v) =>+   Causal.T (Parameter a, v) (Result v)+causal =+   Causal.fromSimpleModifier modifier
src/Synthesizer/Plain/Filter/Recursive/FirstOrderComplex.hs view
@@ -38,15 +38,10 @@ import qualified Algebra.Module                as Module import qualified Algebra.Transcendental        as Trans import qualified Algebra.Algebraic             as Algebraic-import qualified Algebra.Field                 as Field import qualified Algebra.Ring                  as Ring-import qualified Algebra.Additive              as Additive --- import Algebra.Module((*>))- import Control.Monad.Trans.State (State, state, ) -import qualified Prelude as P import NumericPrelude.Base import NumericPrelude.Numeric
src/Synthesizer/Plain/Filter/Recursive/Hilbert.hs view
@@ -39,9 +39,6 @@ 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((*>))  import qualified Number.Complex as Complex import Number.Complex ((+:), )
src/Synthesizer/Plain/Filter/Recursive/Moog.hs view
@@ -28,11 +28,7 @@  import qualified Algebra.Module                as Module import qualified Algebra.Transcendental        as Trans-import qualified Algebra.Field                 as Field import qualified Algebra.Ring                  as Ring-import qualified Algebra.Additive              as Additive--import Algebra.Module((*>))  import Data.Function.HT (nest, )
src/Synthesizer/Plain/Filter/Recursive/SecondOrder.hs view
@@ -29,14 +29,10 @@  import qualified Synthesizer.Causal.Process as Causal --- import qualified Algebra.VectorSpace           as VectorSpace import qualified Algebra.Module                as Module--- import qualified Algebra.Transcendental        as Trans import qualified Algebra.Field                 as Field import qualified Algebra.Ring                  as Ring import qualified Algebra.Additive              as Additive--import Algebra.Module((*>))  import Data.List (zipWith6)
src/Synthesizer/Plain/Filter/Recursive/SecondOrderCascade.hs view
@@ -27,21 +27,14 @@  import qualified Synthesizer.Causal.Process as Causal --- import qualified Algebra.VectorSpace           as VectorSpace import qualified Algebra.Module                as Module--- import qualified Algebra.Transcendental        as Trans--- import qualified Algebra.Field                 as Field import qualified Algebra.Ring                  as Ring--- import qualified Algebra.Additive              as Additive --- import Algebra.Module((*>))- import qualified Control.Monad.Trans.State as MS  import qualified Data.StorableVector as SV import Foreign.Storable (Storable(..)) -import qualified Prelude as P import NumericPrelude.Base import NumericPrelude.Numeric
src/Synthesizer/Plain/Filter/Recursive/Universal.hs view
@@ -31,7 +31,6 @@  import qualified Algebra.Module                as Module import qualified Algebra.Transcendental        as Trans-import qualified Algebra.Field                 as Field import qualified Algebra.Ring                  as Ring import qualified Algebra.Additive              as Additive
src/Synthesizer/Plain/IO.hs view
@@ -23,8 +23,6 @@ import qualified Data.ByteString.Lazy as B import qualified Data.Binary.Builder as Builder -import qualified Algebra.Ring      as Ring- import Data.Char (ord, )  import qualified Prelude as P98
src/Synthesizer/Plain/Interpolation.hs view
@@ -33,7 +33,6 @@ import qualified Algebra.Ring      as Ring import qualified Algebra.Additive  as Additive -import Algebra.Additive(zero) import Data.Maybe (fromMaybe) import qualified Data.List.HT as ListHT
src/Synthesizer/Plain/Miscellaneous.hs view
@@ -1,25 +1,26 @@ {-# LANGUAGE NoImplicitPrelude #-}-module Synthesizer.Plain.Miscellaneous where+module Synthesizer.Plain.Miscellaneous+   {-# DEPRECATED "Use Synthesizer.Causal.Spatial instead" #-} where +import qualified Synthesizer.Causal.Spatial as Spatial+import qualified Synthesizer.Causal.Process as Causal+ import qualified Algebra.NormedSpace.Euclidean as Euc import qualified Algebra.Field                 as Field--- import qualified Algebra.Ring                  as Ring--- import qualified Algebra.Additive              as Additive -import qualified Prelude as P import NumericPrelude.Base-import NumericPrelude.Numeric+-- import NumericPrelude.Numeric   {- * Spatial effects -} -{-| simulate an moving sounding object-   convert the way of the object through 3D space-   into a delay and attenuation information,-   sonicDelay is the reciprocal of the sonic velocity -}+{-|+simulate an moving sounding object++convert the way of the object through 3D space+into a delay and attenuation information,+sonicDelay is the reciprocal of the sonic velocity+-} receive3Dsound :: (Field.C a, Euc.C a v) => a -> a -> v -> [v] -> ([a],[a])-receive3Dsound att sonicDelay ear way =-   let dists   = map (Euc.norm) (map (subtract ear) way)-       phase   = map (sonicDelay*) dists-       volumes = map (\x -> 1/(att+x)^2) dists-   in  (phase, volumes)+receive3Dsound att sonicDelay ear =+   unzip . Causal.apply (Spatial.receive3Dsound att sonicDelay ear)
src/Synthesizer/Plain/Noise.hs view
@@ -11,7 +11,6 @@  import Data.List.HT (sliceVertical, ) -import qualified Prelude as P import NumericPrelude.Base import NumericPrelude.Numeric
src/Synthesizer/Plain/Oscillator.hs view
@@ -23,26 +23,13 @@  import Synthesizer.Plain.ToneModulation (freqsToPhases, ) -{--import qualified Algebra.RealTranscendental    as RealTrans-import qualified Algebra.Module                as Module-import qualified Algebra.VectorSpace           as VectorSpace--import Algebra.Module((*>))--} import qualified Algebra.Transcendental        as Trans import qualified Algebra.RealField             as RealField import qualified Algebra.RealRing              as RealRing-import qualified Algebra.Ring                  as Ring-import qualified Algebra.Additive              as Additive --- import qualified Number.NonNegative       as NonNeg- import Data.Tuple.HT (mapFst, mapSnd, )  import NumericPrelude.Numeric---- import qualified Prelude as P import NumericPrelude.Base
src/Synthesizer/Plain/Signal.hs view
@@ -1,7 +1,7 @@ {-# OPTIONS_GHC -fglasgow-exts #-} {- glasgow-exts are for the rules -} {- |-Copyright   :  (c) Henning Thielemann 2008+Copyright   :  (c) Henning Thielemann 2008-2011 License     :  GPL  Maintainer  :  synthesizer@henning-thielemann.de@@ -15,6 +15,7 @@ import qualified Synthesizer.Plain.Modifier as Modifier  import qualified Data.List.Match as ListMatch+import qualified Data.List.HT    as ListHT import qualified Data.List       as List  import Data.Tuple.HT (forcePair, mapFst, mapSnd, )@@ -146,8 +147,9 @@ dropMarginRem :: Int -> Int -> T a -> (Int, T a) dropMarginRem n m =    head .-   dropMargin n m .-   zipWithTails (,) (iterate pred m)+   dropMargin (1+n) m .+   zip (iterate (max 0 . pred) m) .+   ListHT.tails  dropMargin :: Int -> Int -> T a -> T a dropMargin n m xs =@@ -170,7 +172,7 @@ zipWithTails ::    (y0 -> T y1 -> y2) -> T y0 -> T y1 -> T y2 zipWithTails f xs =-   zipWith f xs . init . List.tails+   zipWith f xs . init . ListHT.tails  zipWithRest ::    (y0 -> y0 -> y1) ->
src/Synthesizer/Plain/ToneModulation.hs view
@@ -54,12 +54,9 @@ import qualified Synthesizer.Plain.Signal as Sig import qualified Synthesizer.Plain.Interpolation as Interpolation import Synthesizer.Interpolation (Margin, )--- import qualified Data.Array as Array import Data.Array (Array, (!), listArray, ) --- import qualified Algebra.Transcendental        as Trans import qualified Algebra.RealField             as RealField-import qualified Algebra.Field                 as Field import qualified Algebra.RealRing              as RealRing import qualified Algebra.Ring                  as Ring import qualified Algebra.Additive              as Additive
src/Synthesizer/Plain/Wave.hs view
@@ -10,11 +10,7 @@ import qualified Synthesizer.Plain.Signal as Sig import Data.Array ((!), listArray) --- import qualified Synthesizer.Basic.Phase as Phase- import qualified Algebra.RealField             as RealField-import qualified Algebra.Ring                  as Ring-import qualified Algebra.Additive              as Additive  import NumericPrelude.Numeric
src/Synthesizer/State/Analysis.hs view
@@ -20,16 +20,11 @@ import qualified Algebra.NormedSpace.Euclidean as NormedEuc import qualified Algebra.NormedSpace.Sum       as NormedSum -import qualified Data.Array as Array- import qualified Data.IntMap as IntMap---- import Algebra.Module((*>))+import qualified Data.Array as Array  import Data.Array (accumArray)--- import Data.List (foldl', ) -import qualified Prelude as P import NumericPrelude.Base import NumericPrelude.Numeric @@ -370,7 +365,6 @@ chirpTransform :: Ring.C y =>    y -> Sig.T y -> Sig.T y chirpTransform z xs =-   let powers = Ctrl.curveMultiscaleNeutral (*) z one-       powerPowers =-          Sig.map (\zn -> Ctrl.curveMultiscaleNeutral (*) zn one) powers-   in  Sig.map (scalarProduct xs) powerPowers+   Sig.map (scalarProduct xs) $+   Sig.map (\zn -> Ctrl.curveMultiscaleNeutral (*) zn one) $+   Ctrl.curveMultiscaleNeutral (*) z one
src/Synthesizer/State/Control.hs view
@@ -24,14 +24,8 @@ import qualified Algebra.Ring                  as Ring import qualified Algebra.Additive              as Additive -import Algebra.Module((*>))---- import Number.Complex (cis,real)--- import qualified Number.Complex as Complex- import Data.Ix (Ix, ) -import qualified Prelude as P import NumericPrelude.Base import NumericPrelude.Numeric
src/Synthesizer/State/Cut.hs view
@@ -28,7 +28,7 @@ import qualified Algebra.Additive as Additive  import qualified Data.Array as Array-import Data.Array (Array, Ix, (!), elems, )+import Data.Array (Array, Ix, (!), ) import Control.Applicative (Applicative, ) import Data.Traversable (sequenceA, )
src/Synthesizer/State/Displacement.hs view
@@ -3,13 +3,10 @@  import qualified Synthesizer.State.Signal as Sig --- import qualified Algebra.Module                as Module import qualified Algebra.Transcendental        as Trans--- import qualified Algebra.Field                 as Field import qualified Algebra.Ring                  as Ring import qualified Algebra.Additive              as Additive -import qualified Prelude as P import NumericPrelude.Base import NumericPrelude.Numeric
src/Synthesizer/State/Filter/NonRecursive.hs view
@@ -21,8 +21,6 @@ import qualified Algebra.Ring           as Ring import qualified Algebra.Additive       as Additive -import Algebra.Module( {- linearComb, -} (*>))- import Data.Function.HT (nest, ) import Data.Tuple.HT (mapFst, )
src/Synthesizer/State/Filter/Recursive/Comb.hs view
@@ -20,13 +20,9 @@ import qualified Synthesizer.State.Filter.Delay as Delay  import qualified Algebra.Module                as Module--- import qualified Algebra.Field                 as Field import qualified Algebra.Ring                  as Ring import qualified Algebra.Additive              as Additive -import Algebra.Module((*>))--import qualified Prelude as P import NumericPrelude.Base import NumericPrelude.Numeric
src/Synthesizer/State/Interpolation.hs view
@@ -1,17 +1,12 @@ {-# LANGUAGE NoImplicitPrelude #-} module Synthesizer.State.Interpolation where -import qualified Synthesizer.Interpolation as Interpolation import Synthesizer.Interpolation    (T, offset, number, func, )  import qualified Synthesizer.State.Signal  as Sig --- import qualified Algebra.Module    as Module--- import qualified Algebra.Field     as Field import qualified Algebra.RealRing  as RealRing-import qualified Algebra.Ring      as Ring-import qualified Algebra.Additive  as Additive  import Data.Maybe (fromMaybe)
src/Synthesizer/State/Noise.hs view
@@ -10,7 +10,6 @@ import System.Random (Random, RandomGen, randomR, mkStdGen, ) import qualified System.Random as Rnd -import qualified Prelude as P import NumericPrelude.Base import NumericPrelude.Numeric
src/Synthesizer/State/NoiseCustom.hs view
@@ -16,7 +16,6 @@ import System.Random (Random, RandomGen, ) import qualified System.Random as Rnd -import qualified Prelude as P import NumericPrelude.Base import NumericPrelude.Numeric
src/Synthesizer/State/Signal.hs view
@@ -81,7 +81,7 @@    format = showsPrec  instance Functor T where-   fmap = map+   fmap g (Cons f s) = Cons (fmap g f) s  instance App.Applicative T where    pure = singleton@@ -157,6 +157,13 @@    generate PtrSt.viewL .    PtrSt.cons +{-# INLINE fromStrictStorableSignal #-}+fromStrictStorableSignal ::+   (Storable a) =>+   V.Vector a -> T a+fromStrictStorableSignal xs =+   map (V.index xs) $ take (V.length xs) $ iterate succ zero+ {-# INLINE toStorableSignal #-} toStorableSignal ::    (Storable a) =>@@ -233,7 +240,8 @@  {-# INLINE map #-} map :: (a -> b) -> (T a -> T b)-map f = crochetL (\x _ -> Just (f x, ())) ()+map = fmap+-- map f = crochetL (\x _ -> Just (f x, ())) ()   {- |@@ -266,7 +274,9 @@  {-# INLINE take #-} take :: Int -> T a -> T a-take = crochetL (\x n -> toMaybe (n>zero) (x, pred n))+take n =+   map snd . takeWhile ((>0) . fst) . zip (iterate pred n)+   -- crochetL (\x n -> toMaybe (n>zero) (x, pred n))  {-# INLINE takeWhile #-} takeWhile :: (a -> Bool) -> T a -> T a@@ -340,6 +350,13 @@ foldL :: (acc -> x -> acc) -> acc -> T x -> acc foldL f = foldL' (flip f) +{-# INLINE foldL1 #-}+foldL1 :: (x -> x -> x) -> T x -> x+foldL1 f =+   switchL+      (error "State.Signal.foldL1: empty signal")+      (foldL f)+ {-# INLINE length #-} length :: T a -> Int length = foldL' (const succ) zero@@ -486,8 +503,8 @@ dropMarginRem :: Int -> Int -> T a -> (Int, T a) dropMarginRem n m =    switchL (error $ "StateSignal.dropMaringRem: length xs < " ++ show n) const .-   dropMargin n m .-   zipWithTails (,) (iterate pred m)+   dropMargin (succ n) m .+   zipWithTails1 (,) (iterate (max 0 . pred) m)  {-# INLINE dropMargin #-} dropMargin :: Int -> Int -> T a -> T a@@ -607,9 +624,9 @@ append :: T a -> T a -> T a append xs ys =    generate-      (\(b,xs0) ->+      (\(b,xys) ->           mplus-             (fmap (mapSnd ((,) b)) $ viewL xs0)+             (fmap (mapSnd ((,) b)) $ viewL xys)              (if b                 then Nothing                 else fmap (mapSnd ((,) True)) $ viewL ys))@@ -781,7 +798,7 @@             return (f xs, ys))) -} --- only non-empty suffixes are processed+-- | only non-empty suffixes are processed {-# INLINE zipWithTails #-} zipWithTails ::    (y0 -> T y1 -> y2) -> T y0 -> T y1 -> T y2@@ -801,6 +818,26 @@          (xs2,ys2))) -} +-- | in contrast to 'zipWithTails' it also generates the empty suffix (once)+{-# INLINE zipWithTails1 #-}+zipWithTails1 ::+   (y0 -> T y1 -> y2) -> T y0 -> T y1 -> T y2+zipWithTails1 f xs ys =+   generate (\(xs0,ys0) ->+      do (x,xs1) <- viewL xs0+         ys1 <- ys0+         return (f x ys1, (xs1, fmap snd $ viewL ys1)))+      (xs, Just ys)++-- | in contrast to 'zipWithTails' it appends infinitely many empty suffixes+{-# INLINE zipWithTailsInf #-}+zipWithTailsInf ::+   (y0 -> T y1 -> y2) -> T y0 -> T y1 -> T y2+zipWithTailsInf f =+   curry $ generate (\(xs0,ys0) ->+      do (x,xs) <- viewL xs0+         return (f x ys0, (xs, switchL empty (flip const) ys0)))+ {- This can hardly be implemented in an efficient way. But this means, we cannot implement the Generic.Transform class.@@ -912,3 +949,17 @@          fmap (\b -> (b, (Nothing, s0))) carry `mplus`          fmap (\((a,b),s1) -> (a, (Just b, s1))) (next s0))       (Nothing,t))++interleave, interleaveAlt ::+   T y -> T y -> T y+interleave xs ys =+   runViewL xs (\nextX sx ->+   runViewL ys (\nextY sy ->+   unfoldR+      (\(select,(sx0,sy0)) ->+         case select of+            False -> fmap (mapSnd (\sx1 -> (True,  (sx1,sy0)))) $ nextX sx0+            True  -> fmap (mapSnd (\sy1 -> (False, (sx0,sy1)))) $ nextY sy0)+      (False, (sx,sy))))++interleaveAlt xs ys = flattenPairs $ zip xs ys
src/Synthesizer/Storable/Cut.hs view
@@ -12,16 +12,13 @@ import qualified Data.EventList.Relative.TimeMixed as EventListTM import qualified Data.EventList.Absolute.TimeBody  as AbsEventList import Control.Monad.Trans.State (runState, modify, gets, put, )--- import Control.Monad (mapM, ) import Data.Tuple.HT (mapSnd, ) --- import qualified Algebra.RealRing     as RealRing import qualified Algebra.Additive as Additive import qualified Number.NonNegative as NonNeg  import Foreign.Storable (Storable) -import NumericPrelude.Base import NumericPrelude.Numeric
src/Synthesizer/Storable/Filter/NonRecursive.hs view
@@ -30,8 +30,6 @@ import Foreign.Storable (Storable, ) import Foreign.Storable.Tuple () -import Algebra.Module( {- linearComb, -} (*>), )- import Control.Monad (mplus, ) import Data.Maybe.HT (toMaybe, ) import Data.Maybe (fromMaybe, )@@ -43,7 +41,6 @@  import NumericPrelude.Base import NumericPrelude.Numeric as NP-import qualified Prelude as P   {- |
src/Synthesizer/Storable/Oscillator.hs view
@@ -28,13 +28,9 @@ import Algebra.Module((*>)) -} import qualified Algebra.Transcendental        as Trans--- import qualified Algebra.Field                 as Field import qualified Algebra.RealRing              as RealRing-import qualified Algebra.Additive              as Additive  import NumericPrelude.Numeric--import qualified Prelude as P import NumericPrelude.Base
src/Test/Main.hs view
@@ -6,9 +6,16 @@ import qualified Test.Sound.Synthesizer.Plain.Interpolation  as Interpolation import qualified Test.Sound.Synthesizer.Plain.Oscillator     as Oscillator import qualified Test.Sound.Synthesizer.Plain.Wave           as Wave+import qualified Test.Sound.Synthesizer.Basic.NumberTheory   as NumberTheory import qualified Test.Sound.Synthesizer.Basic.ToneModulation as ToneModulation import qualified Test.Sound.Synthesizer.Plain.ToneModulation as ToneModulationL import qualified Test.Sound.Synthesizer.Generic.ToneModulation as ToneModulationG+import qualified Test.Sound.Synthesizer.Generic.Permutation as Permutation+import qualified Test.Sound.Synthesizer.Generic.Fourier as Fourier+import qualified Test.Sound.Synthesizer.Generic.FourierInteger as FourierInteger+import qualified Test.Sound.Synthesizer.Generic.Filter  as FilterG+import qualified Test.Sound.Synthesizer.Generic.Cut  as CutG+import qualified Test.Sound.Synthesizer.Causal.Analysis as AnalysisC import qualified Test.Sound.Synthesizer.Storable.Cut as Cut  import Data.Tuple.HT (mapFst, )@@ -29,7 +36,14 @@       prefix "Plain.Oscillator"     Oscillator.tests :       prefix "Plain.Wave"           Wave.tests :       prefix "Storable.Cut"         Cut.tests :+      prefix "Generic.Cut"          CutG.tests :       prefix "Basic.ToneModulation" ToneModulation.tests :       prefix "Plain.ToneModulation" ToneModulationL.tests :       prefix "Generic.ToneModulation" ToneModulationG.tests :+      prefix "Generic.Permutation"    Permutation.tests :+      prefix "Generic.Fourier"        Fourier.tests :+      prefix "Basic.NumberTheory"     NumberTheory.tests :+      prefix "Generic.FourierInteger" FourierInteger.tests :+      prefix "Generic.Filter"         FilterG.tests :+      prefix "Causal.Analysis"        AnalysisC.tests :       []
+ src/Test/Sound/Synthesizer/Basic/NumberTheory.hs view
@@ -0,0 +1,119 @@+{-# LANGUAGE NoImplicitPrelude #-}+module Test.Sound.Synthesizer.Basic.NumberTheory (tests) where++import Synthesizer.Basic.NumberTheory (Order(Order), )+import qualified Synthesizer.Basic.NumberTheory as NT+import qualified Data.Set as Set++import Test.QuickCheck (Testable, Arbitrary, arbitrary, quickCheck, )++import qualified Algebra.Absolute              as Absolute++import NumericPrelude.Numeric+import NumericPrelude.Base+import Prelude ()+++newtype Cardinal a = Cardinal a+   deriving (Show)++instance (Absolute.C a, Arbitrary a) => Arbitrary (Cardinal a) where+   arbitrary = fmap (Cardinal . abs) arbitrary+++newtype Positive a = Positive a+   deriving (Show)++instance (Absolute.C a, Arbitrary a) => Arbitrary (Positive a) where+   arbitrary = fmap (Positive . (1+) . abs) arbitrary+++simple ::+   (Testable t, Arbitrary (wrapper Integer), Show (wrapper Integer)) =>+   (wrapper Integer -> t) -> IO ()+simple = quickCheck++tests :: [(String, IO ())]+tests =+   ("primitiveRootsOfUnity naive vs. power",+      simple $ \(Cardinal m) order ->+         NT.primitiveRootsOfUnityNaive m order+         ==+         NT.primitiveRootsOfUnityPower m order) :+   ("primitiveRootsOfUnity naive vs. fullorbit",+      simple $ \(Cardinal m) order ->+         NT.primitiveRootsOfUnityNaive m order+         ==+         (Set.toAscList $ Set.fromList $+          NT.primitiveRootsOfUnityFullOrbit m order)) :+   ("Carmichael theorem",+      simple $ \(Positive a) (Positive b) ->+         NT.getOrder (NT.maximumOrderOfPrimitiveRootsOfUnity (lcm a b))+         ==+         lcm+            (NT.getOrder (NT.maximumOrderOfPrimitiveRootsOfUnity a))+            (NT.getOrder (NT.maximumOrderOfPrimitiveRootsOfUnity b))) :+   ("maximumOrderOfPrimitiveRootsOfUnity naive vs. integer",+      simple $ \(Positive m) ->+         NT.maximumOrderOfPrimitiveRootsOfUnityNaive m+         ==+         NT.maximumOrderOfPrimitiveRootsOfUnityInteger m) :+   ("number of rootsOfUnityPower, lcm",+      simple $ \(Positive m) ao@(Order a) bo@(Order b) ->+         let g = length . NT.rootsOfUnityPower m+         in  g (Order $ lcm a b) == lcm (g ao) (g bo)) :+   ("ringsWithPrimitiveRootsOfUnityAndUnits: minimal modulus",+      quickCheck $ \order@(Order expo) ->+         (head $ NT.ringsWithPrimitiveRootOfUnityAndUnit order)+         ==+         (head $ NT.ringsWithPrimitiveRootsOfUnityAndUnitsNaive+            [order] [expo])) :+   ("combine two rings with primitive roots of certain orders",+      quickCheck $ \m n ->+         let r = lcm+                   (head (NT.ringsWithPrimitiveRootOfUnityAndUnit m))+                   (head (NT.ringsWithPrimitiveRootOfUnityAndUnit n))+         in  NT.hasPrimitiveRootOfUnityInteger r m+             &&+             NT.hasPrimitiveRootOfUnityInteger r n) :+   ("combine many rings with primitive roots of certain orders",+      quickCheck $ \n0 ns0 ->+         let ns = take 3 $ map (\n -> 1 + mod n 10) (n0:ns0)+             order = NT.lcmMulti ns+         in  take 3 (NT.ringsWithPrimitiveRootsOfUnityAndUnitsNaive+                       (map Order ns) ns)+             ==+             take 3 (NT.ringsWithPrimitiveRootsOfUnityAndUnitsNaive+                       [Order order] [order])) :+{-+Unfortunately rings with certain units cannot be combined+while maintaining these elements as units.++Counterexample:+   ringsWithPrimitiveRootOfUnityAndUnit 2 = 3:...+   ringsWithPrimitiveRootOfUnityAndUnit 3 = 7:...+   But in Z_{3·7} the number 3 is no unit.++   ("combine rings with certain units",+      quickCheck $ \(Positive m) (Positive n) ->+         let r = fromIntegral $ lcm+                (head (NT.ringsWithPrimitiveRootOfUnityAndUnit m))+                (head (NT.ringsWithPrimitiveRootOfUnityAndUnit n))+         in  PID.coprime r m && PID.coprime r n) :+-}+   ("number of roots of unity lcm",+      quickCheck $ \(Positive n) (Positive k) (Positive l) ->+         let orders = NT.ordersOfRootsOfUnityInteger !! (n-1)+         in  lcm (orders!!(k-1)) (orders!!(l-1))+             ==+             orders !! (lcm k l - 1)) :+   ("number of roots of unity vs. primitive roots",+      quickCheck $ \(Positive n) (Positive k) ->+         (sum $ map snd $+          filter (flip divides k . fst) $+          zip+             [1..]+             (NT.ordersOfPrimitiveRootsOfUnityInteger !! (n-1)))+         ==+         NT.ordersOfRootsOfUnityInteger !! (n-1) !! (k-1)) :+   []
src/Test/Sound/Synthesizer/Basic/ToneModulation.hs view
@@ -12,14 +12,9 @@ -- import Test.Utility  import qualified Number.NonNegative       as NonNeg--- import qualified Number.NonNegativeChunky as Chunky --- import qualified Algebra.RealTranscendental    as RealTrans--- 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 NumericPrelude.Numeric
+ src/Test/Sound/Synthesizer/Causal/Analysis.hs view
@@ -0,0 +1,32 @@+module Test.Sound.Synthesizer.Causal.Analysis (tests) where++import qualified Synthesizer.Causal.Analysis as AnaC+import qualified Synthesizer.Causal.Process as Causal+import qualified Synthesizer.Plain.Analysis as Ana++import Control.Arrow ((<<<), )++import qualified Data.List.Match as Match++import Test.QuickCheck (quickCheck, )++import NumericPrelude.Numeric+import NumericPrelude.Base+import Prelude ()+++tests :: [(String, IO ())]+tests =+   ("deltaSigmaModulation",+      quickCheck $ \xs ->+         Match.take xs (Ana.deltaSigmaModulation xs)+         ==+         Causal.apply AnaC.deltaSigmaModulation (xs::[Rational])) :+   ("deltaSigmaModulationPositive",+      quickCheck $ \threshold xs ->+         Match.take xs (Ana.deltaSigmaModulationPositive threshold xs)+         ==+         Causal.apply+            (AnaC.deltaSigmaModulationPositive <<<+             Causal.feedConstFst threshold) (xs::[Rational])) :+   []
+ src/Test/Sound/Synthesizer/Generic/Cut.hs view
@@ -0,0 +1,104 @@+module Test.Sound.Synthesizer.Generic.Cut (tests) where++import qualified Synthesizer.Generic.Cut as CutG+import qualified Synthesizer.Generic.Signal as SigG++import qualified Synthesizer.Storable.Signal as SigSt++import qualified Synthesizer.ChunkySize.Signal as SigChunky+import qualified Synthesizer.ChunkySize as ChunkySize++import qualified Data.StorableVector as SV+import qualified Data.StorableVector.Lazy.Pattern as SVP++import qualified Synthesizer.State.Signal as SigS++import qualified Data.EventList.Relative.BodyTime as EventList++import qualified Number.NonNegative as NonNeg+import qualified Number.NonNegativeChunky as Chunky++import qualified Numeric.NonNegative.Wrapper as NonNeg98++import Data.Tuple.HT (mapSnd, )++import Test.QuickCheck (quickCheck, )++import NumericPrelude.Numeric+import NumericPrelude.Base+import Prelude ()+++dropMarginRemLength :: NonNeg.Int -> NonNeg.Int -> [Int] -> Bool+dropMarginRemLength nn nm xs =+   let n = NonNeg.toNumber nn+       m = NonNeg.toNumber nm+       (k,ys) = CutG.dropMarginRem n m xs+   in  length xs - m == length ys - k++dropMarginRemState :: NonNeg.Int -> NonNeg.Int -> [Int] -> Bool+dropMarginRemState nn nm xs =+   let n = NonNeg.toNumber nn+       m = NonNeg.toNumber nm+   in  CutG.dropMarginRem n m (SigS.fromList xs)+       ==+       mapSnd SigS.fromList (CutG.dropMarginRem n m xs)++dropMarginRemSV :: NonNeg.Int -> NonNeg.Int -> [Int] -> Bool+dropMarginRemSV nn nm xs =+   let n = NonNeg.toNumber nn+       m = NonNeg.toNumber nm+   in  CutG.dropMarginRem n m (SV.pack xs)+       ==+       mapSnd SV.pack (CutG.dropMarginRem n m xs)++dropMarginRemSVL :: NonNeg.Int -> NonNeg.Int -> ChunkySize.T -> [Int] -> Bool+dropMarginRemSVL nn nm pat xs =+   let n = NonNeg.toNumber nn+       m = NonNeg.toNumber nm+   in  CutG.dropMarginRem n m+          (CutG.take (CutG.length pat) xs)+       ==+       mapSnd SigG.toList+          (CutG.dropMarginRem n m+             (SigChunky.fromState pat $+              SigG.toState xs :: SigSt.T Int))++dropMarginRemChunkySize ::+   NonNeg.Int -> NonNeg.Int -> ChunkySize.T -> Int -> Bool+dropMarginRemChunkySize nn nm pat x =+   let n = NonNeg.toNumber nn+       m = NonNeg.toNumber nm+   in  CutG.dropMarginRem n m pat+       ==+       mapSnd+          (ChunkySize.fromStorableVectorSize . SVP.length)+          (CutG.dropMarginRem n m+             (SVP.replicate (ChunkySize.toStorableVectorSize pat) x))++dropMarginRemPiecewise ::+   NonNeg.Int -> NonNeg.Int -> ChunkySize.T -> Int -> Bool+dropMarginRemPiecewise nn nm pat x =+   let n = NonNeg.toNumber nn+       m = NonNeg.toNumber nm+   in  CutG.dropMarginRem n m pat+       ==+       mapSnd+          (Chunky.fromChunks .+           map (\size -> SigG.LazySize $ NonNeg98.toNumber size) .+           EventList.getTimes)+          (CutG.dropMarginRem n m+             (EventList.fromPairList $ map ((,) x) $+              map (\(SigG.LazySize size) -> NonNeg98.fromNumber size) $+              Chunky.toChunks pat))+++tests :: [(String, IO ())]+tests =+   ("dropMarginRemLength", quickCheck dropMarginRemLength) :+   ("dropMarginRemState", quickCheck dropMarginRemState) :+   ("dropMarginRemSV", quickCheck dropMarginRemSV) :+   ("dropMarginRemSVL", quickCheck dropMarginRemSVL) :+   ("dropMarginRemChunkySize", quickCheck dropMarginRemChunkySize) :+   ("dropMarginRemPiecewise", quickCheck dropMarginRemPiecewise) :+   []
+ src/Test/Sound/Synthesizer/Generic/Filter.hs view
@@ -0,0 +1,64 @@+{-# LANGUAGE NoImplicitPrelude #-}+module Test.Sound.Synthesizer.Generic.Filter (tests) where++import qualified Synthesizer.Generic.Filter.NonRecursive as FiltNRG+import qualified Synthesizer.Generic.Cyclic as Cyclic+import qualified Synthesizer.Generic.Signal as SigG+import qualified Synthesizer.Generic.Cut as CutG+import qualified Synthesizer.Plain.Signal as Sig++import qualified Test.Sound.Synthesizer.Plain.NonEmpty as NonEmpty++import Test.QuickCheck (Testable, quickCheck, )++-- import qualified Algebra.Ring                  as Ring++import qualified Algebra.Laws                  as Law++import NumericPrelude.Numeric+import NumericPrelude.Base+++simple ::+   (Testable t) =>+   (Sig.T Int -> t) -> IO ()+simple = quickCheck++(=|=) ::+   (Eq sig, CutG.Transform sig) =>+   sig -> sig -> Bool+x =|= y =+   CutG.take 100 x == CutG.take 100 y++tests :: [(String, IO ())]+tests =+   ("identity",+      simple $ Law.identity FiltNRG.generic $ SigG.singleton one) :+   ("commutativity",+      simple $ Law.commutative FiltNRG.generic) :+   ("distributivity",+      simple $ Law.leftDistributive FiltNRG.generic SigG.mix) :+   ("karatsuba finite",+      simple $ \x y -> FiltNRG.generic x y == FiltNRG.karatsubaFinite (*) x y) :+   ("karatsuba finite-infinite",+      simple $ \x y -> FiltNRG.generic x y == FiltNRG.karatsubaFiniteInfinite (*) x y) :+   ("karatsuba infinite",+      simple $ \x y -> FiltNRG.generic x y == FiltNRG.karatsubaInfinite (*) x y) :+   ("karatsuba finite-infinite cycle",+      simple $ \x yn ->+         case NonEmpty.toInfiniteList yn of+            y -> FiltNRG.generic x y =|= FiltNRG.karatsubaFiniteInfinite (*) x y) :+   ("karatsuba infinite cycle",+      simple $ \x yn ->+         case NonEmpty.toInfiniteList yn of+            y -> FiltNRG.generic x y =|= FiltNRG.karatsubaInfinite (*) x y) :+   ("convolve triple",+      quickCheck $ \x y ->+         Cyclic.sumAndConvolveTriple x y ==+         Cyclic.sumAndConvolveTripleAlt x (y :: Cyclic.Triple Integer)) :+   ("periodic summation",+      simple $ \x y n ->+         let periodic = Cyclic.fromSignal SigG.defaultLazySize (1 + abs n)+         in  Cyclic.convolve (periodic x) (periodic y) ==+             periodic (FiltNRG.generic x y)) :+   []
+ src/Test/Sound/Synthesizer/Generic/Fourier.hs view
@@ -0,0 +1,152 @@+{-# LANGUAGE NoImplicitPrelude #-}+module Test.Sound.Synthesizer.Generic.Fourier (tests) where++import qualified Synthesizer.Generic.Fourier as Fourier+import qualified Synthesizer.Generic.Cyclic as Cyclic+import qualified Synthesizer.Generic.Filter.NonRecursive as FiltNRG+import qualified Synthesizer.Generic.Analysis as AnaG+import qualified Synthesizer.Generic.Signal as SigG+import qualified Synthesizer.Generic.Signal2 as SigG2+import qualified Synthesizer.Generic.Cut as CutG+import qualified Synthesizer.Storable.Signal as SigSt+import qualified Synthesizer.State.Signal as SigS++import Test.QuickCheck (Testable, Arbitrary, arbitrary, quickCheck, )+import Test.Utility (approxEqualAbs, approxEqualComplexAbs, )++import qualified Number.Complex as Complex++import qualified Algebra.Ring                  as Ring+import qualified Algebra.Additive              as Additive++import Control.Monad (liftM2, )++import NumericPrelude.Numeric+import NumericPrelude.Base+import Prelude ()+++tolerance :: Double+tolerance = 1e-10++normalize ::+   SigSt.T (Complex.T Double) -> SigSt.T (Complex.T Double)+normalize xs =+   FiltNRG.amplifyVector+      (recip $ max (0.1::Double) $ AnaG.volumeVectorMaximum xs) xs++newtype Normed = Normed (SigSt.T (Complex.T Double))+   deriving (Show)++instance Arbitrary Normed where+   arbitrary = fmap (Normed . normalize) arbitrary+++data Normed2 =+      Normed2+         (SigSt.T (Complex.T Double))+         (SigSt.T (Complex.T Double))+   deriving (Show)++instance Arbitrary Normed2 where+   arbitrary =+      liftM2+         (\x y ->+            let len = min (CutG.length x) (CutG.length y)+            in  Normed2+                   (normalize $ CutG.take len x)+                   (normalize $ CutG.take len y))+         arbitrary+         arbitrary+++-- could be moved to NumericPrelude+class Complex a where+   conjugate :: a -> a++instance (Additive.C a) => Complex (Complex.T a) where+   conjugate = Complex.conjugate++scalarProduct ::+   (SigG.Read sig y, Ring.C y, Complex y) =>+   sig y -> sig y -> y+scalarProduct xs ys =+   SigS.sum $+   SigS.zipWith (*)+      (SigG.toState xs)+      (SigS.map conjugate $ SigG.toState ys)++(=~=) ::+   SigSt.T (Complex.T Double) ->+   SigSt.T (Complex.T Double) ->+   Bool+(=~=) xs ys =+   SigG.length xs == SigG.length ys &&+   (SigG.foldR (&&) True $+    SigG2.zipWith (approxEqualComplexAbs tolerance) xs ys)++simple ::+   (Testable t) =>+   (SigSt.T (Complex.T Double) -> t) -> IO ()+simple = quickCheck++tests :: [(String, IO ())]+tests =+   ("fourier inverse",+      quickCheck $ \(Normed x) ->+         x =~=+         (FiltNRG.amplify (recip $ fromIntegral $ SigG.length x) $+          Fourier.transformBackward $ Fourier.transformForward x)) :+   ("double fourier = reverse",+      quickCheck $ \(Normed x) ->+         x =~=+         (Cyclic.reverse $+          FiltNRG.amplify (recip $ fromIntegral $ SigG.length x) $+          Fourier.transformForward $+          Fourier.transformForward x)) :+   ("fourier of reverse",+      quickCheck $ \(Normed x) ->+         Cyclic.reverse (Fourier.transformForward x) =~=+         Fourier.transformForward (Cyclic.reverse x)) :+   ("fourier of conjugate",+      quickCheck $ \(Normed x) ->+         (SigG.map Complex.conjugate $ Fourier.transformForward x)+         =~=+         (Fourier.transformForward $+          SigG.map Complex.conjugate $ Cyclic.reverse x)) :+   ("additivity",+      quickCheck $ \(Normed2 x y) ->+         SigG.mix (Fourier.transformForward x) (Fourier.transformForward y)+         =~=+         Fourier.transformForward (SigG.mix x y)) :+   ("isometry",+      simple $ \xs x0 ->+         let x = normalize (SigG.cons x0 xs)+         in  approxEqualAbs tolerance+                (AnaG.volumeVectorEuclideanSqr $ Fourier.transformForward x)+                (fromIntegral (SigG.length x) *+                 AnaG.volumeVectorEuclideanSqr x)) :+   ("unitarity",+      quickCheck $ \(Normed2 x y) ->+         approxEqualComplexAbs tolerance+            (scalarProduct+               (Fourier.transformForward x) (Fourier.transformForward y))+            (fromIntegral (SigG.length x) * scalarProduct x y)) :+   ("convolution",+      quickCheck $ \(Normed2 x y) ->+         SigG.zipWith (*)+            (Fourier.transformForward x)+            (Fourier.transformForward y)+         =~=+         Fourier.transformForward (Cyclic.convolve x y)) :+   ("convolution cyclic",+      quickCheck $ \(Normed2 x y) ->+         Fourier.convolveCyclic x y+         =~=+         Cyclic.convolve x y) :+   ("convolution long",+      quickCheck $ \(Normed x) (Normed y) ->+         FiltNRG.karatsubaFinite (*) x y+         =~=+         Fourier.convolveWithWindow (Fourier.window x) y) :+   []
+ src/Test/Sound/Synthesizer/Generic/FourierInteger.hs view
@@ -0,0 +1,178 @@+{-# LANGUAGE NoImplicitPrelude #-}+module Test.Sound.Synthesizer.Generic.FourierInteger (tests) where++import qualified Synthesizer.Generic.Fourier as Fourier+import qualified Synthesizer.Generic.Cyclic as Cyclic+import qualified Synthesizer.Generic.Filter.NonRecursive as FiltNRG+import qualified Synthesizer.Generic.Signal as SigG+import qualified Synthesizer.Generic.Cut as CutG+import qualified Synthesizer.State.Signal as SigS+import qualified Synthesizer.Plain.Signal as Sig++import Test.QuickCheck (Testable, Arbitrary, arbitrary, quickCheck, )++import qualified Synthesizer.Basic.NumberTheory as NT++import qualified Number.ResidueClass.Check as RC+import Number.ResidueClass.Check ((/:), )++import qualified Algebra.ToInteger             as ToInteger+import qualified Algebra.IntegralDomain        as Integral+import qualified Algebra.Ring                  as Ring++import Control.Monad (liftM2, )++import NumericPrelude.Numeric+import NumericPrelude.Base+import Prelude ()+++newtype Modulus a = Modulus a+   deriving (Show)++instance Ring.C a => Arbitrary (Modulus a) where+   arbitrary = fmap (Modulus . (2+) . fromInteger) arbitrary+++data ModularSignal =+      ModularSignal (Modulus Integer) (Sig.T (RC.T Integer))+   deriving (Show)++instance Arbitrary ModularSignal where+   arbitrary =+      fmap (uncurry ModularSignal . signal) arbitrary+++data ModularSignal2 =+      ModularSignal2+         (Modulus Integer) (Sig.T (RC.T Integer)) (Sig.T (RC.T Integer))+   deriving (Show)++instance Arbitrary ModularSignal2 where+   arbitrary =+      liftM2+         (\x y ->+            let len = min (CutG.length x) (CutG.length y)+                m = NT.fastFourierRing len+            in  ModularSignal2+                   (Modulus m)+                   (fmap (/: m) $ CutG.take len x)+                   (fmap (/: m) $ CutG.take len y))+         arbitrary+         arbitrary++scalarProduct ::+   Modulus Integer ->+   Sig.T (RC.T Integer) -> Sig.T (RC.T Integer) ->+   RC.T Integer+scalarProduct (Modulus m) xs ys =+   SigS.foldL (+) (RC.zero m) $+   SigS.zipWith (*)+      (SigG.toState xs)+      (SigG.toState ys)++{-+signal ::+   Integral.C a =>+   Modulus a -> Sig.T a -> Sig.T (RC.T a)+signal (Modulus a) = fmap (/: a)+-}++signal ::+   Sig.T Integer -> (Modulus Integer, Sig.T (RC.T Integer))+signal xs =+   let m = NT.fastFourierRing $ length xs+   in  (Modulus m, fmap (/: m) xs)++modular ::+   (Integral.C a, ToInteger.C b) =>+   Modulus a -> b -> RC.T a+modular (Modulus m) =+   RC.fromRepresentative m . fromIntegral+++simple ::+   (Testable t) =>+   (Sig.T Integer -> t) -> IO ()+simple = quickCheck++tests :: [(String, IO ())]+tests =+   ("fourier inverse",+      quickCheck $ \(ModularSignal m x) ->+         (Fourier.transformBackward $ Fourier.transformForward x)+         ==+         FiltNRG.amplify (modular m $ length x) x) :+   ("double fourier = reverse",+      quickCheck $ \(ModularSignal m x) ->+         (Cyclic.reverse $+          Fourier.transformForward $+          Fourier.transformForward x)+         ==+         FiltNRG.amplify (modular m $ length x) x) :+   ("fourier of reverse",+      quickCheck $ \(ModularSignal _m x) ->+         Cyclic.reverse (Fourier.transformForward x) ==+         Fourier.transformForward (Cyclic.reverse x)) :+   ("homogenity",+      quickCheck $ \(ModularSignal m x) y ->+         (FiltNRG.amplify (modular m (y::Integer)) $+          Fourier.transformForward x)+         ==+         (Fourier.transformForward $+          FiltNRG.amplify (modular m y) x)) :+   ("additivity",+      quickCheck $ \(ModularSignal2 _m x y) ->+         SigG.mix (Fourier.transformForward x) (Fourier.transformForward y)+         ==+         Fourier.transformForward (SigG.mix x y)) :+{-+   ("isometry",+      simple $ \xs x0 ->+         let (m,x) = signal (SigG.cons x0 xs)+         in  (AnaG.volumeVectorEuclideanSqr $ Fourier.transformForward x)+             ==+             (modular m (SigG.length x) *+              AnaG.volumeVectorEuclideanSqr x)) :+-}+   ("unitarity",+      quickCheck $ \(ModularSignal2 m x y) ->+         {-+         since there is no equivalent of a complex conjugate+         we have to take the scalar product with the backwards transform.+         -}+         scalarProduct m+            (Fourier.transformForward x) (Fourier.transformBackward y)+         ==+         modular m (length x) * scalarProduct m x y) :+   ("convolution",+      quickCheck $ \(ModularSignal2 _m x y) ->+         SigG.zipWith (*)+            (Fourier.transformForward x)+            (Fourier.transformForward y)+         ==+         Fourier.transformForward (Cyclic.convolve x y)) :+   ("convolution cyclic",+      quickCheck $ \(ModularSignal2 _m x y) ->+         Fourier.convolveCyclic x y+         ==+         Cyclic.convolve x y) :+   ("convolution long",+      simple $ \x0 y0 ->+         let m = Modulus $ NT.fastFourierRing $+                 2 * (NT.ceilingPowerOfTwo $ length x0)+             x = fmap (modular m) x0+             y = fmap (modular m) y0+         in  fmap (modular m) (FiltNRG.karatsubaFinite (*) x0 y0)+             ==+             Fourier.convolveWithWindow (Fourier.window x) y) :+   ("convolution long modular",+      simple $ \x0 y0 ->+         let m = Modulus $ NT.fastFourierRing $+                 2 * (NT.ceilingPowerOfTwo $ length x0)+             x = fmap (modular m) x0+             y = fmap (modular m) (y0 :: Sig.T Integer)+         in  FiltNRG.karatsubaFinite (*) x y+             ==+             Fourier.convolveWithWindow (Fourier.window x) y) :+   []
+ src/Test/Sound/Synthesizer/Generic/Permutation.hs view
@@ -0,0 +1,45 @@+{-+wish list:+ - custom Permutation type with Arbitrary instance+-}+{-# LANGUAGE NoImplicitPrelude #-}+module Test.Sound.Synthesizer.Generic.Permutation (tests) where++import qualified Synthesizer.Generic.Permutation as Permutation++import Test.QuickCheck (quickCheck, )++import NumericPrelude.Numeric+import NumericPrelude.Base+import Prelude ()+++tests :: [(String, IO ())]+tests =+   ("inverse transposition",+      quickCheck $ \n0 m0 ->+         let n = mod n0 100+             m = mod m0 100+         in  Permutation.inverse (Permutation.transposition n m)+             ==+             Permutation.transposition m n) :+   ("inverse skewGrid",+      quickCheck $ \n0 m0 ->+         let g = gcd n0 m0+             (n,m) = if g==0 then (0,0) else (abs (div n0 g), abs (div m0 g))+         in  Permutation.inverse (Permutation.skewGrid n m)+             ==+             Permutation.skewGridInv n m) :+   ("inverse skewGridCRT",+      quickCheck $ \n0 m0 ->+         let g = gcd n0 m0+             (n,m) = if g==0 then (0,0) else (abs (div n0 g), abs (div m0 g))+         in  Permutation.inverse (Permutation.skewGridCRT n m)+             ==+             Permutation.skewGridCRTInv n m) :+   {-+   reverse (multiplicative (generator n) n)+   ==+   multiplicative (recip $ generator n) n+   -}+   []
src/Test/Sound/Synthesizer/Generic/ToneModulation.hs view
@@ -35,12 +35,7 @@  import qualified Number.NonNegative       as NonNeg --- import qualified Algebra.RealTranscendental    as RealTrans--- 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 Data.List.HT (viewL, takeWhileJust, )
src/Test/Sound/Synthesizer/Plain/Analysis.hs view
@@ -6,8 +6,6 @@ import qualified Algebra.RealField             as RealField import qualified Algebra.Field                 as Field import qualified Algebra.RealRing              as RealRing-import qualified Algebra.Ring                  as Ring-import qualified Algebra.Additive              as Additive  import qualified Algebra.NormedSpace.Maximum   as NormedMax import qualified Algebra.NormedSpace.Euclidean as NormedEuc
src/Test/Sound/Synthesizer/Plain/Oscillator.hs view
@@ -8,16 +8,8 @@ -- import qualified Test.Sound.Synthesizer.Plain.Interpolation as InterpolationTest  import Test.QuickCheck (test, {- Property, (==>), -} )--- import Test.Utility --- import qualified Number.NonNegative       as NonNeg---- import qualified Algebra.RealTranscendental    as RealTrans--- 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 NumericPrelude.Numeric
src/Test/Sound/Synthesizer/Plain/ToneModulation.hs view
@@ -20,7 +20,7 @@ import qualified Test.Sound.Synthesizer.Plain.NonEmpty as NonEmpty import qualified Test.Sound.Synthesizer.Plain.Interpolation as InterpolationTest -import Test.QuickCheck (quickCheck, Property, (==>), Arbitrary, arbitrary, )+import Test.QuickCheck (quickCheck, Property, (==>), ) import Test.Utility (ArbChar, )  import qualified Number.NonNegative       as NonNeg@@ -29,11 +29,8 @@ import qualified Algebra.RealTranscendental    as RealTrans 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 Control.Monad (liftM2, ) import Data.List.HT (isAscending, ) import Data.Ord.HT (limit, ) import Data.Tuple.HT (mapPair, mapSnd, )
src/Test/Sound/Synthesizer/Plain/Wave.hs view
@@ -9,11 +9,7 @@ import qualified Number.NonNegative       as NonNeg  import qualified Algebra.RealTranscendental    as RealTrans--- 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 Control.Monad (liftM, liftM2, ) import System.Random (Random)
src/Test/Sound/Synthesizer/Storable/Cut.hs view
@@ -17,9 +17,6 @@ import Test.QuickCheck (quickCheck, ) import Test.Utility (equalList, ) --- import qualified Algebra.Ring                  as Ring--- import qualified Algebra.Additive              as Additive- import NumericPrelude.Numeric import NumericPrelude.Base import Prelude ()
src/Test/Utility.hs view
@@ -3,9 +3,11 @@  import Test.QuickCheck (Arbitrary(arbitrary)) -import qualified Algebra.RealRing                  as RealRing-import qualified Algebra.Ring                  as Ring+import qualified Number.Complex as Complex +import qualified Algebra.RealRing              as RealRing++import qualified Data.List.HT as ListHT import qualified Data.Char as Char  import NumericPrelude.Base@@ -14,14 +16,17 @@  equalList :: Eq a => [a] -> Bool equalList xs =-   -- 'drop 1' instead of 'take' for suppression of error-   and (zipWith (==) xs (drop 1 xs))+   and (ListHT.mapAdjacent (==) xs)   approxEqual :: (RealRing.C a) => a -> a -> a -> Bool approxEqual eps x y =    2 * abs (x-y) <= eps * (abs x + abs y) +approxEqualAbs :: (RealRing.C a) => a -> a -> a -> Bool+approxEqualAbs eps x y =+   abs (x-y) <= eps+ approxEqualListRel :: (RealRing.C a) => a -> [a] -> Bool approxEqualListRel eps xs =    let n = fromIntegral $ length xs@@ -32,6 +37,20 @@    let n = fromIntegral $ length xs        s = sum xs    in  sum (map (\x -> abs (n*x-s)) xs)  <=  eps+++approxEqualComplex ::+   (RealRing.C a) =>+   a -> Complex.T a -> Complex.T a -> Bool+approxEqualComplex eps x y =+   2 * Complex.magnitudeSqr (x-y)+      <= eps^2 * (Complex.magnitudeSqr x + Complex.magnitudeSqr y)++approxEqualComplexAbs ::+   (RealRing.C a) =>+   a -> Complex.T a -> Complex.T a -> Bool+approxEqualComplexAbs eps x y =+   Complex.magnitudeSqr (x-y) <= eps^2   -- see event-list
synthesizer-core.cabal view
@@ -1,5 +1,5 @@ Name:           synthesizer-core-Version:        0.4.0.4+Version:        0.4.1 License:        GPL License-File:   LICENSE Author:         Henning Thielemann <haskell@henning-thielemann.de>@@ -18,6 +18,10 @@ -- the Overview module does not really fit into one of the part packages --   For an overview of the organization of the package --   and the discussion of various design issues see "Synthesizer.Overview".+   .+   Functions:+     Oscillators, Noise generators, Frequency filters,+     Fast Fourier transform for computation of frequency spectrum Stability:      Experimental Tested-With:    GHC==6.4.1, GHC==6.8.2, GHC==6.10.4 Cabal-Version:  >=1.6@@ -48,7 +52,7 @@   Source-Repository this-  Tag:         0.4.0.4+  Tag:         0.4.1   Type:        darcs   Location:    http://code.haskell.org/synthesizer/core/ @@ -63,7 +67,8 @@     transformers >=0.2 && <0.3,     event-list >=0.1 && <0.2,     non-negative >=0.1 && <0.2,-    numeric-prelude >=0.2 && <0.3,+    explicit-exception >=0.1.6 && <0.2,+    numeric-prelude >=0.2.1 && <0.3,     numeric-quest >=0.1 && <0.2,     utility-ht >=0.0.5 && <0.1,     filepath >=1.1 && <1.2,@@ -85,7 +90,7 @@       Build-Depends: base >= 3 && <4     Build-Depends:       array >=0.1 && <0.4,-      containers >=0.1 && <0.4,+      containers >=0.1 && <0.5,       random >=1.0 && <2.0,       process >=1.0 && <1.1   Else@@ -184,23 +189,27 @@     Synthesizer.State.ToneModulation     Synthesizer.Causal.Process     Synthesizer.Causal.Arrow+    Synthesizer.Causal.Analysis     Synthesizer.Causal.Cut     Synthesizer.Causal.Displacement     Synthesizer.Causal.Interpolation     Synthesizer.Causal.Oscillator     Synthesizer.Causal.Oscillator.Core     Synthesizer.Causal.ToneModulation+    Synthesizer.Causal.Spatial     Synthesizer.Causal.Filter.NonRecursive     Synthesizer.Causal.Filter.Recursive.Integration     Synthesizer.Generic.Analysis     Synthesizer.Generic.Cut     Synthesizer.Generic.Control+    Synthesizer.Generic.Cyclic     Synthesizer.Generic.Displacement     Synthesizer.Generic.Filter.NonRecursive     Synthesizer.Generic.Filter.Delay     Synthesizer.Generic.Filter.Recursive.Integration     Synthesizer.Generic.Filter.Recursive.MovingAverage     Synthesizer.Generic.Filter.Recursive.Comb+    Synthesizer.Generic.Fourier     Synthesizer.Generic.Interpolation     Synthesizer.Generic.Loop     Synthesizer.Generic.Noise@@ -219,6 +228,9 @@    Other-Modules:     Synthesizer.Basic.ComplexModule+    Synthesizer.Basic.NumberTheory+    Synthesizer.Generic.Permutation+    Synthesizer.Generic.LengthSignal   Executable test@@ -243,9 +255,37 @@     Test.Sound.Synthesizer.Plain.ToneModulation     Test.Sound.Synthesizer.Plain.Wave     Test.Sound.Synthesizer.Basic.ToneModulation+    Test.Sound.Synthesizer.Basic.NumberTheory+    Test.Sound.Synthesizer.Generic.Cut     Test.Sound.Synthesizer.Generic.ToneModulation+    Test.Sound.Synthesizer.Generic.Permutation+    Test.Sound.Synthesizer.Generic.Fourier+    Test.Sound.Synthesizer.Generic.FourierInteger+    Test.Sound.Synthesizer.Generic.Filter     Test.Sound.Synthesizer.Storable.Cut+    Test.Sound.Synthesizer.Causal.Analysis   Main-Is: Test/Main.hs+++Executable fouriertest+  If flag(buildProfilers)+    Build-Depends:+      storablevector >=0.2.7 && <0.3,+      utility-ht >=0.0.5 && <0.1,+      storable-tuple >=0.0.1 && <0.1,+      numeric-prelude >=0.2 && <0.3,+      timeit >=1.0 && <1.1,+      base >=4 && <5+  Else+    Buildable: False++  GHC-Options:    -Wall -auto-all+  Hs-Source-Dirs: speedtest, src+  If flag(category)+    Hs-Source-Dirs: src-4+  Else+    Hs-Source-Dirs: src-3+  Main-Is:        Fourier.hs  Executable speedtest   If !flag(buildProfilers)