DistortionVisualise.hs

import Graphics.Dynamic.Plot.R2

import Data.Function

main :: IO ()
main = do
   plotWindow $
      [ plotLatest
         [plotMultiple
           [ signalPlot (fmap (*μ) $ initSignal) & legendName "clean"
           , signalPlot (lp12 20 . distortion . fmap (*μ) $ initSignal)
                & legendName "distortion" ]
         | μ <- (/2) . (1.5-) . cos <$> [0, 0.1 ..] ]
      , timeRange, amplitudeRange]
   return ()

type Signal = [Double]

sRate = 441

fuzz :: Signal -> Signal
fuzz = fmap $ (\u -> if abs u<1 then u else signum u) . (*6)

tubeOD :: Signal -> Signal
tubeOD = fmap $ (\u -> tanh u - 0.1/(0.6+u^2)) . (*6)

distortion :: Signal -> Signal
distortion = fuzz . hp6 20

lp6 :: Double -> Signal -> Signal
lp6 ν₀ = go 0
 where go _ [] = []
       go uSvf (u:us) = uSvf : go (uSvf + c₀*(u-uSvf)) us
       c₀ = pi * ν₀ / sRate

hp6 :: Double -> Signal -> Signal
hp6 ν₀ sig = zipWith (-) sig $ lp6 ν₀ sig

lp12 :: Double -> Signal -> Signal
lp12 ν₀ = go (0,0)
 where go _ [] = []
       go (uSvf,uS') (u:us) = uSvf : go (uSvf + uS'*c₀, uS'*c₁ + (u-uSvf)*c₀) us
       c₀ = sqrt $ ν₀ / sRate
       c₁ = 1 - 2*c₀

initSignal :: Signal
initSignal = [ sum [ μ*sin (n*f₀*2*pi*t + φ)
                   | (n,(μ,φ)) <- zip [1..]
                       [(0.6, 0), (0.1, 1), (0.1, 2), (0.1, 1), (0.06, 0), (0.06, 5)]
                       ]
             | t <- timeDomain ]
 where f₀ = 5

timeDomain :: [Double]
timeDomain = [0, 1/sRate .. 2]

signalPlot :: Signal -> DynamicPlottable
signalPlot = lineSegPlot . zip timeDomain 

timeRange = forceXRange (t₀ + l/3, t₁ - l/3)
 where t₀ = minimum timeDomain
       t₁ = maximum timeDomain
       l = t₁ - t₀

amplitudeRange = yInterval (-1, 1)