{-# OPTIONS -O2 -fno-implicit-prelude #-} {- | Copyright : (c) Henning Thielemann 2006 License : GPL Maintainer : synthesizer@henning-thielemann.de Stability : provisional Portability : requires multi-parameter type classes Avoid importing this module. Better use functions from "Synthesizer.Plain.Oscillator" and "Synthesizer.Basic.Wave" Input data is interpreted as samples of data on a cylinder in the following form: > |* | > | * | > | * | > | * | > | * | > | * | > | * | > | *| > | * | > | * | > | * | > ----------- > * > * > * > * > * > * > * > * > * > * > * > ----------- We have to interpolate in the parallelograms. -} module Synthesizer.Plain.ToneModulation where import qualified Synthesizer.Basic.Phase as Phase import qualified Synthesizer.Plain.Interpolation as Interpolation -- 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.Real as Real import qualified Algebra.Ring as Ring import qualified Algebra.Additive as Additive import qualified Number.NonNegative as NonNeg import qualified Number.NonNegativeChunky as Chunky import Synthesizer.Utility (viewListL, viewListR, clip, mapPair, ) import Control.Monad (guard) import qualified Data.List as List import NumericPrelude.List (replicateMatch, takeMatch, ) import NumericPrelude -- import qualified Prelude as P import PreludeBase -- * general helpers interpolateCell :: Interpolation.T a y -> Interpolation.T b y -> (a, b) -> [[y]] -> y interpolateCell ipLeap ipStep (qLeap,qStep) = Interpolation.func ipStep qStep . map (Interpolation.func ipLeap qLeap) {- | Convert from the (shape,phase) parameter pair to the index within a wave (step) and the index of a wave (leap) in the sampled prototype tone. -} untangleShapePhase :: (Field.C a) => Int -> a -> (a, a) -> (a, a) untangleShapePhase periodInt period (shape,phase) = let leap = shape/period - phase step = shape - leap * fromIntegral periodInt in (leap, step) untangleShapePhaseAnalytic :: (Field.C a) => Int -> a -> (a, a) -> (a, a) untangleShapePhaseAnalytic periodInt period (shape,phase) = let periodRound = fromIntegral periodInt vLeap = (periodRound, periodRound-period) vStep = (1,1) in solveSLE2 (vLeap,vStep) (shape,period*phase) {- Cramer's rule see HTam/Numerics/ZeroFinder/Root, however the matrix is transposed -} solveSLE2 :: Field.C a => ((a,a), (a,a)) -> (a,a) -> (a,a) solveSLE2 a@(a0,a1) b = let det = det2 a in (det2 (b, a1) / det, det2 (a0, b) / det) det2 :: Ring.C a => ((a,a), (a,a)) -> a det2 ((a00,a10),(a01,a11)) = a00*a11 - a10*a01 {- transpose :: ((a,a), (a,a)) -> ((a,a), (a,a)) transpose ((a00,a10),(a01,a11)) = ((a00,a01),(a10,a11)) -} flattenShapePhase :: RealField.C a => Int -> a -> (a, Phase.T a) -> (Int, (a, a)) flattenShapePhase periodInt period (shape,phase) = let (xShape,xWave) = untangleShapePhase periodInt period (shape, Phase.toRepresentative phase) (nLeap,qLeap) = splitFraction xShape (nStep,qStep) = splitFraction xWave {- reverse solveSLE2 for the shape parameter with respect to the rounded (wave,shape) coordinates -} n = nStep + nLeap * periodInt in (n,(qLeap,qStep)) shapeLimits :: Ring.C t => Interpolation.T a v -> Interpolation.T a v -> Int -> t -> (t, t) shapeLimits ipLeap ipStep periodInt len = let minShape = fromIntegral $ interpolationOffset ipLeap ipStep periodInt + periodInt maxShape = minShape + len - fromIntegral (Interpolation.number ipStep + Interpolation.number ipLeap * periodInt) in (minShape, maxShape) interpolationOffset :: Interpolation.T a v -> Interpolation.T a v -> Int -> Int interpolationOffset ipLeap ipStep periodInt = Interpolation.offset ipStep + Interpolation.offset ipLeap * periodInt -- * array based shape variable wave data Prototype a v = Prototype { protoIpLeap, protoIpStep :: Interpolation.T a v, protoIpOffset :: Int, protoPeriod :: a, protoPeriodInt :: Int, protoShapeLimits :: (a,a), protoArray :: Array Int v } makePrototype :: (RealField.C a) => Interpolation.T a v -> Interpolation.T a v -> a -> [v] -> Prototype a v makePrototype ipLeap ipStep period tone = let periodInt = round period ipOffset = interpolationOffset ipLeap ipStep periodInt len = length tone (lower,upper) = shapeLimits ipLeap ipStep periodInt len limits = if lower > upper then error "min>max" else (fromIntegral lower, fromIntegral upper) arr = listArray (0, pred len) tone in Prototype { protoIpLeap = ipLeap, protoIpStep = ipStep, protoIpOffset = ipOffset, protoPeriod = period, protoPeriodInt = periodInt, protoShapeLimits = limits, protoArray = arr } sampledToneCell :: (RealField.C a) => Prototype a v -> a -> Phase.T a -> ((a,a),[[v]]) sampledToneCell p shape phase = let (n, q) = flattenShapePhase (protoPeriodInt p) (protoPeriod p) (uncurry clip (protoShapeLimits p) shape, phase) in (q, map (map (protoArray p ! ) . iterate (protoPeriodInt p +)) $ enumFrom (n - protoIpOffset p)) sampledToneAltCell :: (RealField.C a) => Prototype a v -> a -> Phase.T a -> ((a,a),[[v]]) sampledToneAltCell p shape phase = let (n, q) = flattenShapePhase (protoPeriodInt p) (protoPeriod p) (uncurry clip (protoShapeLimits p) shape, phase) in (q, iterate (drop (protoPeriodInt p)) $ map (protoArray p ! ) (enumFrom (n - protoIpOffset p))) {- M = ((1,1)^T, (periodRound, period-periodRound)^T) equation for the line 0 = (nStep - offset ipStep) + (nLeap - offset ipLeap) * periodInt <(1,periodInt), (offset ipStep, offset ipLeap)> = <(1,periodInt), (nStep,nLeap)> d = <a,x> = <a,M^-1*M*x> = <(M^-T)*a,M*x> = <(M^-T)*a,y> b = (M^-T)*a required: y0 such that y1=0 y0 such that y1=period The line {x : d = <a,x>} converted to (shape,phase) coordinates has constant shape and meets all phases. -} -- * lazy oscillator oscillatorCells :: (RealField.C t) => Interpolation.T t y -> Interpolation.T t y -> t -> [y] -> (t,[t]) -> (Phase.T t,[t]) -> [((t,t),[[y]])] oscillatorCells ipLeap ipStep period sampledTone shapes freqs = let periodInt = round period ptrs = List.transpose $ takeWhile (not . null) $ iterate (drop periodInt) sampledTone ipOffset = interpolationOffset ipLeap ipStep periodInt (skip:skips,coords) = unzip $ oscillatorCoords periodInt period (limitRelativeShapes ipLeap ipStep periodInt sampledTone shapes) freqs in zipWith -- n will be zero within the data, it's only needed for extrapolation (\(k,q) (n,ptr) -> if n>0 then error "ToneModulation.oscillatorCells: limit of shape parameter is buggy" else (q, drop (periodInt+k) ptr)) coords $ tail $ scanl {- since we clip the coordinates before calling oscillatorCells we do not need 'dropRem', since 'drop' would never go beyond the list end -} (\ (n,ptr0) d0 -> dropRem (n+d0) ptr0) (0,ptrs) ((skip - ipOffset - periodInt) : skips) dropFrac :: RealField.C i => i -> [a] -> (Int, i, [a]) dropFrac = let recurse acc n xt = if n>=1 then case xt of _:xs -> recurse (succ acc) (n-1) xs [] -> (acc, n, []) else (acc,n,xt) in recurse 0 dropFrac' :: RealField.C i => i -> [a] -> (Int, i, [a]) dropFrac' = let recurse acc n xt = maybe (acc,n,xt) (recurse (succ acc) (n-1) . snd) (guard (n>=1) >> viewListL xt) in recurse 0 propDropFrac :: (RealField.C i, Eq a) => i -> [a] -> Bool propDropFrac n xs = dropFrac n xs == dropFrac' n xs dropRem :: Int -> [a] -> (Int, [a]) dropRem = let recurse n xt = if n>0 then case xt of _:xs -> recurse (pred n) xs [] -> (n, []) else (n,xt) in recurse dropRem' :: Int -> [a] -> (Int, [a]) dropRem' = let recurse n xt = maybe (n,xt) (recurse (pred n) . snd) (guard (n>0) >> viewListL xt) in recurse propDropRem :: (Eq a) => Int -> [a] -> Bool propDropRem n xs = dropRem n xs == dropRem' n xs {- *Synthesizer.Plain.ToneModulation> Test.QuickCheck.quickCheck (\n xs -> propDropRem n (xs::[Int])) OK, passed 100 tests. *Synthesizer.Plain.ToneModulation> Test.QuickCheck.quickCheck (\n xs -> propDropFrac (n::Rational) (xs::[Int])) OK, passed 100 tests. -} oscillatorCoords :: (RealField.C t) => Int -> t -> (t,[t]) -> (Phase.T t, [t]) -> [(Int,(Int,(t,t)))] oscillatorCoords periodInt period (shape0, shapes) (phase, freqs) = let shapeOffsets = scanl (\(_,s) c -> splitFraction (s+c)) (splitFraction shape0) shapes phases = let (s:ss) = map (\(n,_) -> fromIntegral n / period) shapeOffsets in freqToPhase (Phase.increment (-s) phase) -- phase - s (zipWith (-) freqs ss) in zipWith -- (\(d,s) p -> (d, (s,p))) (\(d,s) p -> (d, flattenShapePhase periodInt period (s,p))) shapeOffsets phases {- mapM print $ take 30 $ let period = 1/0.07::Double in oscillatorCoords (round period) period 0 0 (repeat 0.1) (repeat 0.01) *Synthesizer.Plain.Oscillator> mapM print $ take 30 $ let period = 1/0.07::Rational in oscillatorCoords (round period) period 0 0 (repeat 1) (repeat 0.07) *Synthesizer.Plain.Oscillator> mapM print $ take 30 $ let period = 1/0.07::Rational in oscillatorCoords (round period) period 0 0 (repeat 0.25) (repeat 0.0175) -} -- this function fits better in the Oscillator module {- | Convert a list of phase steps into a list of momentum phases phase is a number in the interval [0,1) freq contains the phase steps -} freqToPhase :: RealField.C a => Phase.T a -> [a] -> [Phase.T a] freqToPhase phase freq = scanl (flip Phase.increment) phase freq limitRelativeShapes :: (RealField.C t) => Interpolation.T t y -> Interpolation.T t y -> Int -> [y] -> (t,[t]) -> (t,[t]) limitRelativeShapes ipLeap ipStep periodInt sampledTone = let -- len = List.genericLength sampledTone len = Chunky.fromChunks (replicateMatch sampledTone one) (minShape, maxShape) = shapeLimits ipLeap ipStep periodInt len fromChunky = NonNeg.toNumber . Chunky.toNumber toChunky = Chunky.fromNumber . NonNeg.fromNumber in mapPair (fromChunky, map fromChunky) . uncurry (limitMaxRelativeValuesNonNeg maxShape) . mapPair (toChunky, map toChunky) . uncurry (limitMinRelativeValues (fromChunky minShape)) {- *Synthesizer.Plain.Oscillator> let ip = Interpolation.linear in limitRelativeShapes ip ip 13 (take 100 $ iterate (1+) (0::Double)) (0::Double, cycle [0.5,1.5]) (13.0,[0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,1.0,0.5,1.5,0.5,1.5,0.5,1.5,0.5,1.5,0.5,1.5,0.5,1.5,0.5,1.5,0.5,1.5,0.5,1.5,0.5,1.5,0.5,1.5,0.5,1.5,0.5,1.5,0.5,1.5,0.5,1.5,0.5,1.5,0.5,1.5,0.5,1.5,0.5,1.5,0.5,1.5,0.5,1.5,0.5,1.5,0.5,1.5,0.5,1.5,0.5,1.5,0.5,1.5,0.5,1.5,0.5,1.5,0.5,1.5,0.5,1.5,0.5,1.5,0.5,1.5,0.5,1.5,0.5,1.5,0.5,1.5,0.5*** Exception: Numeric.NonNegative.Chunky.-: negative number -} limitMinRelativeValues :: (Additive.C a, Ord a) => a -> a -> [a] -> (a, [a]) limitMinRelativeValues xMin x0 xs = let (ys,zs) = span ((<zero).fst) (zip (scanl (+) (x0-xMin) xs) (x0:xs)) in case ys of [] -> (x0,xs) (_:yr) -> (xMin, replicateMatch yr zero ++ case zs of [] -> [] (z:zr) -> fst z : map snd zr) limitMaxRelativeValues :: (Additive.C a, Ord a) => a -> a -> [a] -> (a, [a]) limitMaxRelativeValues xMax x0 xs = let (ys,zs) = span (>zero) (scanl (-) (xMax-x0) xs) in maybe (xMax, replicateMatch xs zero) (\ ~(yl,yr) -> (x0, takeMatch yl xs ++ takeMatch zs (yr : repeat zero))) (viewListR ys) {- | Avoids negative numbers and thus can be used with Chunky numbers. -} limitMaxRelativeValuesNonNeg :: (Additive.C a, Ord a) => a -> a -> [a] -> (a, [a]) limitMaxRelativeValuesNonNeg xMax x0 xs = let (ys,zs) = span fst (scanl (\(_,acc) d -> safeSub acc d) (safeSub xMax x0) xs) in maybe (xMax, replicateMatch xs zero) (\ ~(yl, ~(_,yr)) -> (x0, takeMatch yl xs ++ takeMatch zs (yr : repeat zero))) (viewListR ys) {- *Synthesizer.Plain.Oscillator> limitMaxRelativeValuesNonNeg (let inf = 1+inf in inf) (0::Chunky.T NonNeg.Rational) (repeat 2.5) -} safeSub :: (Additive.C a, Ord a) => a -> a -> (Bool, a) safeSub a b = (a>=b, a-b)