hommage-0.0.5: Sound/Hommage/DFTFilter.hs
-- This file is the new version for C-code FFT algorithm with Nyquist allready on position 1.
-- | This module imports the C-functions for a Fast Fourier Transform and encapsulates them in
-- list-mapping functions. (At first it was a DFT algoritm which I changed, but I didn'd change
-- all the names of the functions, so it is called dft...)
-- The functionality of this module is used for filtering audio data, i. e. for building
-- audio filters.
module Sound.Hommage.DFTFilter
(
-- | Filtering is done with an FFT algorithm written in C (in the file dft.c).
-- This module imports the c-functions for fft ('analyseDFT') and inverse fft ('syntheseDFT'),
-- which operate on Arrays with 1024 Double values.
-- The array-operations are embedded in a mechanism that does buffering, overlapping and windowing
-- and works like a black box that maps one input value to one output value in the IO monad
-- ('mkAnalyse', 'mkSynthese', 'mkFilterBuffered').
-- Around these IO's some list mapping functions are wrapped ('dftanalyse', 'dftsynthese', 'dftfilter',
-- 'dftfilterBy') for convenient use in a Haskell programm.
--
-- Mapping the Fourier coefficients during the filter process means filtering the data.
-- The result depends on the way the coefficients are mapped resp. modified while mapping.
-- The function to do this has the type of 'CoeffMap'.
--
-- /Description of the filtering process:/
--
-- Filtering is done by mapping the coefficients with 'CoeffMap'.
-- For this the wave data is split into parts of 1024 values, where every part overlaps with the
-- preceeding and succeeding part by 512 values. Via the fourier-transform (analyse) every part
-- is mapped to its frequency spectrum, i. e. its fourier coefficients.
--
-- The 'CoeffMap' action is applied to every frequency spectrum (see 'CoeffArr') in sequence.
-- It reads the coefficients from the first 'CoeffArr' where the input is stored and writes it
-- (or some data derived from it) to the second one. This is the point where the \'real\' filtering
-- happens. Via the inverse fourier-transform (synthese) the resulting coefficients are mapped
-- to their real signal values. These signal values are multiplied with a half cosinus curve for
-- fading in and out (windowing). Then they are mixed in the way their source data was split:
-- Every part overlaps for 512 values with the last and the next one.
-- * User-Level Filter Functions
dftanalyse
, dftsynthese
, dftfilter
, dftfilterBy
, dftfilterBy'
-- * CoeffArr
, CoeffArr
, readCoeffArr
-- , readCoeffArr'
, writeCoeffArr
-- , writeCoeffArr'
, storeCoeff
, unstoreCoeff
-- * CoeffMap
, CoeffMap
, mkCoeffMap
, coeffmap
-- * IO Wrapper
, mkAnalyse
, mkSynthese
, mkFilterBuffered
-- * Interface to C-code
, analyseDFT
, syntheseDFT
, kurveDFT
)
where
import GHC.Weak
import GHC.IOBase
import Foreign.ForeignPtr
import Foreign.Ptr
import Foreign.Storable
import Data.Array.Storable
import Data.Complex
import Data.Array.IO
import Sound.Hommage.Misc
---------------------------------------------------------------------------------------------------
-- ([Complex Double] -> [Complex Double]) -> IO CoeffMap
-- | (Fast) Fourier Transformation and Inverse with a buffer-mapping action for filtering.
dftfilterBy :: IO (CoeffMap ()) -> [Double] -> [Double]
dftfilterBy cm = drop 1536 . inList (mkFilterBuffered cm) . (++ replicate 512 0.0)
-- | (Fast) Fourier Transformation and Inverse with a buffer-mapping action for filtering.
dftfilterBy' :: s -> IO (s -> CoeffMap s) -> [Double] -> [Double]
dftfilterBy' s cm = drop 1536 . inList (mkFilterBuffered' s cm) . (++ replicate 512 0.0)
-- | (Fast) Fourier Transformation and Inverse with an additional argument with sublists of
-- factors to weight the coefficients. See 'mkCoeffMap' for a decription of that argument.
dftfilter :: [Double] -> [Double] -> [Double]
dftfilter ds = dftfilterBy (mkCoeffMap (replicate 512 0 ++ ds))
-- | Inverse (Fast) Fourier Transformation
dftsynthese :: [Complex Double] -> [Double]
dftsynthese = inList mkSynthese
-- | (Fast) Fourier Transformation.
dftanalyse :: [Double] -> [(Complex Double)]
dftanalyse = inList mkAnalyse
---------------------------------------------------------------------------------------------------
foreign import ccall "cinit" initfou :: Int
foreign import ccall "canalyse" cAnalyse :: Ptr Double -> Ptr Double -> IO ()
foreign import ccall "csynthese" cSynthese :: Ptr Double -> Ptr Double -> IO ()
foreign import ccall "ckurve" cKurve :: Ptr Double -> IO ()
---------------------------------------------------------------------------------------------------
-- | Maps an array with 1024 real values to its frequency spectrum with 512 complex values.
-- The spectrum is stored in the second array with 1024 Doubles which are the real and imaginary
-- parts of the complex values. See 'CoeffArr' for a detailed frequency spectrum
-- desription.
analyseDFT :: StorableArray Int Double -> StorableArray Int Double -> IO ()
analyseDFT wave coeffs = seq initfou $
withStorableArray wave $ \ptr_w ->
withStorableArray coeffs $ \ptr_c ->
cAnalyse ptr_w ptr_c
-- | Maps a frequency spectrum to its signal value. Both arrays have a length of 1024; the first one
-- contains the 512 complex values of the spectrum. The real signal value will be stored in
-- the second array. See 'CoeffArr' for a detailed frequency spectrum
-- desription.
syntheseDFT :: StorableArray Int Double -> StorableArray Int Double -> IO ()
syntheseDFT coeffs wave = seq initfou $
withStorableArray wave $ \ptr_w ->
withStorableArray coeffs $ \ptr_c ->
cSynthese ptr_c ptr_w
-- | Fades a signal which is stored in an array of length 1024 in and out.
-- The signal is multiplied with a cosinus curve with range 0..pi.
kurveDFT :: StorableArray Int Double -> IO ()
kurveDFT wave = seq initfou $
withStorableArray wave $ \ptr_w ->
cKurve ptr_w
---------------------------------------------------------------------------------------------------
--
---------------------------------------------------------------------------------------------------
-- | Represents a frequency spectrum.
-- Range is from 0 to 1023 (N=1024).
-- Index 0 is constant value, index 1 is Nyquest frequency (real part of the N\/2 frequency).
-- Index 2 and 3 are the complex value for the basefrequency, index 4 and 5
-- are the complex value for the double basefrequency, 6 and 7 for 3 * basefrequency and so on:
--
-- > arr [0] = real part of zero frequency (imaginary part is allways 0)
-- > arr [1] = real part of Nyquest frequency (N/2) (imaginary part is allways 0)
-- > arr [2] = real part of base frequency
-- > arr [3] = imag part of base frequency
-- > arr [4] = real part of base frequency * 2
-- > arr [5] = imag part of base frequency * 2
-- > arr [6] = real part of base frequency * 3
-- > arr [7] = imag part of base frequency * 3
-- > ...
-- > arr [N-4] = real part of base frequency * (N/2 - 2)
-- > arr [N-3] = imag part of base frequency * (N/2 - 2)
-- > arr [N-2] = real part of base frequency * (N/2 - 1)
-- > arr [N-1] = imag part of base frequency * (N/2 - 1)
--
-- Pseudocode:
--
-- > c (0) = (arr [0] :+ 0)
-- > c (512) = (arr [1] :+ 0)
--
-- and for all i=[1..511]:
--
-- > c (i) = (arr [i*2] :+ arr [i*2+1])
--
type CoeffArr = StorableArray Int Double
{-
-- | Reads a Complex value of a CoeffArr as described at 'CoeffArr'.
readCoeffArr :: CoeffArr -> Int -> IO (Complex Double)
readCoeffArr arr n = if n == 0 then readArray arr 0 >>= \x -> readArray arr 1023 >>= \y -> return (x :+ y)
else let j = 2*n
j' = j-1
in readArray arr j' >>= \x -> readArray arr j >>= \y -> return (x :+ y)
-- | Writes a Complex value of a CoeffArr as described at 'CoeffArr'.
writeCoeffArr :: CoeffArr -> Int -> Complex Double -> IO ()
writeCoeffArr arr n c = if n == 0 then writeArray arr 0 (realPart c) >> writeArray arr 1023 (imagPart c)
else let j = 2*n
j' = j-1
in writeArray arr j' (realPart c) >> writeArray arr j (imagPart c)
-}
-- | Reads a Complex value of a CoeffArr (with Nyquest as imaginary part of coeff 0).
readCoeffArr :: CoeffArr -> Int -> IO (Complex Double)
readCoeffArr arr n =
let j = 2*n
j' = j+1
in readArray arr j >>= \x -> readArray arr j' >>= \y -> return (x :+ y)
-- | Writes a Complex value of a CoeffArr (with Nyquest as imaginary part of coeff 0).
writeCoeffArr :: CoeffArr -> Int -> Complex Double -> IO ()
writeCoeffArr arr n c =
let j = 2*n
j' = j+1
in writeArray arr j (realPart c) >> writeArray arr j' (imagPart c)
storeCoeff :: IOArray Int (Complex Double) -> CoeffArr -> IO ()
storeCoeff arr brr = do
for 0 (<512) (+1) $ \i -> let i1 = i * 2
i2 = i1 + 1
in do (x :+ y) <- readArray arr i
writeArray brr i1 x
writeArray brr i2 y
unstoreCoeff :: CoeffArr -> IOArray Int (Complex Double) -> IO ()
unstoreCoeff brr arr = do
for 0 (<512) (+1) $ \i -> let i1 = i * 2
i2 = i1 + 1
in do x <- readArray brr i1
y <- readArray brr i2
writeArray arr i (x :+ y)
---------------------------------------------------------------------------------------------------
-- | An action that maps the fourier coefficients from the first array to the second.
-- Modifying the data while mapping means filtering. See 'CoeffArr' for a description of the arrays.
type CoeffMap a = CoeffArr -> CoeffArr -> IO a
-- | Constructs a Coeffmap.
-- The sublists contain the 512 real values with which the basefrequency and the 511 complex fourier coefficients
-- are multiplied (the Nyquist frequency will be zero).
-- The first value is the factor for the constant coefficient, the second for the base frequency, the next one
-- for the double base freq and so on.
-- If the sublists have more than 512 elements, these elements are thrown away.
-- If it is shorter, the array will be filled up with zeros.
{-
mkCoeffMap :: [[Double]] -> IO (CoeffMap ())
mkCoeffMap cs = do
r <- newIORef cs
b <- newArray (0,511) 0.0
return $ \a1 a2 -> let loop i (x:xs) = if i >= 512 then return () else
writeArray b i x >> loop (i+1) xs
loop i [] = for i (<512) (+1) (\i -> writeArray b i 0.0)
in do xs <- readIORef r
let (xi,xr) | null xs = ([], [])
| otherwise = (head xs, tail xs)
writeIORef r xr
loop 0 xi
coeffmap b a1 a2
-}
mkCoeffMap :: [Double] -> IO (CoeffMap ())
mkCoeffMap cs = do
r <- newIORef cs
b <- newArray (0,511) 0.0
return $ \a1 a2 -> let loop i (x:xs) = if i >= 512
then writeIORef r (x:xs)
else writeArray b i x >> loop (i+1) xs
loop i [] = for i (<512) (+1) (\i -> writeArray b i 0.0)
>> writeIORef r []
in do xs <- readIORef r
loop 0 xs
coeffmap b a1 a2
-- | 'coeffmap' takes an Array with 512 values and uses them for a weigted map of the fourier
-- coefficients. The Nyquist frequency is muliplied with 0.0.
coeffmap :: StorableArray Int Double -> CoeffMap ()
coeffmap filt coeffs coeffs' = do
for 0 (<512) (+1) (\i -> let j = 2*i
j' = j + 1
in
readArray filt i >>= \c ->
readArray coeffs j >>= \c1 ->
readArray coeffs j' >>= \c2 ->
writeArray coeffs' j (c * c1) >>
writeArray coeffs' j' (c * c2) )
writeArray coeffs' 1 0.0
{-
coeffmap' :: ([Complex Double] -> [Complex Double]) -> IO (CoeffMap ())
coeffmap' f = do
let toDest arr 1024 _ = return ()
toDest arr n (x:xs) = writeArray arr n x >> toDest arr (n+1) xs
toDest arr n [] = writeArray arr n 0 >> toDest arr (n+1) []
return $ \arr brr -> getElems arr >>= toDest brr 0 . f
-}
---------------------------------------------------------------------------------------------------
-- | Constructs an action that maps wave-data to coefficient-data.
-- Has a delay of 512, i. e. the first 512 elements are zero and the result has (these)
-- 512 elemets more than the input.
mkAnalyse :: IO (Maybe Double -> IO (Maybe (Complex Double)))
mkAnalyse = do
warr <- newArray (0, 1023) 0.0
carr <- newArray (0, 1023) 0.0
oarr <- newArray (0, 511) (0.0 :+ 0.0)
-- SIGNALSTREAM reada closea <- openSignalStream sa
rpos <- newIORef 0
rcnt <- newIORef Nothing
let read ma = checkpos >> readIORef rcnt >>= maybe (more ma) rest
more ma = maybe irest snext ma
rest k | k <= 0 = return Nothing
| otherwise = do writeIORef rcnt $ Just (k-1)
pos <- readIORef rpos
x <- readArray oarr pos
writeIORef rpos (pos + 1)
return $ Just x
irest = do writeIORef rcnt (Just 511)
pos <- readIORef rpos
fillrest pos
x <- readArray oarr pos
writeIORef rpos (pos + 1)
return $ Just x
snext a = do pos <- readIORef rpos
x <- readArray oarr pos
writeArray warr (pos+512) a
writeIORef rpos (pos + 1)
return $ Just x
checkpos = readIORef rpos >>= \p -> if p < 512 then return () else do
writeIORef rpos 0
analyseDFT warr carr
unstoreCoeff carr oarr
for 0 (<512) (+1) (\i -> readArray warr (i+512) >>= writeArray warr i)
fillrest k = for (512+k) (<1024) (+1) (\i -> writeArray warr i 0.0)
return read
---------------------------------------------------------------------------------------------------
-- | Constructs an action that maps coefficient-data to wave-data.
-- Has a delay of 512, i. e. the first 512 elements are zero and the result has (these)
-- 512 elemets more than the input.
mkSynthese :: IO (Maybe (Complex Double) -> IO (Maybe Double))
mkSynthese = do
iarr <- newArray (0, 511) (0.0 :+ 0.0)
carr <- newArray (0, 1023) 0.0
warr <- newArray (0, 1023) 0.0
(oarr :: IOArray Int Double) <- newArray (0, 511) 0.0
-- SIGNALSTREAM reada closea <- openSignalStream sa
rpos <- newIORef 0
rcnt <- newIORef Nothing
let read ma = checkpos >> readIORef rcnt >>= maybe (more ma) rest
more ma = maybe irest snext ma
rest k | k <= 0 = return Nothing
| otherwise = do writeIORef rcnt $ Just (k-1)
pos <- readIORef rpos
x <- readArray oarr pos
writeIORef rpos (pos + 1)
return $ Just x
irest = do writeIORef rcnt (Just 511)
pos <- readIORef rpos
fillrest pos
x <- readArray oarr pos
writeIORef rpos (pos + 1)
return $ Just x
snext a = do pos <- readIORef rpos
x <- readArray oarr pos
writeArray iarr pos a
writeIORef rpos (pos + 1)
return $ Just x
checkpos = readIORef rpos >>= \p -> if p < 512 then return () else do
writeIORef rpos 0
storeCoeff iarr carr
for 0 (<512) (+1) (\i -> readArray warr (i+512) >>= writeArray oarr i)
syntheseDFT carr warr
kurveDFT warr
for 0 (<512) (+1) (\i -> readArray warr i >>= \x ->
readArray oarr i >>= \y -> writeArray oarr i (x+y))
fillrest k = for k (<512) (+1) (\i -> writeArray iarr i (0.0 :+ 0.0))
return read
---------------------------------------------------------------------------------------------------
-- | Constructs an action that maps wave-data to wave-data via a Fast Foutrier Transform and inverse,
-- filtered by the given 'CoeffMap'.
-- Has a delay of 1024.
mkFilterBuffered :: IO (CoeffMap ()) -> IO (Maybe Double -> IO (Maybe Double))
mkFilterBuffered mkf = do
f <- mkf
warr <- newArray (0, 1023) 0.0
carr <- newArray (0, 1023) 0.0
carr' <- newArray (0, 1023) 0.0
warr' <- newArray (0, 1023) 0.0
(oarr :: IOArray Int Double) <- newArray (0, 511) 0.0
-- SIGNALSTREAM reada closea <- openSignalStream sa
rpos <- newIORef 0
rcnt <- newIORef Nothing
let read ma = checkpos >> readIORef rcnt >>= maybe (more ma) rest
more ma = maybe irest snext ma
rest k | k <= 0 = return Nothing
| otherwise = do writeIORef rcnt $ Just (k-1)
pos <- readIORef rpos
x <- readArray oarr pos
writeIORef rpos (pos + 1)
return $ Just x
irest = do writeIORef rcnt (Just 511)
pos <- readIORef rpos
fillrest pos
x <- readArray oarr pos
writeIORef rpos (pos + 1)
return $ Just x
snext a = do pos <- readIORef rpos
x <- readArray oarr pos
writeArray warr (pos+512) a
writeIORef rpos (pos + 1)
return $ Just x
checkpos = readIORef rpos >>= \p -> if p < 512 then return () else do
writeIORef rpos 0
analyseDFT warr carr
for 0 (<512) (+1) (\i -> readArray warr (i+512) >>= writeArray warr i >>
readArray warr' (i+512) >>= writeArray oarr i)
f carr carr'
syntheseDFT carr' warr'
kurveDFT warr'
for 0 (<512) (+1) (\i -> readArray warr' i >>= \x ->
readArray oarr i >>= \y -> writeArray oarr i (x+y))
fillrest k = for (512+k) (<1024) (+1) (\i -> writeArray warr i 0.0)
return read
---------------------------------------------------------------------------------------------------
-- | Constructs an action that maps wave-data to wave-data via a Fast Foutrier Transform and inverse,
-- filtered by the given 'CoeffMap'.
-- Has a delay of 1024.
mkFilterBuffered' :: s -> IO (s -> CoeffMap s) -> IO (Maybe Double -> IO (Maybe Double))
mkFilterBuffered' st mkf = do
stref <- newIORef st
f <- mkf
warr <- newArray (0, 1023) 0.0
carr <- newArray (0, 1023) 0.0
carr' <- newArray (0, 1023) 0.0
warr' <- newArray (0, 1023) 0.0
(oarr :: IOArray Int Double) <- newArray (0, 511) 0.0
-- SIGNALSTREAM reada closea <- openSignalStream sa
rpos <- newIORef 0
rcnt <- newIORef Nothing
let read ma = checkpos >> readIORef rcnt >>= maybe (more ma) rest
more ma = maybe irest snext ma
rest k | k <= 0 = return Nothing
| otherwise = do writeIORef rcnt $ Just (k-1)
pos <- readIORef rpos
x <- readArray oarr pos
writeIORef rpos (pos + 1)
return $ Just x
irest = do writeIORef rcnt (Just 511)
pos <- readIORef rpos
fillrest pos
x <- readArray oarr pos
writeIORef rpos (pos + 1)
return $ Just x
snext a = do pos <- readIORef rpos
x <- readArray oarr pos
writeArray warr (pos+512) a
writeIORef rpos (pos + 1)
return $ Just x
checkpos = readIORef rpos >>= \p -> if p < 512 then return () else do
writeIORef rpos 0
analyseDFT warr carr
for 0 (<512) (+1) (\i -> readArray warr (i+512) >>= writeArray warr i >>
readArray warr' (i+512) >>= writeArray oarr i)
s <- readIORef stref
s' <- f s carr carr'
writeIORef stref s'
syntheseDFT carr' warr'
kurveDFT warr'
for 0 (<512) (+1) (\i -> readArray warr' i >>= \x ->
readArray oarr i >>= \y -> writeArray oarr i (x+y))
fillrest k = for (512+k) (<1024) (+1) (\i -> writeArray warr i 0.0)
return read
---------------------------------------------------------------------------------------------------
{-
-- | Reads a Complex value of a CoeffArr
readCoeffArr :: CoeffArr -> Int -> IO (Complex Double)
readCoeffArr arr n = if n == 0 then readArray arr 0 >>= \x -> readArray arr 1023 >>= \y -> return (x :+ y)
else let j = 2*n
j' = j-1
in readArray arr j' >>= \x -> readArray arr j >>= \y -> return (x :+ y)
-- | Writes a Complex value of a CoeffArr
writeCoeffArr :: CoeffArr -> Int -> Complex Double -> IO ()
writeCoeffArr arr n c = if n == 0 then writeArray arr 0 (realPart c) >> writeArray arr 1023 (imagPart c)
else let j = 2*n
j' = j-1
in writeArray arr j' (realPart c) >> writeArray arr j (imagPart c)
-- | Copies coefficient data from an IOArray with Complex values (512 elements) to a
-- StorableArray with Doubles (1024 elements)
storeCoeff :: IOArray Int (Complex Double) -> CoeffArr -> IO ()
storeCoeff arr brr = do
for 1 (<512) (+1) $ \i -> let i2 = i * 2
i1 = i2 - 1
in do (x :+ y) <- readArray arr i
writeArray brr i1 x
writeArray brr i2 y
(x :+ y) <- readArray arr 0
writeArray brr 0 x
writeArray brr 1023 y
-- | Copies coefficient data from a StorableArray with Doubles (1024 elements) to
-- an IOArray with Complex values (512 elements).
unstoreCoeff :: CoeffArr -> IOArray Int (Complex Double) -> IO ()
unstoreCoeff brr arr = do
for 1 (<512) (+1) $ \i -> let i2 = i * 2
i1 = i2 - 1
in do x <- readArray brr i1
y <- readArray brr i2
writeArray arr i (x :+ y)
x <- readArray brr 0
y <- readArray brr 1023
writeArray arr 0 (x :+ y)
coeffmap :: StorableArray Int Double -> CoeffMap ()
coeffmap filt coeffs coeffs' = do
for 1 (<512) (+1) (\i -> let j' = 2*i
j = j' - 1
in
readArray filt i >>= \c ->
readArray coeffs j >>= \c1 ->
readArray coeffs j' >>= \c2 ->
writeArray coeffs' j (c * c1) >>
writeArray coeffs' j' (c * c2) )
readArray filt 0 >>= \c -> readArray coeffs 0 >>= \c' -> writeArray coeffs' 0 (c * c')
writeArray coeffs' 1023 0.0
-}