classify-frog-0.2.3: src/SpectralDistributionTest.hs

module Main where

import qualified SpectralDistribution as SD
import qualified SignalProcessing as SP

import qualified Synthesizer.Plain.Filter.Recursive.Universal as UniFilt
import qualified Synthesizer.Plain.Filter.Recursive.FirstOrder as Filt1
import qualified Synthesizer.Plain.Filter.Recursive as FiltRec
import qualified Synthesizer.Generic.Signal as SigG
import qualified Synthesizer.State.Filter.NonRecursive as FiltNR
import qualified Synthesizer.State.Signal as SigS
import qualified Synthesizer.State.Oscillator as Osci
import qualified Synthesizer.State.Displacement as Disp
import qualified Synthesizer.State.Control as Ctrl
import qualified Synthesizer.State.Noise as Noise
import qualified Synthesizer.Causal.Process as Causal

import Synthesizer.Causal.Class (($<), ($*), )

import qualified Data.StorableVector.Lazy as SVL

import Control.Arrow ((&&&), (<<<), (^<<), (<<^), )
import Control.Monad (liftM2, )

import Data.Foldable (forM_, )

import NumericPrelude.Numeric
import NumericPrelude.Base
import Prelude ()


frequency :: Float
frequency = 0.01

duration :: Int
duration = 50000

toneEnvelope :: SigS.T Float
toneEnvelope =
   FiltNR.envelope (Ctrl.exponential 20000 1) $
   Osci.staticSine zero frequency

toneChirp :: SigS.T Float
toneChirp =
   Osci.freqModSine zero $
   Ctrl.exponential 10000 (0.1::Float)

{- |
@sqrt@ asserts that the spectral centroid stays constant.
This works, because our time-domain definition of spectral centroid
computes the quadratic mean of frequencies.
-}
toneBroaden :: SigS.T Float
toneBroaden =
   FiltNR.amplify 0.5 $
   Disp.mix
      (Osci.freqModSine zero $ fmap sqrt $
       Ctrl.line duration (0.01, 0.02::Float))
      (Osci.freqModSine zero $ fmap sqrt $
       Ctrl.line duration (0.01, 0::Float))


toneMix :: Float -> Float -> SigS.T Float
toneMix freq0 freq1 =
   FiltNR.amplify 0.5 $
   Disp.mix
      (Osci.staticSine zero freq0)
      (Osci.staticSine zero freq1)



noiseEnvelope :: SigS.T Float
noiseEnvelope =
   FiltNR.envelope (Ctrl.exponential 20000 1) Noise.white

noiseChirp :: SigS.T Float
noiseChirp =
   (UniFilt.bandpass ^<< UniFilt.causal)
    $< SigS.map
          (UniFilt.parameter . FiltRec.Pole 10)
          (Ctrl.exponential 10000 (0.1::Float))
    $* Noise.white

noiseBroaden :: SigS.T Float
noiseBroaden =
   (UniFilt.bandpass ^<< UniFilt.causal)
    $< SigS.map
          (UniFilt.parameter . flip FiltRec.Pole frequency)
          (Ctrl.exponential 10000 (100::Float))
    $* FiltNR.amplify 0.5 Noise.white



smooth :: Causal.T Float Float
smooth =
   Filt1.lowpass_
   ^<<
   Filt1.causal
   <<<
   Causal.feedConstFst (Filt1.parameter (0.0002::Float))


volume :: SigS.T Float -> Float
volume  =  SigS.sum . SigS.map abs

followEnvelope :: Causal.T Float Float
followEnvelope  =  smooth <<^ abs

spectralDistribution1 :: Causal.T Float (SD.T Float)
spectralDistribution1 =
   (\(d0,(d1,d2)) ->
      SD.mapSpread SD.signedSqrt $ SD.spectralDistribution1 d0 d1 d2)
   ^<<
   (followEnvelope <<< SP.zerothMoment) &&&
   (followEnvelope <<< SP.firstMoment) &&&
   (followEnvelope <<< SP.secondMoment)


volumeSquare :: SigS.T Float -> Float
volumeSquare  =  SigS.sum . SigS.map (^2)

meanSquare :: Causal.T Float Float
meanSquare  =  smooth <<^ (^2)

spectralDistribution2 :: Causal.T Float (SD.T Float)
spectralDistribution2 =
   (\(d0,(d1,d2)) ->
      SD.mapSpread SD.signedSqrt $ SD.spectralDistribution2 d0 d1 d2)
   ^<<
   (meanSquare <<< SP.zerothMoment) &&&
   (meanSquare <<< SP.firstMoment) &&&
   (meanSquare <<< SP.secondMoment)


main :: IO ()
main = do
   forM_ [(0.01, 0.01), (0.01, 0.01*sqrt 2), (0.01, 0.02), (0.01, 0.04)] $
         \(freq0,freq1) ->
      let sig = SigS.take 100000 $ toneMix freq0 freq1
          d0 = volume $ Causal.apply SP.zerothMoment sig
          d1 = volume $ Causal.apply SP.firstMoment  sig
          d2 = volume $ Causal.apply SP.secondMoment sig
          (SD.Cons centroid1 spread1) = SD.spectralDistribution1 d0 d1 d2
          s0 = volumeSquare $ Causal.apply SP.zerothMoment sig
          s1 = volumeSquare $ Causal.apply SP.firstMoment  sig
          s2 = volumeSquare $ Causal.apply SP.secondMoment sig
          (SD.Cons centroid2 spread2) = SD.spectralDistribution2 s0 s1 s2
          r1 = s1/s0
      in  do
          putStrLn $ "\nfreqs: " ++ show (freq0,freq1)
          print (s0, s1, s2, s1/s0, s2/s1)
          print (s2/s0, r1^2, s2/s0-r1^2)
          print
             (sqrt ((freq0^2+freq1^2)/2),
              (centroid1 / (2*pi), sqrt spread1 / (2*pi)),
              (centroid2 / (2*pi), sqrt spread2 / (2*pi)))

   let signals =
          ("tone-env", toneEnvelope) :
          ("tone-chirp", toneChirp) :
          ("tone-broaden", toneBroaden) :
          ("noise-env", noiseEnvelope) :
          ("noise-chirp", noiseChirp) :
          ("noise-broaden", noiseBroaden) :
          []
       write process filename =
          SVL.writeFile filename .
          SigG.fromState SigG.defaultLazySize . SigS.take duration . process
       renderers =
          ("", write id) :
          ("-distribution1", write $ Causal.apply spectralDistribution1) :
          ("-distribution2", write $ Causal.apply spectralDistribution2) :
          []

   sequence_ $
      liftM2
         (\(featureName,render) (signalName,signal) ->
            render ("/tmp/" ++ signalName ++ featureName ++ ".f32") signal)
         renderers signals