diff --git a/LICENSE b/LICENSE
new file mode 100644
--- /dev/null
+++ b/LICENSE
@@ -0,0 +1,26 @@
+Copyright (c) 2021 Henning Thielemann
+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.
+3. Neither the name of the University nor the names of its contributors
+   may be used to endorse or promote products derived from this software
+   without specific prior written permission.
+
+THIS SOFTWARE IS PROVIDED BY THE REGENTS 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 REGENTS 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/Makefile b/Makefile
new file mode 100644
--- /dev/null
+++ b/Makefile
@@ -0,0 +1,28 @@
+CABALOPTS =
+
+run-test:
+	runhaskell Setup configure --user
+	runhaskell Setup build
+	runhaskell Setup haddock
+
+	runhaskell Setup configure --user -fblas
+	runhaskell Setup build
+
+align-audio.deb:
+	cabal v2-install $(CABALOPTS) \
+	  --installdir=bin --overwrite-policy=always --install-method=copy
+	fpm -s dir -t deb -n align-audio -v 0.0 -p $@ \
+	  -d libsox3 \
+	  -d libfftw3-single3 -d libfftw3-double3 \
+	  --description='Find relative time displacement of two recordings of the same music' \
+	  bin=/usr
+
+align-audio-blas.deb:
+	cabal v2-install -fblas $(CABALOPTS) \
+	  --installdir=bin --overwrite-policy=always --install-method=copy
+	fpm -s dir -t deb -n align-audio -v 0.0 -p $@ \
+	  -d libblas3 -d liblapack3 \
+	  -d libsox3 \
+	  -d libfftw3-single3 -d libfftw3-double3 \
+	  --description='Find relative time displacement of two recordings of the same music' \
+	  bin=/usr
diff --git a/Setup.lhs b/Setup.lhs
new file mode 100644
--- /dev/null
+++ b/Setup.lhs
@@ -0,0 +1,3 @@
+#! /usr/bin/env runhaskell
+> import Distribution.Simple
+> main = defaultMain
diff --git a/align-audio.cabal b/align-audio.cabal
new file mode 100644
--- /dev/null
+++ b/align-audio.cabal
@@ -0,0 +1,73 @@
+Cabal-Version:  2.2
+Name:           align-audio
+Version:        0.0
+License:        BSD-3-Clause
+License-File:   LICENSE
+Author:         Henning Thielemann <haskell@henning-thielemann.de>
+Maintainer:     Henning Thielemann <haskell@henning-thielemann.de>
+Category:       Sound
+Synopsis:       Find relative time displacement of two recordings of the same music
+Description:
+  Say, you have a video with some background music
+  and a clean recording of the background music.
+  You want to know exact displacement of the background music.
+  This program should find it.
+  .
+  > align-audio orig.wav video.wav
+  .
+  The program actually performs a simple correlation.
+  .
+  The program can handle sources of different sample rates
+  if the sample rates are integers.
+Tested-With:    GHC==8.6.5
+Build-Type:     Simple
+Extra-Source-Files:
+  Makefile
+
+Flag blas
+  Description: Build executable using BLAS support
+  Manual: True
+  Default: False
+
+Source-Repository this
+  Tag:         0.0
+  Type:        darcs
+  Location:    http://hub.darcs.net/thielema/align-audio/
+
+Source-Repository head
+  Type:        darcs
+  Location:    http://hub.darcs.net/thielema/align-audio/
+
+Executable align-audio
+  Build-Depends:
+    comfort-fftw >=0.0 && <0.1,
+    comfort-array >=0.5 && <0.6,
+    netlib-ffi >=0.1.1 && <0.2,
+    containers >=0.2 && <0.7,
+    Stream >=0.4.7 && <0.5,
+    storablevector >=0.2 && <0.3,
+    synthesizer-core >=0.7 && <0.9,
+    soxlib >=0.0.1 && <0.1,
+    numeric-prelude >=0.4.1 && <0.5,
+    shell-utility >=0.0 && <0.2,
+    optparse-applicative >=0.11 && <0.17,
+    utility-ht >=0.0.12 && <0.1,
+    base >= 3 && <5
+
+  If flag(blas)
+    Hs-Source-Dirs: src/blas
+    Build-Depends:
+      lapack >=0.4 && <0.5,
+  Else
+    Hs-Source-Dirs: src/plain
+
+  Default-Language: Haskell98
+  GHC-Options:    -Wall
+  Hs-Source-Dirs: src
+  Main-Is:        Main.hs
+  Other-Modules:
+    Correlate
+    CorrelateResample
+    Common
+    Vector
+    Option
diff --git a/src/Common.hs b/src/Common.hs
new file mode 100644
--- /dev/null
+++ b/src/Common.hs
@@ -0,0 +1,19 @@
+module Common where
+
+import qualified Data.Array.Comfort.Storable as Array
+import qualified Data.Array.Comfort.Shape as Shape
+import qualified Data.StorableVector.Lazy as SVL
+import qualified Data.StorableVector as SV
+import Data.Array.Comfort.Storable (Array)
+
+import Foreign.Storable (Storable)
+
+
+cyclicFromVector :: (Storable a) => SVL.Vector a -> Array (Shape.Cyclic Int) a
+cyclicFromVector =
+   Array.mapShape (\(Shape.ZeroBased n) -> Shape.Cyclic n) .
+   Array.fromStorableVector . SV.concat . SVL.chunks
+
+pad :: (Storable a, Num a) => Int -> SVL.Vector a -> SVL.Vector a
+pad n xs =
+   SVL.append xs $ SVL.replicate SVL.defaultChunkSize (n - SVL.length xs) 0
diff --git a/src/Correlate.hs b/src/Correlate.hs
new file mode 100644
--- /dev/null
+++ b/src/Correlate.hs
@@ -0,0 +1,54 @@
+module Correlate where
+
+import Vector (findPeak)
+import Common (cyclicFromVector, pad)
+
+import qualified Numeric.FFTW.Rank1 as Trafo1
+import qualified Numeric.Netlib.Class as Class
+
+import qualified Data.Array.Comfort.Storable as Array
+import qualified Data.Array.Comfort.Shape as Shape
+import qualified Data.StorableVector.Lazy as SVL
+import Data.Array.Comfort.Storable (Array)
+
+import qualified Data.Stream as Stream
+import qualified Data.Complex as Complex
+
+
+{- |
+Round to next higher number of the form m*2^n with m<16.
+This size contains almost only factor 2 and at most one ugly factor < 16
+and we add at most 7% pad data.
+
+>>> map ceilingFFTSize [1..40]
+[1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,18,18,20,20,22,22,24,24,26,26,28,28,30,30,32,32,36,36,36,36,40,40,40,40]
+prop> \n -> n <= ceilingFFTSize n
+prop> \n -> ceilingFFTSize n <= div (n * 107) 100
+-}
+ceilingFFTSize :: Int -> Int
+ceilingFFTSize n =
+   case Stream.dropWhile ((>=16).fst) $
+        Stream.zip
+            (Stream.iterate (\k -> div (k+1) 2) n)
+            (Stream.iterate (2*) 1) of
+      Stream.Cons (m,p) _ -> m*p
+
+
+correlate ::
+   (Class.Real a) =>
+   SVL.Vector a -> SVL.Vector a -> Array (Shape.Cyclic Int) a
+correlate xs ys =
+   let nx = SVL.length xs
+       ny = SVL.length ys
+       n = ceilingFFTSize (nx+ny)
+   in Trafo1.fourierCR $
+      Array.zipWith
+         (\x y -> x * Complex.conjugate y)
+         (Trafo1.fourierRC (cyclicFromVector (pad n xs)))
+         (Trafo1.fourierRC (cyclicFromVector (pad n ys)))
+
+determineLag :: (Class.Real a) => SVL.Vector a -> SVL.Vector a -> Int
+determineLag xs ys =
+   let zs = correlate xs ys
+   in case findPeak zs of
+         (neg,nonNeg) -> if nonNeg >= SVL.length xs then neg else nonNeg
diff --git a/src/CorrelateResample.hs b/src/CorrelateResample.hs
new file mode 100644
--- /dev/null
+++ b/src/CorrelateResample.hs
@@ -0,0 +1,68 @@
+module CorrelateResample where
+
+import Vector (findPeak)
+import Common (cyclicFromVector, pad)
+
+import qualified Numeric.FFTW.Shape as Spectrum
+import qualified Numeric.FFTW.Rank1 as Trafo1
+import qualified Numeric.Netlib.Class as Class
+
+import qualified Data.Array.Comfort.Storable as Array
+import qualified Data.Array.Comfort.Shape as Shape
+import qualified Data.StorableVector.Lazy as SVL
+import Data.Array.Comfort.Storable.Unchecked (Array)
+import Data.Array.Comfort.Shape ((::+)((::+)))
+
+import qualified Data.Complex as Complex
+
+import Algebra.IntegralDomain (divUp)
+
+
+adjust ::
+   (Shape.C sh0, Shape.C sh1, Class.Floating a) =>
+   sh1 -> Array sh0 a -> Array sh1 a
+adjust sh1 xs =
+   let n0 = Shape.size $ Array.shape xs
+       n1 = Shape.size sh1
+   in case compare n0 n1 of
+         EQ -> Array.reshape sh1 xs
+         GT ->
+            Array.takeLeft $ Array.reshape (sh1::+Shape.ZeroBased (n0-n1)) xs
+         LT ->
+            Array.reshape sh1 $ Array.append xs $
+            Array.fromAssociations 0 (Shape.ZeroBased (n1-n0)) []
+
+{- |
+Resample and correlate input vectors.
+We pad the inputs to a multiple of whole seconds.
+We perform resampling implicitly by cutting the spectrum.
+-}
+correlate ::
+   (Class.Real a) =>
+   Integer ->
+   (Integer, SVL.Vector a) -> (Integer, SVL.Vector a) ->
+   Array (Shape.Cyclic Int) a
+correlate dstRate (xrate,xs) (yrate,ys) =
+   let nx = fromIntegral $ SVL.length xs
+       ny = fromIntegral $ SVL.length ys
+       seconds = fromInteger $ divUp (nx*yrate + ny*xrate) (xrate*yrate)
+       paddedAtRate r = seconds * fromInteger r
+       shz = Spectrum.Half $ paddedAtRate dstRate
+   in Trafo1.fourierCR $
+      Array.zipWith
+         (\x y -> x * Complex.conjugate y)
+         (adjust shz $ Trafo1.fourierRC $
+          cyclicFromVector $ pad (paddedAtRate xrate) xs)
+         (adjust shz $ Trafo1.fourierRC $
+          cyclicFromVector $ pad (paddedAtRate yrate) ys)
+
+determineLag ::
+   (Class.Real a) =>
+   Integer -> (Integer, SVL.Vector a) -> (Integer, SVL.Vector a) -> Int
+determineLag dstRate xt@(xrate,xs) yt =
+   let zs = correlate dstRate xt yt
+   in case findPeak zs of
+         (neg,nonNeg) ->
+            if fromIntegral nonNeg * xrate >=
+                  fromIntegral (SVL.length xs) * dstRate
+               then neg else nonNeg
diff --git a/src/Main.hs b/src/Main.hs
new file mode 100644
--- /dev/null
+++ b/src/Main.hs
@@ -0,0 +1,73 @@
+module Main where
+
+import qualified CorrelateResample
+import qualified Correlate
+import qualified Option
+
+import qualified Synthesizer.Basic.Binary as Bin
+import qualified Data.StorableVector.Lazy as SVL
+import qualified Sound.SoxLib as SoxLib
+import Foreign.Storable (peek)
+
+import qualified Options.Applicative as OP
+import Shell.Utility.Exit (exitFailureMsg)
+
+import Text.Printf (printf)
+
+import Control.Monad (when)
+import Control.Applicative (liftA2, (<*>))
+import Data.Maybe.HT (toMaybe)
+import Data.Maybe (fromMaybe)
+
+
+withFirstChannel ::
+   FilePath ->
+   (Double -> SVL.Vector Float -> IO ()) ->
+   IO ()
+withFirstChannel src act =
+   SoxLib.withRead SoxLib.defaultReaderInfo src $ \fmtInPtr -> do
+      fmtIn <- peek fmtInPtr
+      let numChan = fromMaybe 1 $ SoxLib.channels $ SoxLib.signalInfo fmtIn
+      rate <-
+         case SoxLib.rate $ SoxLib.signalInfo fmtIn of
+            Nothing -> exitFailureMsg "no sample rate found"
+            Just rate -> return rate
+      act rate .
+         SVL.sieve numChan . SVL.map Bin.toCanonical =<<
+         SoxLib.readStorableVectorLazy fmtInPtr (SVL.ChunkSize 65536)
+
+maybeInteger :: Double -> Maybe Integer
+maybeInteger nf =
+   let n = round nf
+   in toMaybe (nf == fromInteger n) n
+
+absLag :: Int -> (String, Int)
+absLag n =
+   if n >= 0
+      then ("second", n)
+      else ("first", -n)
+
+
+main :: IO ()
+main = SoxLib.formatWith $ do
+   opt <-
+      OP.execParser $
+      OP.info (OP.helper <*> Option.parseFlags) Option.desc
+   withFirstChannel (Option.input0 opt) $ \rate0 audio0 ->
+      withFirstChannel (Option.input1 opt) $ \rate1 audio1 ->
+      case liftA2 (,) (maybeInteger rate0) (maybeInteger rate1) of
+         Nothing -> do
+            when (rate0/=rate1) $
+               exitFailureMsg "sample rates are fractional and differ"
+            let (ref,lag) = absLag $ Correlate.determineLag audio0 audio1
+            printf
+               "you must delay the %s source by %d samples at %f, i.e. %.6f seconds\n"
+               ref lag rate0 (fromIntegral lag / rate0)
+         Just (irate0,irate1) -> do
+            let dstRate = Option.rate opt
+            let (ref,lag) = absLag $
+                  CorrelateResample.determineLag dstRate
+                     (irate0,audio0) (irate1,audio1)
+            printf
+               "you must delay the %s source by %d samples at %d, i.e. %.3f seconds\n"
+               ref lag dstRate (fromIntegral lag / fromInteger dstRate :: Double)
diff --git a/src/Option.hs b/src/Option.hs
new file mode 100644
--- /dev/null
+++ b/src/Option.hs
@@ -0,0 +1,49 @@
+module Option where
+
+import qualified Shell.Utility.Verbosity as Verbosity
+import Shell.Utility.ParseArgument (parseNumber)
+import Shell.Utility.Verbosity (Verbosity)
+
+import qualified Options.Applicative as OP
+import Options.Applicative (Parser, long, short, help, metavar, value)
+
+import Control.Applicative (pure, (<*>))
+import Data.Monoid ((<>))
+
+
+data T =
+   Cons {
+      verbosity :: Verbosity,
+      rate :: Integer,
+      input0, input1 :: FilePath
+   } deriving (Show)
+
+
+parseFlags :: Parser T
+parseFlags =
+   pure Cons
+   <*> OP.option (OP.eitherReader Verbosity.parse)
+          ( value Verbosity.normal
+         <> short 'v'
+         <> long "verbose"
+         <> metavar "0..3"
+         <> help "verbosity" )
+   <*> OP.option
+         (OP.eitherReader $ parseNumber "sample rate" (0<=) "non-negative")
+          ( value 1000
+         <> short 'r'
+         <> long "rate"
+         <> metavar "RATE"
+         <> help "destination sample rate (default: 1000)" )
+   <*> OP.strArgument (metavar "SRC0")
+   <*> OP.strArgument (metavar "SRC1")
+
+
+desc :: OP.InfoMod a
+desc =
+   OP.fullDesc
+   <>
+   OP.progDesc "Find relative time displacement of two recordings of the same music"
+
+info :: OP.ParserInfo T
+info = OP.info (OP.helper <*> parseFlags) desc
diff --git a/src/blas/Vector.hs b/src/blas/Vector.hs
new file mode 100644
--- /dev/null
+++ b/src/blas/Vector.hs
@@ -0,0 +1,20 @@
+module Vector where
+
+import qualified Numeric.LAPACK.Vector as BlasVector
+import qualified Numeric.Netlib.Class as Class
+import Data.Complex (Complex)
+
+import qualified Data.Array.Comfort.Storable as Array
+import qualified Data.Array.Comfort.Shape as Shape
+import Data.Array.Comfort.Storable (Array)
+
+
+findPeak :: (Class.Real a) => Array (Shape.Cyclic Int) a -> (Int, Int)
+findPeak xs =
+   let (k,_) = BlasVector.argAbs1Maximum xs
+   in (k - Shape.size (Array.shape xs), k)
+
+mulConj ::
+   (Shape.C sh, Eq sh, Class.Real a) =>
+   Array sh (Complex a) -> Array sh (Complex a) -> Array sh (Complex a)
+mulConj = flip BlasVector.mulConj
diff --git a/src/plain/Vector.hs b/src/plain/Vector.hs
new file mode 100644
--- /dev/null
+++ b/src/plain/Vector.hs
@@ -0,0 +1,20 @@
+module Vector where
+
+import qualified Numeric.Netlib.Class as Class
+import qualified Data.Complex as Complex
+import Data.Complex (Complex)
+
+import qualified Data.Array.Comfort.Storable as Array
+import qualified Data.Array.Comfort.Shape as Shape
+import Data.Array.Comfort.Storable (Array)
+
+
+findPeak :: (Class.Real a) => Array (Shape.Cyclic Int) a -> (Int, Int)
+findPeak xs =
+   let (k,_) = Array.argMaximum $ Array.map abs xs
+   in (k - Shape.size (Array.shape xs), k)
+
+mulConj ::
+   (Shape.C sh, Eq sh, Class.Real a) =>
+   Array sh (Complex a) -> Array sh (Complex a) -> Array sh (Complex a)
+mulConj = Array.zipWith (\x y -> x * Complex.conjugate y)