Hi all, I'm sorry if this topic has been rehashed a lot, but I poked around in the mailing list archive, and didn't find what I was looking for. I currently have some free time on my hands, and have been implementing some digital signal processing and spectral/frequency estimation algorithms along with the needed matrix routines in Haskell. For those unfamiliar with this field, most algorthms are defined in textbooks in terms of indexing through discrete sequences. For example the implementation of cross-correlation
rxy x y k | k >= 0 = sum [ x!(i+k) * (conjugate (y!i)) | i <- [0..(n-1-k)] ] | k < 0 = conjugate (rxy y x (-k)) where n = snd (bounds x) + 1
looks very similar to the textbook definition if one uses arrays. A definition using lists would probably use drop, map, and zipWith, and look nothing like the definitions found in the standard texts. Many spectral estimation routines are defined in terms of special matrices (ie, Toeplitz, etc). Arrays defined recursively by list comprehensions make it easy to implement algorithms like Levinson-Durbin recursion, and they look very similar to the mathematical definitions:
levinson r p = (array (1,p) [ (k, a!(p,k)) | k <- [1..p] ], realPart (rho!p)) where a = array ((1,1),(p,p)) [ ((k,i), ak k i) | k <- [1..p], i <- [1..k] ] rho = array (1,p) [ (k, rhok k) | k <- [1..p] ] ak 1 1 = -r!1 / r!0 ak k i | k==i = -(r!k + sum [ a!(k-1,l) * r!(k-l) | l <- [1..(k-1)] ]) / rho!(k-1) | otherwise = a!(k-1,i) + a!(k,k) * (conjugate (a!(k-1,k-i))) rhok 1 = (1 - (abs (a!(1,1)))^2) * r!0 rhok k = (1 - (abs (a!(k,k)))^2) * rho!(k-1)
This could be rewritten for lists, but would probably need to be defined in terms of an aux. recursive function, which destroys the simplicity of the above definition. OK, my question then has to do with the efficiency of lists versus arrays. Do the latest compilers handle handle arrays efficiently, or are lists really the way to go? If there is a performace difference, is it generally big enough to warrant rewriting algorithms? A related question is how is specilization handled in arrays with lazy evaluation. In the definition of levinson above, the a array is defined in terms of the ak function. By doing this, you save some horizontal space, but it also unburdens the programmer from tracking the recursive dependencies. a!(k,k) is needed before a!(i,j) can be calculated, but lazy evaluation takes care of this. If the above function is specialized for r::Array Int (Complex Double) and p::Int, would I be correct to say that the final value of the function would be unboxed, but all intermediate values wouldn't? Now, in some cases, a user may need all of the model orders from 1..p. This is handled easilly enough by just changing the first line. to
levinson r p = (a, fmap realPart rho)
Would the a matrix in the tuple be unboxed with specilization? If anyone is interesting in what I have put together, I will be making everything public sometime next week. I have a lot of algorithms implemented, but I need to clean up the documentation a bit (well, a lot). Thanks. -- Matthew Donadio (m.p.donadio@ieee.org)