[Haskell-cafe] Audio signal processing in Haskell
Henning Thielemann
iakd0 at clusterf.urz.uni-halle.de
Sat Apr 24 15:46:32 EDT 2004
On Sat, 3 Apr 2004, David Roundy wrote:
> On Sat, Apr 03, 2004 at 07:04:18PM +0200, Henning Thielemann wrote:
> >
> > What is the most promising strategy to speed up the computation without
> > loosing the elegance of Haskell's lists?
> > 1. Upgrading GHC from 6.0 to the current version?
> > 2. Use of simplification rules? How can I find out what new rules may be
> > useful? The several dumps that GHC can generate are a bit hard to
> > understand. :-( Is there some more detailed information about what rules
> > GHC already applies than what one can get with -ddump-simpl-stats and
> > -ddump-rules?
> > 3. More strictness? Is a custom list structure with strict entries and
> > maybe with only one constructor (thus always infinite) a considerable
> > alternative? This would disallow the use of lots of functions that deal
> > with lists, including State monad functions.
> > 4. A list of arrays? Then signal processes with feedback like an echo are
> > much more difficult.
>
> Lists are slow, and laziness is slow, but on the other hand, computers are
> fast, so before trying any of 1-4 (all of which seem likely to uglify your
> code a bit)
For 1. a change wouldn't be necessary at all,
for 2. the rules can be kept away from the actual implementation.
Thus these variants would uglify the code only minimally.
> I assume that if your highest priority was writing
> blindingly fast code, you'd be writing in C.
No, not _that_ ugly! :-)
I have already coded some of the algorithms in
Motorola DSP 56000 assembly language.
> I guess one thing you could try (if you haven't already done so) is
> writing a simple program that just reads and writes and does nothing in
> between (except convert to a lazy list), to see how slow that is.
I have attached such an example together with the generated profile.
It generates two signals:
1. 100000 zero samples,
2. 100000 samples of a saw oscillator with exponentially
decaying frequency.
Together it needs 1.0 seconds on my machine with 1.6 GHz which is very
slow compared to the computational effort. The conversion from floating
point numbers to 16 bit signed integers is quite cumbersome, maybe you
know of a better one. Btw. is there a version of 'properFraction' that
evaluates 'properFraction (-1.5)' to (-2,0.5) rather than (-1,-0.5) ?
> Is laziness really an important issue? Are you interested in it mostly
> just for memory consumption reasons, to avoid holding an entire file in
> memory, or are you perhaps thinking to make the code work with real-time
> input?
I don't want to exclude real-time processing from the beginning. It's an
unfortunate connection that laziness is perfectly suited for real-time
processing but is much too slow for it. ;-)
> On the other hand, with ForeignPtr's you leave yourself open to the
> possibility of calling an FFT program through the FFI, and you
> definitely *don't* want to write an FFT in haskell.
It's already done:
http://haskelldsp.sourceforge.net/doc/Numeric.Transform.Fourier.FFT.html
:-)
-------------- next part --------------
Sat Apr 24 13:38 2004 Time and Allocation Profiling Report (Final)
a.out +RTS -p -RTS
total time = 0.98 secs (49 ticks @ 20 ms)
total alloc = 133,367,448 bytes (excludes profiling overheads)
COST CENTRE MODULE %time %alloc
floatToBin Main 49.0 50.1
freqToPhase Main 28.6 29.7
signalToBinaryMono Main 6.1 9.0
osciModSaw Main 6.1 3.6
writeSignalMono Main 4.1 0.0
main Main 4.1 3.6
exponential Main 2.0 3.9
individual inherited
COST CENTRE MODULE no. entries %time %alloc %time %alloc
MAIN MAIN 1 0 0.0 0.0 100.0 100.0
main Main 142 1 0.0 0.0 49.0 42.1
writeSignalMono Main 143 2 4.1 0.0 49.0 42.1
signalToBinaryMono Main 149 1 4.1 4.5 44.9 42.0
floatToBin Main 152 100000 40.8 37.5 40.8 37.5
clip Main 153 100000 0.0 0.0 0.0 0.0
CAF Main 136 11 0.0 0.0 51.0 57.9
main Main 144 0 4.1 3.6 51.0 57.9
exponential Main 154 1 2.0 3.9 2.0 3.9
osciModSaw Main 150 1 6.1 3.6 34.7 33.3
freqToPhase Main 151 1 28.6 29.7 28.6 29.7
writeSignalMono Main 145 0 0.0 0.0 10.2 17.1
signalToBinaryMono Main 146 1 2.0 4.5 10.2 17.1
floatToBin Main 147 100000 8.2 12.6 8.2 12.6
clip Main 148 100000 0.0 0.0 0.0 0.0
CAF GHC.Float 111 9 0.0 0.0 0.0 0.0
CAF GHC.Handle 102 6 0.0 0.0 0.0 0.0
CAF System.Posix.Internals 78 9 0.0 0.0 0.0 0.0
-------------- next part --------------
module Main
where
{-
ghc -prof -auto-all -O -fvia-C SpeedTest.hs
a.out +RTS -p -RTS
-}
import GHC.IOBase
{- saw tooth oscillator with modulated frequency -}
osciModSaw :: RealFrac a => a -> [a] -> [a]
osciModSaw phase freq = map (\x -> 2*x-1) (freqToPhase phase freq)
{- Convert a list of phase steps into a list of momentum phases
phase is a number in the interval [0,1)
freq contains the phase steps -}
freqToPhase :: RealFrac a => a -> [a] -> [a]
freqToPhase phase freq =
scanl (\phase dif -> snd (properFraction (phase+dif))) phase freq
exponential :: Floating a => a -> a -> [a]
exponential halfLife y0 = iterate (*0.5**(1/halfLife)) y0
-- write the signal as binary file containing 1 6bit words
writeSignalMono :: (Floating a, RealFrac a) =>
FilePath -> [a] -> IO ()
writeSignalMono fileName signal =
writeFile fileName (signalToBinaryMono signal)
signalToBinaryMono :: (Floating a, RealFrac a) => [a] -> String
signalToBinaryMono = concat.(map floatToBin)
clip :: Ord a => a -> a -> a -> a
clip lower upper = (max lower).(min upper)
-- work-around the problem, that properFraction doesn't preserve
-- that the fractional part is betten 0 and 1
splitFrac :: (RealFrac a, Integral b) => a -> (b,a)
splitFrac x = if x>=0
then properFraction x
else let (n,f) = properFraction x in (n-1,f+1)
floatToBin :: (Floating a, RealFrac a) => a -> String
floatToBin x = let int = round (32767 * (clip (-1) 1 x))
(hi,lo) = divMod int 256
in [toEnum lo, toEnum (mod hi 256)]
main =
do writeSignalMono "zero.sw" (replicate 100000 0)
writeSignalMono "saw.sw" (take 100000 (osciModSaw 0 (exponential 10000 0.1)))
More information about the Haskell-Cafe
mailing list