[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