packages feed

align-audio (empty) → 0.0

raw patch · 11 files changed

+433/−0 lines, 11 filesdep +Streamdep +basedep +comfort-arraysetup-changed

Dependencies added: Stream, base, comfort-array, comfort-fftw, containers, lapack, netlib-ffi, numeric-prelude, optparse-applicative, shell-utility, soxlib, storablevector, synthesizer-core, utility-ht

Files

+ LICENSE view
@@ -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.
+ Makefile view
@@ -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
+ Setup.lhs view
@@ -0,0 +1,3 @@+#! /usr/bin/env runhaskell+> import Distribution.Simple+> main = defaultMain
+ align-audio.cabal view
@@ -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
+ src/Common.hs view
@@ -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
+ src/Correlate.hs view
@@ -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
+ src/CorrelateResample.hs view
@@ -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
+ src/Main.hs view
@@ -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)
+ src/Option.hs view
@@ -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
+ src/blas/Vector.hs view
@@ -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
+ src/plain/Vector.hs view
@@ -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)