{-# 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 =
=
= <(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 = } 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 (( (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)