packages feed

vocoder (empty) → 0.1.0.0

raw patch · 7 files changed

+532/−0 lines, 7 filesdep +basedep +randomdep +vectorsetup-changed

Dependencies added: base, random, vector, vector-fftw

Files

+ ChangeLog.md view
@@ -0,0 +1,4 @@+## 0.1.0.0++Initial version.+
+ LICENSE view
@@ -0,0 +1,26 @@+Copyright (c) 2020, Marek Materzok+All rights reserved.++Redistribution and use in source and binary forms, with or without+modification, are permitted provided that the following conditions are+met:++1. Redistributions of source code must retain the above copyright+   notice, this list of conditions and the following disclaimer.++2. Redistributions in binary form must reproduce the above copyright+   notice, this list of conditions and the following disclaimer in the+   documentation and/or other materials provided with the+   distribution.++THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS+"AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT+LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR+A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT+OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,+SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT+LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,+DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY+THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT+(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE+OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+ Setup.hs view
@@ -0,0 +1,2 @@+import Distribution.Simple+main = defaultMain
+ src/Vocoder.hs view
@@ -0,0 +1,223 @@+{-| +    Module      : Vocoder+    Description : Phase vocoder+    Copyright   : (c) Celina Pawlińska, 2020+                      Marek Materzok, 2021+    License     : BSD2++This module implements the phase vocoder algorithms. +The implementation is designed to be used directly or to be integrated+into some convenient abstraction (streaming or FRP).+-}+module Vocoder (+      Moduli,+      Phase,+      PhaseInc,+      Frame,+      Window,+      HopSize,+      Length,+      STFTFrame,+      FFTOutput,+      VocoderParams,+      vocoderParams,+      vocFrameLength,+      vocInputFrameLength,+      vocFreqFrameLength,+      vocHopSize,+      vocWindow,+      doFFT,+      doIFFT,+      analysisBlock,+      analysisStep,+      analysisStage,+      synthesisBlock,+      synthesisStep,+      synthesisStage,+      zeroPhase,+      volumeCoeff,+      frameFromComplex,+      frameToComplex,+      addFrames+    ) where++import Data.List+import Data.Complex+import Data.Fixed+import Data.Tuple+import Control.Arrow+import Numeric.FFT.Vector.Invertible as FFT+import Numeric.FFT.Vector.Plan as FFTp+import qualified Data.Vector.Storable as V++-- | Complex moduli of FFT frames. Represent signal amplitudes.+type Moduli = V.Vector Double++-- | Complex arguments of FFT frames. Represent signal phases.+type Phase = V.Vector Double++-- | Phase increments. Represent the deviation of the phase difference+-- between successive frames from the expected difference for the center+-- frequencies of the FFT bins.+type PhaseInc = V.Vector Double++-- | Time domain frame.+type Frame = V.Vector Double++-- | Sampled STFT window function.+type Window = Frame++-- | Offset between successive STFT frames, in samples.+type HopSize = Int++-- | Size in samples.+type Length = Int++-- | STFT processing unit.+type STFTFrame = (Moduli, PhaseInc)++-- | Frequency domain frame.+type FFTOutput = V.Vector (Complex Double)++-- | Type of FFT plans for real signals.+type FFTPlan = FFTp.Plan Double (Complex Double)++-- | Type of IFFT plans for real signals.+type IFFTPlan = FFTp.Plan (Complex Double) Double++-- | Configuration parameters for the phase vocoder algorithm.+data VocoderParams = VocoderParams{+    -- | FFT plan used in analysis stage.+    vocFFTPlan  :: FFTPlan,+    -- | FFT plan used in synthesis stage.+    vocIFFTPlan :: IFFTPlan,+    -- | STFT hop size.+    vocHopSize :: HopSize,+    -- | Window function used during analysis and synthesis.+    vocWindow :: Window+    -- TODO thread safety?+}++-- | FFT frequency frame length.+vocFreqFrameLength :: VocoderParams -> Length+vocFreqFrameLength par = planOutputSize $ vocFFTPlan par++-- | FFT frame length. Can be larger than `vocInputFrameLength` for zero-padding.+vocFrameLength :: VocoderParams -> Length+vocFrameLength par = planInputSize $ vocFFTPlan par++-- | STFT frame length.+vocInputFrameLength :: VocoderParams -> Length+vocInputFrameLength par = V.length $ vocWindow par++-- | Create a vocoder configuration.+vocoderParams :: Length -> HopSize -> Window -> VocoderParams+vocoderParams len hs wnd = VocoderParams (plan dftR2C len) (plan dftC2R len) hs wnd++-- | Apply a window function on a time domain frame.+applyWindow :: Window -> Frame -> Frame+applyWindow = V.zipWith (*)++-- | Change the vector indexing so that the sample at the middle has the number 0.+-- This is done so that the FFT of the window has zero phase, and therefore does not+-- introduce phase shifts in the signal.+rewind :: (V.Storable a) => V.Vector a -> V.Vector a+rewind vec = uncurry (V.++) $ swap $ V.splitAt (V.length vec `div` 2) vec++-- | Zero-pad the signal symmetrically from both sides.+addZeroPadding :: Length+    -> Frame+    -> Frame+addZeroPadding len v+    | diff < 0  = error $ "addZeroPadding: input is " ++ (show diff) ++ " samples longer than target length"+    | diff == 0 = v+    | otherwise = res+    where+    l = V.length v+    diff = len - l+    halfdiff = diff - (diff `div` 2)+    res = (V.++) ((V.++) (V.replicate halfdiff 0) v) (V.replicate (diff-halfdiff) 0)++-- | Perform FFT processing, which includes the actual FFT, rewinding, zero-paddding+-- and windowing.+doFFT :: VocoderParams -> Frame -> FFTOutput+doFFT par =+    FFT.execute (vocFFTPlan par) . rewind . addZeroPadding (vocFrameLength par) . applyWindow (vocWindow par)++-- | Perform analysis on a sequence of frames. This consists of FFT processing+-- and performing analysis on frequency domain frames.+analysisStage :: Traversable t => VocoderParams -> Phase -> t Frame -> (Phase, t STFTFrame)+analysisStage par ph = mapAccumL (analysisBlock par) ph++-- | Perform FFT transform and frequency-domain analysis.+analysisBlock :: VocoderParams -> Phase -> Frame -> (Phase, STFTFrame)+analysisBlock par prev_ph vec = analysisStep (vocHopSize par) (vocFrameLength par) prev_ph (doFFT par vec)++-- | Analyze a frequency domain frame. Phase from a previous frame must be supplied.+-- It returns the phase of the analyzed frame and the result.+analysisStep :: HopSize -> Length -> Phase -> FFTOutput -> (Phase, STFTFrame)+analysisStep h eN prev_ph vec =+    (ph,(mag,ph_inc))+    where+    (mag, ph) = frameFromComplex vec+    ph_inc = V.imap (calcPhaseInc eN h) $ V.zipWith (-) ph prev_ph++-- | Wraps an angle (in radians) to the range [-pi : pi].+wrap :: Double -> Double+wrap e = (e+pi) `mod'` (2*pi) - pi++calcPhaseInc :: Length -> HopSize -> Int -> Double -> Double+calcPhaseInc eN hop k ph_diff =+    (omega + wrap (ph_diff - omega)) / fromIntegral hop+    where+    omega = (2*pi*fromIntegral k*fromIntegral hop) / fromIntegral eN++-- | Perform synthesis on a sequence of frames. This consists of performing+-- synthesis and IFFT processing.+synthesisStage :: Traversable t => VocoderParams -> Phase -> t STFTFrame -> (Phase, t Frame)+synthesisStage par ph frs = mapAccumL (synthesisBlock par) ph frs++-- | Perform frequency-domain synthesis and IFFT transform.+synthesisBlock :: VocoderParams -> Phase -> STFTFrame -> (Phase, Frame)+synthesisBlock par ph fr = (id *** doIFFT par) $ synthesisStep (vocHopSize par) ph fr++-- | Synthesize a frequency domain frame. Phase from the previously synthesized frame+-- must be supplied. It returns the phase of the synthesized frame and the result.+synthesisStep :: HopSize -> Phase -> STFTFrame -> (Phase, FFTOutput)+synthesisStep hop ph (mag, ph_inc) =+    (new_ph, frameToComplex (mag, new_ph))+    where+    new_ph = V.zipWith (+) ph $ V.map (* fromIntegral hop) ph_inc++-- | Perform IFFT processing, which includes the actual IFFT, rewinding, removing padding+-- and windowing.+doIFFT :: VocoderParams -> FFTOutput -> Frame+doIFFT par =+    applyWindow (vocWindow par) . cutCenter (vocInputFrameLength par) . rewind . FFT.execute (vocIFFTPlan par)++-- | Cut the center of a time domain frame, discarding zero padding.+cutCenter :: (V.Storable a) => Length -> V.Vector a -> V.Vector a+cutCenter len vec = V.take len $ V.drop ((V.length vec - len) `div` 2) vec++-- | Zero phase for a given vocoder configuration.+-- Can be used to initialize the synthesis stage.+zeroPhase :: VocoderParams -> Phase+zeroPhase par = V.replicate (vocFreqFrameLength par) 0++-- | An amplitude change coefficient for the processing pipeline.+-- Can be used to ensure that the output has the same volume as the input.+volumeCoeff :: VocoderParams -> Double+volumeCoeff par = fromIntegral (vocHopSize par) / V.sum (V.map (**2) $ vocWindow par)++-- | Converts frame representation to complex numbers.+frameToComplex :: STFTFrame -> FFTOutput+frameToComplex = uncurry $ V.zipWith mkPolar++-- | Converts frame representation to magnitude and phase.+frameFromComplex :: FFTOutput -> STFTFrame+frameFromComplex = V.map magnitude &&& V.map phase++-- | Adds STFT frames.+addFrames :: STFTFrame -> STFTFrame -> STFTFrame+addFrames f1 f2 = frameFromComplex $ V.zipWith (+) (frameToComplex f1) (frameToComplex f2)+
+ src/Vocoder/Filter.hs view
@@ -0,0 +1,158 @@+{-# LANGUAGE TupleSections #-}+{-| +    Module      : Vocoder.Filter+    Description : Frequency-domain filters+    Copyright   : (c) Marek Materzok, 2021+    License     : BSD2++This module defines some useful frequency-domain filters for use in+the vocoder framework.+-}+module Vocoder.Filter (+    FreqStep,+    Filter,+    composeFilters,+    addFilters,+    idFilter,+    amplitudeFilter,+    linearAmplitudeFilter,+    amplify,+    lowpassBrickwall,+    highpassBrickwall,+    bandpassBrickwall,+    bandstopBrickwall,+    lowpassButterworth,+    highpassButterworth,+    bandpassButterworth,+    bandstopButterworth,+    pitchShiftInterpolate,+    convolution,+    convolutionFilter,+    envelope,+    envelopeFilter,+    randomPhaseFilter+    ) where++import Vocoder+import Vocoder.Window+import Control.Monad+import Control.Monad.IO.Class+import System.Random+import qualified Data.Vector.Storable as V++-- | A frequency step is a coefficient relating physical frequency (in Hz)+--   to FFT bin numbers. It is used to define filters independently of the+--   FFT window size.+type FreqStep = Double++-- | The type of frequency-domain filters. A frequency-domain filter is+--   a function transforming STFT frames which can depend on the+--   frequency step.+type Filter m = FreqStep -> STFTFrame -> m STFTFrame++-- | Sequential composition of filters.+composeFilters :: Monad m => Filter m -> Filter m -> Filter m+composeFilters f1 f2 step = f1 step >=> f2 step++-- | Addition of filters.+addFilters :: Monad m => Filter m -> Filter m -> Filter m+addFilters f1 f2 step fr = addFrames <$> f1 step fr <*> f2 step fr++-- | Identity filter.+idFilter :: Monad m => Filter m+idFilter _ = return++-- | Creates a filter which transforms only amplitudes, leaving phase+--   increments unchanged.+amplitudeFilter :: Monad m => (FreqStep -> Moduli -> Moduli) -> Filter m+amplitudeFilter f step (mag, ph_inc) = return (f step mag, ph_inc)++-- | Creates a filter which transforms amplitudes and zeroes the phase +--   increments.+amplitudeFilter0 :: Monad m => (FreqStep -> Moduli -> Moduli) -> Filter m+amplitudeFilter0 f step (mag, ph_inc) = return (f step mag, V.replicate (V.length ph_inc) 0)++-- | Creates a filter which scales amplitudes depending on frequency.+linearAmplitudeFilter :: Monad m => (Double -> Double) -> Filter m+linearAmplitudeFilter f = amplitudeFilter $ \step mag -> V.zipWith (*) mag $ V.generate (V.length mag) $ \k -> f (step * fromIntegral k)++-- | Creates an "amplifier" which scales all frequencies.+amplify :: Monad m => Double -> Filter m+amplify k = linearAmplitudeFilter (const k)++-- | Creates a brickwall lowpass filter.+lowpassBrickwall :: Monad m => Double -> Filter m+lowpassBrickwall t = linearAmplitudeFilter $ \x -> if x <= t then 1.0 else 0.0 ++-- | Creates a brickwall highpass filter.+highpassBrickwall :: Monad m => Double -> Filter m+highpassBrickwall t = linearAmplitudeFilter $ \x -> if x >= t then 1.0 else 0.0 ++-- | Creates a brickwall bandpass filter.+bandpassBrickwall :: Monad m => Double -> Double -> Filter m+bandpassBrickwall t u = linearAmplitudeFilter $ \x -> if x >= t && x <= u then 1.0 else 0.0 ++-- | Creates a brickwall bandstop filter.+bandstopBrickwall :: Monad m => Double -> Double -> Filter m+bandstopBrickwall t u = linearAmplitudeFilter $ \x -> if x <= t || x >= u then 1.0 else 0.0 ++butterworthGain :: Double -> Double -> Double -> Double+butterworthGain n t x = 1 / sqrt (1 + (x / t)**(2 * n))++-- | Creates an n-th degree Butterworth-style lowpass filter.+lowpassButterworth :: Monad m => Double -> Double -> Filter m+lowpassButterworth n t = linearAmplitudeFilter $ butterworthGain n t++-- | Creates an n-th degree Butterworth-style highpass filter.+highpassButterworth :: Monad m => Double -> Double -> Filter m+highpassButterworth n t = linearAmplitudeFilter $ butterworthGain (-n) t++-- | Creates an n-th degree Butterworth-style bandpass filter.+bandpassButterworth :: Monad m => Double -> Double -> Double -> Filter m+bandpassButterworth n t u = linearAmplitudeFilter $ \x -> butterworthGain n u x * butterworthGain (-n) t x++-- | Creates an n-th degree Butterworth-style bandstop filter.+bandstopButterworth :: Monad m => Double -> Double -> Double -> Filter m+bandstopButterworth n t u = linearAmplitudeFilter $ \x -> butterworthGain (-n) t x + butterworthGain n u x++interpolate :: Double -> V.Vector Double -> V.Vector Double+interpolate n v = V.generate (V.length v) f +    where+    f x | i + 1 >= V.length v = 0+        | otherwise = (1-k) * v V.! i + k * v V.! (i+1) where+            x' = n * fromIntegral x+            i = floor x'+            k = x' - fromIntegral i++-- | Creates an interpolative pitch-shifting filter.+pitchShiftInterpolate :: Monad m => Double -> Filter m+pitchShiftInterpolate n _ (mag, ph_inc) = return (interpolate n mag, V.map (/n) $ interpolate n ph_inc)++-- | Convolves the amplitude spectrum using a kernel.+convolution :: V.Vector Double -> Moduli -> Moduli+convolution ker mag = V.generate (V.length mag) $ \k -> V.sum $ flip V.imap ker $ \i v -> v * gmag V.! (i + k) / s+    where+    h = V.length ker `div` 2+    gmag = V.replicate h 0 V.++ mag V.++ V.replicate h 0+    s = V.sum ker++-- | Creates a filter which convolves the spectrum using a kernel.+convolutionFilter :: Monad m => V.Vector Double -> Filter m+convolutionFilter ker = amplitudeFilter0 $ \_ -> convolution ker++-- | Calculates the envelope of an amplitude spectrum using convolution.+envelope :: Length -> Moduli -> Moduli+envelope ksize = V.map ((+(-ee)) . exp) . convolution ker . V.map (log . (+ee))+    where+    ee = 2**(-24)+    ker = if ksize <= 3 then boxWindow ksize else blackmanWindow ksize++-- | Creates a filter which replaces the amplitudes with their envelope.+envelopeFilter :: Monad m => Length -> Filter m+envelopeFilter ksize = amplitudeFilter0 $ \_ -> envelope ksize++-- | Sets the phase increments so that the bins have horizontal consistency.+--   This erases the phase information, introducing "phasiness".+randomPhaseFilter :: MonadIO m => Filter m+randomPhaseFilter _ (mag, ph_inc) = (mag, ) <$> V.replicateM (V.length ph_inc) (randomRIO (0, 2*pi))+
+ src/Vocoder/Window.hs view
@@ -0,0 +1,86 @@+{-| +    Module      : Vocoder.Window+    Description : Window functions+    Copyright   : (c) Marek Materzok, 2021+    License     : BSD2++This module defines popular window functions for use in the vocoder+framework.+-}+module Vocoder.Window (+    makeWindow,+    boxWindow,+    triangleWindow,+    hammingWindow,+    hannWindow,+    generalizedBlackmanWindow,+    blackmanWindow,+    exactBlackmanWindow,+    lanczosWindow,+    flatTopWindow+    ) where++import Vocoder+import qualified Data.Vector.Storable as V++-- | Creates a window of given length by sampling a function+--   on the interval [0,1].+makeWindow :: (Double -> Double) -> Length -> Window+makeWindow f n = V.generate n $ \k -> f (fromIntegral k / (fromIntegral n - 1))++-- | Creates a box window.+boxWindow :: Length -> Window+boxWindow = makeWindow $ const 1++-- | Creates a triangular window.+triangleWindow :: Length -> Window+triangleWindow = makeWindow $ \x -> 2 * (0.5 - abs (x - 0.5))++-- | Creates a Hamming window.+hammingWindow :: Length -> Window+hammingWindow = makeWindow $ \x -> alpha - beta * cos (2 * pi * x)+    where+    alpha = 25.0/46.0+    beta = 21.0/46.0++-- | Creates a Hann window.+hannWindow :: Length -> Window+hannWindow = makeWindow $ \x -> 0.5 * (1 - cos (2 * pi * x))++-- | Creates a generalized Blackman window for a given alpha value.+generalizedBlackmanWindow :: Double -> Length -> Window+generalizedBlackmanWindow a = makeWindow $ \x -> let p = 2 * pi * x in a0 - a1 * cos p + a2 * cos (2 * p)+    where+    a0 = (1 - a) / 2+    a1 = 0.5+    a2 = a / 2++-- | Creates a Blackman window (with alpha=0.16).+blackmanWindow :: Length -> Window+blackmanWindow = generalizedBlackmanWindow 0.16++-- | Creates an exact Blackman window.+exactBlackmanWindow :: Length -> Window+exactBlackmanWindow = makeWindow $ \x -> let p = 2 * pi * x in a0 - a1 * cos p + a2 * cos (2 * p)+    where+    a0 = 7938.0/18608.0+    a1 = 9240.0/18608.0+    a2 = 1430.0/18608.0++-- | Creates a Lanczos window.+lanczosWindow :: Length -> Window+lanczosWindow = makeWindow $ \x -> sinc $ 2 * x - 1+    where+    sinc 0 = 1+    sinc x = sin (pi*x) / (pi*x)++-- | Creates a flat top window.+flatTopWindow :: Length -> Window+flatTopWindow = makeWindow $ \x -> a0 - a1 * cos (2 * pi * x) + a2 * cos (4 * pi * x) - a3 * cos (6 * pi * x) + a4 * cos (8 * pi * x)+    where+    a0 = 0.21557895+    a1 = 0.41663158+    a2 = 0.277263158+    a3 = 0.083578947+    a4 = 0.006947368+
+ vocoder.cabal view
@@ -0,0 +1,33 @@+name:                vocoder+version:             0.1.0.0+homepage:            https://github.com/tilk/vocoder+synopsis:            Phase vocoder+description:+    This package is an implementation of phase vocoder frequency domain+    processing algorithms. It has minimal dependencies on external+    libraries. It can be used directly, but for most uses it's more +    convenient to use a streaming or FRP library wrapper. +    Packages vocoder-conduit and vocoder-dunai are provided for this+    purpose.+license:             BSD2+license-file:        LICENSE+author:              Marek Materzok+maintainer:          tilk@tilk.eu+-- copyright:+category:            Sound+build-type:          Simple+extra-source-files:  ChangeLog.md+cabal-version:       >=1.10++library+  exposed-modules: Vocoder, Vocoder.Window, Vocoder.Filter+  -- other-modules:+  -- other-extensions:+  build-depends:       base >=4.11 && <4.15,+                       vector >= 0.12.1.0 && <0.13,+                       vector-fftw >= 0.1.3.8 && < 0.2,+                       random >= 1.2.0 && < 1.3+  hs-source-dirs:      src+  default-language:    Haskell2010+  ghc-options:         -Wall+