Dependent Iteration

Consider the following code:

// Return RC low-pass filter output samples, given input samples,
// time interval dt, and time constant RC
 function lowpass(real[0..n] x, real dt, real RC)
   var real[0..n] y
   var real alpha := dt / (RC + dt)
       y[0] := x[0]
   for i from 1 to n
       y[i] := alpha * x[i] + (1-alpha) * y[i-1]
   
   return y

This is from the section "Digital Simulation" from the Wikipedia article "Low-pass filters" and it has an element that strikes me as fundamental in the "for loops considered so 20th century" mindset. Each iteration uses the results of the previous iteration.

The example code is a digital simulation of a physical phenomenon, a simple resistor-capacitor (RC) circuit which allows frequencies below the cutoff and blocks frequencies above the cutoff. The y[i-1] is helping model the physical phenomenon of the first derivative of output voltage which the capacitor introduces.

Similar physical reasons, like neighboring physical objects and derivatives, make it very natural to use the procedural paradigm where access to neighboring data and previous calculation results involve array access. The whole area of physical simulation programming, including digital signal processing, is centered around the procedural programming paradigm, and pseudocode and examples are largely in terms of x[i], y[i-1] notation. And when they write for or do loops, they might often be for a million, ten million, or more elements. Speed in everything matters on those kinds of input orders of magnitude.

In light of these requirements, can this type of programming be done in a functional style? Should it?

I helped a customer fix a performance problem in some code very similar to this in the "current iteration depends on previous iteration's results" sense. Access to his output data structure was slowing things down and we solved the problem by using a local variable to remember the previous iteration's output value. But I started wondering what a good general iterating idiom would be for this type of problem: one that lets the dependent iteration be expressed elegantly, but one that also performs well.

A directly recursive version might look like this:

fun lowpass(x,dt,RC) =
  let
     val n = length x
     val alpha = dt/(RC+dt)
     fun lowpass' i yprev acc = 
       if i >= n then
         acc
       else
         let val yi = alpha*sub(x,i)+(1-alpha)*yprev
         in
           lowpass' (i+1) yi (append(acc,yi)
         end
  in
    lowpass' 0 (sub(x,0)) []
  end

With a good native compiler this should perform well, the tail recursion optimized to an assembly jump and the lowpass' function parameters assigned to registers. But is this as intuitive? It certainly looks longer than the procedural version.

Of course this is also a natural for fold, with the previous output value included in the accumulator spot. But I wonder if it's less intuitive than the explicitly recursive version, and I wonder if the performance can possibly come close to the metal lanugages unless the compiler can inline fold's function argument and reconstruct the implicit loop.

I think a list-comprehension version would be more compact, and provide both readability and performance. The Mathematica version using Table is reasonable:

lowPass[x_, dt_?NumberQ, RC_] := 
 Module[{yprev = First@x, alpha = dt/(RC + dt)},
  Table[
    With[{ycurrent = alpha*x[[i]] + (1 - alpha)*yprev},
      yprev = ycurrent;
      ycurrent],
    {i, 1, Length[x]}
  ]
 ]

But then this is not purely functional: it uses a stateful local variable yprev. The fold option may be the best purely functional one, but it still seems like a hybrid has the best chance for high performance while not needing the for/do loop token overhead.