diff --git a/speedtest/Fourier.hs b/speedtest/Fourier.hs new file mode 100644 --- /dev/null +++ b/speedtest/Fourier.hs @@ -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) : + []) +-} diff --git a/src/Synthesizer/Basic/ComplexModule.hs b/src/Synthesizer/Basic/ComplexModule.hs --- a/src/Synthesizer/Basic/ComplexModule.hs +++ b/src/Synthesizer/Basic/ComplexModule.hs @@ -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 #-} diff --git a/src/Synthesizer/Basic/Distortion.hs b/src/Synthesizer/Basic/Distortion.hs --- a/src/Synthesizer/Basic/Distortion.hs +++ b/src/Synthesizer/Basic/Distortion.hs @@ -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, ) diff --git a/src/Synthesizer/Basic/DistortionControlled.hs b/src/Synthesizer/Basic/DistortionControlled.hs --- a/src/Synthesizer/Basic/DistortionControlled.hs +++ b/src/Synthesizer/Basic/DistortionControlled.hs @@ -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, ) diff --git a/src/Synthesizer/Basic/NumberTheory.hs b/src/Synthesizer/Basic/NumberTheory.hs new file mode 100644 --- /dev/null +++ b/src/Synthesizer/Basic/NumberTheory.hs @@ -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=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 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 ( Integer +ceiling5Smooth n = + head $ dropWhile ( 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) diff --git a/src/Synthesizer/Basic/ToneModulation.hs b/src/Synthesizer/Basic/ToneModulation.hs --- a/src/Synthesizer/Basic/ToneModulation.hs +++ b/src/Synthesizer/Basic/ToneModulation.hs @@ -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 diff --git a/src/Synthesizer/Causal/Analysis.hs b/src/Synthesizer/Causal/Analysis.hs new file mode 100644 --- /dev/null +++ b/src/Synthesizer/Causal/Analysis.hs @@ -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) diff --git a/src/Synthesizer/Causal/Displacement.hs b/src/Synthesizer/Causal/Displacement.hs --- a/src/Synthesizer/Causal/Displacement.hs +++ b/src/Synthesizer/Causal/Displacement.hs @@ -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 diff --git a/src/Synthesizer/Causal/Interpolation.hs b/src/Synthesizer/Causal/Interpolation.hs --- a/src/Synthesizer/Causal/Interpolation.hs +++ b/src/Synthesizer/Causal/Interpolation.hs @@ -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 diff --git a/src/Synthesizer/Causal/Oscillator.hs b/src/Synthesizer/Causal/Oscillator.hs --- a/src/Synthesizer/Causal/Oscillator.hs +++ b/src/Synthesizer/Causal/Oscillator.hs @@ -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 diff --git a/src/Synthesizer/Causal/Oscillator/Core.hs b/src/Synthesizer/Causal/Oscillator/Core.hs --- a/src/Synthesizer/Causal/Oscillator/Core.hs +++ b/src/Synthesizer/Causal/Oscillator/Core.hs @@ -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 diff --git a/src/Synthesizer/Causal/Spatial.hs b/src/Synthesizer/Causal/Spatial.hs new file mode 100644 --- /dev/null +++ b/src/Synthesizer/Causal/Spatial.hs @@ -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) diff --git a/src/Synthesizer/ChunkySize/Cut.hs b/src/Synthesizer/ChunkySize/Cut.hs --- a/src/Synthesizer/ChunkySize/Cut.hs +++ b/src/Synthesizer/ChunkySize/Cut.hs @@ -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 diff --git a/src/Synthesizer/ChunkySize/Signal.hs b/src/Synthesizer/ChunkySize/Signal.hs --- a/src/Synthesizer/ChunkySize/Signal.hs +++ b/src/Synthesizer/ChunkySize/Signal.hs @@ -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 diff --git a/src/Synthesizer/Generic/Analysis.hs b/src/Synthesizer/Generic/Analysis.hs --- a/src/Synthesizer/Generic/Analysis.hs +++ b/src/Synthesizer/Generic/Analysis.hs @@ -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(xupper) -{- {- | 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 diff --git a/src/Synthesizer/Generic/Control.hs b/src/Synthesizer/Generic/Control.hs --- a/src/Synthesizer/Generic/Control.hs +++ b/src/Synthesizer/Generic/Control.hs @@ -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 diff --git a/src/Synthesizer/Generic/Cut.hs b/src/Synthesizer/Generic/Cut.hs --- a/src/Synthesizer/Generic/Cut.hs +++ b/src/Synthesizer/Generic/Cut.hs @@ -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 #-} diff --git a/src/Synthesizer/Generic/Cyclic.hs b/src/Synthesizer/Generic/Cyclic.hs new file mode 100644 --- /dev/null +++ b/src/Synthesizer/Generic/Cyclic.hs @@ -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)) diff --git a/src/Synthesizer/Generic/Filter/NonRecursive.hs b/src/Synthesizer/Generic/Filter/NonRecursive.hs --- a/src/Synthesizer/Generic/Filter/NonRecursive.hs +++ b/src/Synthesizer/Generic/Filter/NonRecursive.hs @@ -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)) diff --git a/src/Synthesizer/Generic/Filter/Recursive/Comb.hs b/src/Synthesizer/Generic/Filter/Recursive/Comb.hs --- a/src/Synthesizer/Generic/Filter/Recursive/Comb.hs +++ b/src/Synthesizer/Generic/Filter/Recursive/Comb.hs @@ -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 diff --git a/src/Synthesizer/Generic/Fourier.hs b/src/Synthesizer/Generic/Fourier.hs new file mode 100644 --- /dev/null +++ b/src/Synthesizer/Generic/Fourier.hs @@ -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) +-} diff --git a/src/Synthesizer/Generic/Interpolation.hs b/src/Synthesizer/Generic/Interpolation.hs --- a/src/Synthesizer/Generic/Interpolation.hs +++ b/src/Synthesizer/Generic/Interpolation.hs @@ -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 diff --git a/src/Synthesizer/Generic/LengthSignal.hs b/src/Synthesizer/Generic/LengthSignal.hs new file mode 100644 --- /dev/null +++ b/src/Synthesizer/Generic/LengthSignal.hs @@ -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 diff --git a/src/Synthesizer/Generic/Loop.hs b/src/Synthesizer/Generic/Loop.hs --- a/src/Synthesizer/Generic/Loop.hs +++ b/src/Synthesizer/Generic/Loop.hs @@ -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 diff --git a/src/Synthesizer/Generic/Noise.hs b/src/Synthesizer/Generic/Noise.hs --- a/src/Synthesizer/Generic/Noise.hs +++ b/src/Synthesizer/Generic/Noise.hs @@ -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 diff --git a/src/Synthesizer/Generic/Oscillator.hs b/src/Synthesizer/Generic/Oscillator.hs --- a/src/Synthesizer/Generic/Oscillator.hs +++ b/src/Synthesizer/Generic/Oscillator.hs @@ -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 diff --git a/src/Synthesizer/Generic/Permutation.hs b/src/Synthesizer/Generic/Permutation.hs new file mode 100644 --- /dev/null +++ b/src/Synthesizer/Generic/Permutation.hs @@ -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 diff --git a/src/Synthesizer/Generic/Piece.hs b/src/Synthesizer/Generic/Piece.hs --- a/src/Synthesizer/Generic/Piece.hs +++ b/src/Synthesizer/Generic/Piece.hs @@ -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 diff --git a/src/Synthesizer/Generic/Signal.hs b/src/Synthesizer/Generic/Signal.hs --- a/src/Synthesizer/Generic/Signal.hs +++ b/src/Synthesizer/Generic/Signal.hs @@ -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 diff --git a/src/Synthesizer/Generic/Signal2.hs b/src/Synthesizer/Generic/Signal2.hs --- a/src/Synthesizer/Generic/Signal2.hs +++ b/src/Synthesizer/Generic/Signal2.hs @@ -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 diff --git a/src/Synthesizer/Interpolation/Class.hs b/src/Synthesizer/Interpolation/Class.hs --- a/src/Synthesizer/Interpolation/Class.hs +++ b/src/Synthesizer/Interpolation/Class.hs @@ -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 diff --git a/src/Synthesizer/Interpolation/Custom.hs b/src/Synthesizer/Interpolation/Custom.hs --- a/src/Synthesizer/Interpolation/Custom.hs +++ b/src/Synthesizer/Interpolation/Custom.hs @@ -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 ((+.*), ) diff --git a/src/Synthesizer/Interpolation/Module.hs b/src/Synthesizer/Interpolation/Module.hs --- a/src/Synthesizer/Interpolation/Module.hs +++ b/src/Synthesizer/Interpolation/Module.hs @@ -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, ) diff --git a/src/Synthesizer/Plain/Analysis.hs b/src/Synthesizer/Plain/Analysis.hs --- a/src/Synthesizer/Plain/Analysis.hs +++ b/src/Synthesizer/Plain/Analysis.hs @@ -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 diff --git a/src/Synthesizer/Plain/Builder.hs b/src/Synthesizer/Plain/Builder.hs --- a/src/Synthesizer/Plain/Builder.hs +++ b/src/Synthesizer/Plain/Builder.hs @@ -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 diff --git a/src/Synthesizer/Plain/Control.hs b/src/Synthesizer/Plain/Control.hs --- a/src/Synthesizer/Plain/Control.hs +++ b/src/Synthesizer/Plain/Control.hs @@ -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 diff --git a/src/Synthesizer/Plain/Effect.hs b/src/Synthesizer/Plain/Effect.hs --- a/src/Synthesizer/Plain/Effect.hs +++ b/src/Synthesizer/Plain/Effect.hs @@ -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, diff --git a/src/Synthesizer/Plain/Effect/Fly.hs b/src/Synthesizer/Plain/Effect/Fly.hs --- a/src/Synthesizer/Plain/Effect/Fly.hs +++ b/src/Synthesizer/Plain/Effect/Fly.hs @@ -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 diff --git a/src/Synthesizer/Plain/Effect/Glass.hs b/src/Synthesizer/Plain/Effect/Glass.hs --- a/src/Synthesizer/Plain/Effect/Glass.hs +++ b/src/Synthesizer/Plain/Effect/Glass.hs @@ -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 diff --git a/src/Synthesizer/Plain/File.hs b/src/Synthesizer/Plain/File.hs --- a/src/Synthesizer/Plain/File.hs +++ b/src/Synthesizer/Plain/File.hs @@ -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 diff --git a/src/Synthesizer/Plain/Filter/Delay.hs b/src/Synthesizer/Plain/Filter/Delay.hs --- a/src/Synthesizer/Plain/Filter/Delay.hs +++ b/src/Synthesizer/Plain/Filter/Delay.hs @@ -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 diff --git a/src/Synthesizer/Plain/Filter/Delay/Block.hs b/src/Synthesizer/Plain/Filter/Delay/Block.hs --- a/src/Synthesizer/Plain/Filter/Delay/Block.hs +++ b/src/Synthesizer/Plain/Filter/Delay/Block.hs @@ -18,7 +18,6 @@ import Test.QuickCheck ((==>), Property) -import qualified Prelude as P import NumericPrelude.Base import NumericPrelude.Numeric diff --git a/src/Synthesizer/Plain/Filter/Delay/List.hs b/src/Synthesizer/Plain/Filter/Delay/List.hs --- a/src/Synthesizer/Plain/Filter/Delay/List.hs +++ b/src/Synthesizer/Plain/Filter/Delay/List.hs @@ -8,7 +8,6 @@ import Data.List(tails) -import qualified Prelude as P import NumericPrelude.Base import NumericPrelude.Numeric diff --git a/src/Synthesizer/Plain/Filter/Delay/ST.hs b/src/Synthesizer/Plain/Filter/Delay/ST.hs --- a/src/Synthesizer/Plain/Filter/Delay/ST.hs +++ b/src/Synthesizer/Plain/Filter/Delay/ST.hs @@ -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 diff --git a/src/Synthesizer/Plain/Filter/LinearPredictive.hs b/src/Synthesizer/Plain/Filter/LinearPredictive.hs --- a/src/Synthesizer/Plain/Filter/LinearPredictive.hs +++ b/src/Synthesizer/Plain/Filter/LinearPredictive.hs @@ -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 diff --git a/src/Synthesizer/Plain/Filter/NonRecursive.hs b/src/Synthesizer/Plain/Filter/NonRecursive.hs --- a/src/Synthesizer/Plain/Filter/NonRecursive.hs +++ b/src/Synthesizer/Plain/Filter/NonRecursive.hs @@ -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, ) diff --git a/src/Synthesizer/Plain/Filter/Recursive.hs b/src/Synthesizer/Plain/Filter/Recursive.hs --- a/src/Synthesizer/Plain/Filter/Recursive.hs +++ b/src/Synthesizer/Plain/Filter/Recursive.hs @@ -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 diff --git a/src/Synthesizer/Plain/Filter/Recursive/Allpass.hs b/src/Synthesizer/Plain/Filter/Recursive/Allpass.hs --- a/src/Synthesizer/Plain/Filter/Recursive/Allpass.hs +++ b/src/Synthesizer/Plain/Filter/Recursive/Allpass.hs @@ -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, ) diff --git a/src/Synthesizer/Plain/Filter/Recursive/AllpassPoly.hs b/src/Synthesizer/Plain/Filter/Recursive/AllpassPoly.hs --- a/src/Synthesizer/Plain/Filter/Recursive/AllpassPoly.hs +++ b/src/Synthesizer/Plain/Filter/Recursive/AllpassPoly.hs @@ -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 diff --git a/src/Synthesizer/Plain/Filter/Recursive/Butterworth.hs b/src/Synthesizer/Plain/Filter/Recursive/Butterworth.hs --- a/src/Synthesizer/Plain/Filter/Recursive/Butterworth.hs +++ b/src/Synthesizer/Plain/Filter/Recursive/Butterworth.hs @@ -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 diff --git a/src/Synthesizer/Plain/Filter/Recursive/Chebyshev.hs b/src/Synthesizer/Plain/Filter/Recursive/Chebyshev.hs --- a/src/Synthesizer/Plain/Filter/Recursive/Chebyshev.hs +++ b/src/Synthesizer/Plain/Filter/Recursive/Chebyshev.hs @@ -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 diff --git a/src/Synthesizer/Plain/Filter/Recursive/Comb.hs b/src/Synthesizer/Plain/Filter/Recursive/Comb.hs --- a/src/Synthesizer/Plain/Filter/Recursive/Comb.hs +++ b/src/Synthesizer/Plain/Filter/Recursive/Comb.hs @@ -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 diff --git a/src/Synthesizer/Plain/Filter/Recursive/FirstOrder.hs b/src/Synthesizer/Plain/Filter/Recursive/FirstOrder.hs --- a/src/Synthesizer/Plain/Filter/Recursive/FirstOrder.hs +++ b/src/Synthesizer/Plain/Filter/Recursive/FirstOrder.hs @@ -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 diff --git a/src/Synthesizer/Plain/Filter/Recursive/FirstOrderComplex.hs b/src/Synthesizer/Plain/Filter/Recursive/FirstOrderComplex.hs --- a/src/Synthesizer/Plain/Filter/Recursive/FirstOrderComplex.hs +++ b/src/Synthesizer/Plain/Filter/Recursive/FirstOrderComplex.hs @@ -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 diff --git a/src/Synthesizer/Plain/Filter/Recursive/Hilbert.hs b/src/Synthesizer/Plain/Filter/Recursive/Hilbert.hs --- a/src/Synthesizer/Plain/Filter/Recursive/Hilbert.hs +++ b/src/Synthesizer/Plain/Filter/Recursive/Hilbert.hs @@ -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 ((+:), ) diff --git a/src/Synthesizer/Plain/Filter/Recursive/Moog.hs b/src/Synthesizer/Plain/Filter/Recursive/Moog.hs --- a/src/Synthesizer/Plain/Filter/Recursive/Moog.hs +++ b/src/Synthesizer/Plain/Filter/Recursive/Moog.hs @@ -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, ) diff --git a/src/Synthesizer/Plain/Filter/Recursive/SecondOrder.hs b/src/Synthesizer/Plain/Filter/Recursive/SecondOrder.hs --- a/src/Synthesizer/Plain/Filter/Recursive/SecondOrder.hs +++ b/src/Synthesizer/Plain/Filter/Recursive/SecondOrder.hs @@ -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) diff --git a/src/Synthesizer/Plain/Filter/Recursive/SecondOrderCascade.hs b/src/Synthesizer/Plain/Filter/Recursive/SecondOrderCascade.hs --- a/src/Synthesizer/Plain/Filter/Recursive/SecondOrderCascade.hs +++ b/src/Synthesizer/Plain/Filter/Recursive/SecondOrderCascade.hs @@ -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 diff --git a/src/Synthesizer/Plain/Filter/Recursive/Universal.hs b/src/Synthesizer/Plain/Filter/Recursive/Universal.hs --- a/src/Synthesizer/Plain/Filter/Recursive/Universal.hs +++ b/src/Synthesizer/Plain/Filter/Recursive/Universal.hs @@ -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 diff --git a/src/Synthesizer/Plain/IO.hs b/src/Synthesizer/Plain/IO.hs --- a/src/Synthesizer/Plain/IO.hs +++ b/src/Synthesizer/Plain/IO.hs @@ -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 diff --git a/src/Synthesizer/Plain/Interpolation.hs b/src/Synthesizer/Plain/Interpolation.hs --- a/src/Synthesizer/Plain/Interpolation.hs +++ b/src/Synthesizer/Plain/Interpolation.hs @@ -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 diff --git a/src/Synthesizer/Plain/Miscellaneous.hs b/src/Synthesizer/Plain/Miscellaneous.hs --- a/src/Synthesizer/Plain/Miscellaneous.hs +++ b/src/Synthesizer/Plain/Miscellaneous.hs @@ -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) diff --git a/src/Synthesizer/Plain/Noise.hs b/src/Synthesizer/Plain/Noise.hs --- a/src/Synthesizer/Plain/Noise.hs +++ b/src/Synthesizer/Plain/Noise.hs @@ -11,7 +11,6 @@ import Data.List.HT (sliceVertical, ) -import qualified Prelude as P import NumericPrelude.Base import NumericPrelude.Numeric diff --git a/src/Synthesizer/Plain/Oscillator.hs b/src/Synthesizer/Plain/Oscillator.hs --- a/src/Synthesizer/Plain/Oscillator.hs +++ b/src/Synthesizer/Plain/Oscillator.hs @@ -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 diff --git a/src/Synthesizer/Plain/Signal.hs b/src/Synthesizer/Plain/Signal.hs --- a/src/Synthesizer/Plain/Signal.hs +++ b/src/Synthesizer/Plain/Signal.hs @@ -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) -> diff --git a/src/Synthesizer/Plain/ToneModulation.hs b/src/Synthesizer/Plain/ToneModulation.hs --- a/src/Synthesizer/Plain/ToneModulation.hs +++ b/src/Synthesizer/Plain/ToneModulation.hs @@ -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 diff --git a/src/Synthesizer/Plain/Wave.hs b/src/Synthesizer/Plain/Wave.hs --- a/src/Synthesizer/Plain/Wave.hs +++ b/src/Synthesizer/Plain/Wave.hs @@ -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 diff --git a/src/Synthesizer/State/Analysis.hs b/src/Synthesizer/State/Analysis.hs --- a/src/Synthesizer/State/Analysis.hs +++ b/src/Synthesizer/State/Analysis.hs @@ -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 diff --git a/src/Synthesizer/State/Control.hs b/src/Synthesizer/State/Control.hs --- a/src/Synthesizer/State/Control.hs +++ b/src/Synthesizer/State/Control.hs @@ -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 diff --git a/src/Synthesizer/State/Cut.hs b/src/Synthesizer/State/Cut.hs --- a/src/Synthesizer/State/Cut.hs +++ b/src/Synthesizer/State/Cut.hs @@ -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, ) diff --git a/src/Synthesizer/State/Displacement.hs b/src/Synthesizer/State/Displacement.hs --- a/src/Synthesizer/State/Displacement.hs +++ b/src/Synthesizer/State/Displacement.hs @@ -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 diff --git a/src/Synthesizer/State/Filter/NonRecursive.hs b/src/Synthesizer/State/Filter/NonRecursive.hs --- a/src/Synthesizer/State/Filter/NonRecursive.hs +++ b/src/Synthesizer/State/Filter/NonRecursive.hs @@ -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, ) diff --git a/src/Synthesizer/State/Filter/Recursive/Comb.hs b/src/Synthesizer/State/Filter/Recursive/Comb.hs --- a/src/Synthesizer/State/Filter/Recursive/Comb.hs +++ b/src/Synthesizer/State/Filter/Recursive/Comb.hs @@ -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 diff --git a/src/Synthesizer/State/Interpolation.hs b/src/Synthesizer/State/Interpolation.hs --- a/src/Synthesizer/State/Interpolation.hs +++ b/src/Synthesizer/State/Interpolation.hs @@ -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) diff --git a/src/Synthesizer/State/Noise.hs b/src/Synthesizer/State/Noise.hs --- a/src/Synthesizer/State/Noise.hs +++ b/src/Synthesizer/State/Noise.hs @@ -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 diff --git a/src/Synthesizer/State/NoiseCustom.hs b/src/Synthesizer/State/NoiseCustom.hs --- a/src/Synthesizer/State/NoiseCustom.hs +++ b/src/Synthesizer/State/NoiseCustom.hs @@ -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 diff --git a/src/Synthesizer/State/Signal.hs b/src/Synthesizer/State/Signal.hs --- a/src/Synthesizer/State/Signal.hs +++ b/src/Synthesizer/State/Signal.hs @@ -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 diff --git a/src/Synthesizer/Storable/Cut.hs b/src/Synthesizer/Storable/Cut.hs --- a/src/Synthesizer/Storable/Cut.hs +++ b/src/Synthesizer/Storable/Cut.hs @@ -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 diff --git a/src/Synthesizer/Storable/Filter/NonRecursive.hs b/src/Synthesizer/Storable/Filter/NonRecursive.hs --- a/src/Synthesizer/Storable/Filter/NonRecursive.hs +++ b/src/Synthesizer/Storable/Filter/NonRecursive.hs @@ -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 {- | diff --git a/src/Synthesizer/Storable/Oscillator.hs b/src/Synthesizer/Storable/Oscillator.hs --- a/src/Synthesizer/Storable/Oscillator.hs +++ b/src/Synthesizer/Storable/Oscillator.hs @@ -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 diff --git a/src/Test/Main.hs b/src/Test/Main.hs --- a/src/Test/Main.hs +++ b/src/Test/Main.hs @@ -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 : [] diff --git a/src/Test/Sound/Synthesizer/Basic/NumberTheory.hs b/src/Test/Sound/Synthesizer/Basic/NumberTheory.hs new file mode 100644 --- /dev/null +++ b/src/Test/Sound/Synthesizer/Basic/NumberTheory.hs @@ -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)) : + [] diff --git a/src/Test/Sound/Synthesizer/Basic/ToneModulation.hs b/src/Test/Sound/Synthesizer/Basic/ToneModulation.hs --- a/src/Test/Sound/Synthesizer/Basic/ToneModulation.hs +++ b/src/Test/Sound/Synthesizer/Basic/ToneModulation.hs @@ -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 diff --git a/src/Test/Sound/Synthesizer/Causal/Analysis.hs b/src/Test/Sound/Synthesizer/Causal/Analysis.hs new file mode 100644 --- /dev/null +++ b/src/Test/Sound/Synthesizer/Causal/Analysis.hs @@ -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])) : + [] diff --git a/src/Test/Sound/Synthesizer/Generic/Cut.hs b/src/Test/Sound/Synthesizer/Generic/Cut.hs new file mode 100644 --- /dev/null +++ b/src/Test/Sound/Synthesizer/Generic/Cut.hs @@ -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) : + [] diff --git a/src/Test/Sound/Synthesizer/Generic/Filter.hs b/src/Test/Sound/Synthesizer/Generic/Filter.hs new file mode 100644 --- /dev/null +++ b/src/Test/Sound/Synthesizer/Generic/Filter.hs @@ -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)) : + [] diff --git a/src/Test/Sound/Synthesizer/Generic/Fourier.hs b/src/Test/Sound/Synthesizer/Generic/Fourier.hs new file mode 100644 --- /dev/null +++ b/src/Test/Sound/Synthesizer/Generic/Fourier.hs @@ -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) : + [] diff --git a/src/Test/Sound/Synthesizer/Generic/FourierInteger.hs b/src/Test/Sound/Synthesizer/Generic/FourierInteger.hs new file mode 100644 --- /dev/null +++ b/src/Test/Sound/Synthesizer/Generic/FourierInteger.hs @@ -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) : + [] diff --git a/src/Test/Sound/Synthesizer/Generic/Permutation.hs b/src/Test/Sound/Synthesizer/Generic/Permutation.hs new file mode 100644 --- /dev/null +++ b/src/Test/Sound/Synthesizer/Generic/Permutation.hs @@ -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 + -} + [] diff --git a/src/Test/Sound/Synthesizer/Generic/ToneModulation.hs b/src/Test/Sound/Synthesizer/Generic/ToneModulation.hs --- a/src/Test/Sound/Synthesizer/Generic/ToneModulation.hs +++ b/src/Test/Sound/Synthesizer/Generic/ToneModulation.hs @@ -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, ) diff --git a/src/Test/Sound/Synthesizer/Plain/Analysis.hs b/src/Test/Sound/Synthesizer/Plain/Analysis.hs --- a/src/Test/Sound/Synthesizer/Plain/Analysis.hs +++ b/src/Test/Sound/Synthesizer/Plain/Analysis.hs @@ -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 diff --git a/src/Test/Sound/Synthesizer/Plain/Oscillator.hs b/src/Test/Sound/Synthesizer/Plain/Oscillator.hs --- a/src/Test/Sound/Synthesizer/Plain/Oscillator.hs +++ b/src/Test/Sound/Synthesizer/Plain/Oscillator.hs @@ -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 diff --git a/src/Test/Sound/Synthesizer/Plain/ToneModulation.hs b/src/Test/Sound/Synthesizer/Plain/ToneModulation.hs --- a/src/Test/Sound/Synthesizer/Plain/ToneModulation.hs +++ b/src/Test/Sound/Synthesizer/Plain/ToneModulation.hs @@ -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, ) diff --git a/src/Test/Sound/Synthesizer/Plain/Wave.hs b/src/Test/Sound/Synthesizer/Plain/Wave.hs --- a/src/Test/Sound/Synthesizer/Plain/Wave.hs +++ b/src/Test/Sound/Synthesizer/Plain/Wave.hs @@ -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) diff --git a/src/Test/Sound/Synthesizer/Storable/Cut.hs b/src/Test/Sound/Synthesizer/Storable/Cut.hs --- a/src/Test/Sound/Synthesizer/Storable/Cut.hs +++ b/src/Test/Sound/Synthesizer/Storable/Cut.hs @@ -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 () diff --git a/src/Test/Utility.hs b/src/Test/Utility.hs --- a/src/Test/Utility.hs +++ b/src/Test/Utility.hs @@ -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 diff --git a/synthesizer-core.cabal b/synthesizer-core.cabal --- a/synthesizer-core.cabal +++ b/synthesizer-core.cabal @@ -1,5 +1,5 @@ Name: synthesizer-core -Version: 0.4.0.4 +Version: 0.4.1 License: GPL License-File: LICENSE Author: Henning Thielemann @@ -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)