diff --git a/ChangeLog.md b/ChangeLog.md
new file mode 100644
--- /dev/null
+++ b/ChangeLog.md
@@ -0,0 +1,4 @@
+## 0.1.0.0
+
+Initial version.
+
diff --git a/LICENSE b/LICENSE
new file mode 100644
--- /dev/null
+++ b/LICENSE
@@ -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.
diff --git a/Setup.hs b/Setup.hs
new file mode 100644
--- /dev/null
+++ b/Setup.hs
@@ -0,0 +1,2 @@
+import Distribution.Simple
+main = defaultMain
diff --git a/src/Vocoder.hs b/src/Vocoder.hs
new file mode 100644
--- /dev/null
+++ b/src/Vocoder.hs
@@ -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)
+
diff --git a/src/Vocoder/Filter.hs b/src/Vocoder/Filter.hs
new file mode 100644
--- /dev/null
+++ b/src/Vocoder/Filter.hs
@@ -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))
+
diff --git a/src/Vocoder/Window.hs b/src/Vocoder/Window.hs
new file mode 100644
--- /dev/null
+++ b/src/Vocoder/Window.hs
@@ -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
+
diff --git a/vocoder.cabal b/vocoder.cabal
new file mode 100644
--- /dev/null
+++ b/vocoder.cabal
@@ -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
+