Porting C lowpass filter code to Julia, example

Nakul_Tiruviluamala · July 11, 2018, 3:49am

Hi guys,
I am newer to Julia. I am wondering how I could port this code into Julia. In particular, I am having trouble understanding how the for loop can be equivalently executed into Julia. Here is the code off of Julius Smith’s online Stanford signal processing book. Definition of the Simplest Low-Pass

/* C function implementing the simplest lowpass:
 * 
 *      y(n) = x(n) + x(n-1)
 *
 */
double simplp (double *x, double *y, 
               int M, double xm1)
{
  int n;
  y[0] = x[0] + xm1;
  for (n=1; n < M ; n++) {
    y[n] =  x[n]  + x[n-1];
  }
  return x[M-1];
}

(edit: added triple backticks to format code nicer)

Thanks!
Nakul

ssfrr · July 11, 2018, 4:44am

For the most literal translation, you can just shift your indices (Julia generally uses 1-based indices) and use a range for the for loop.

"""
    simplp(x, y, M, xm1)

Simplest lowpass filter, takes `x` as input and puts the output in `y`.
`xm1` is the previous sample of `x` before the first index.

y(n) = x(n) + x(n-1)
"""
function simplp(x, y, M, xm1)
    y[1] = x[1] + xm1
    for n in 2:M
        y[n] = x[n] + x[n-1]
    end
    return x[M]
end

But there are some changes you could make to follow julia conventions:

functions that modify(mutate) their arguments should have a ! at the end of their name
for mutating functions it’s common to put the destination argument first
Julia Arrays know how long they are, so you don’t need to pass the length as an argument
it’s also common to have a non-mutating version without the ! that returns the data rather than filling in a pre-allocated array

So a slightly more julian version would be:

"""
    simplp!(y, x, xm1)

Simplest lowpass filter, takes `x` as input and puts the output in `y`.
`xm1` is the previous sample of `x` before the first index.

y(n) = x(n) + x(n-1)
"""
function simplp!(y, x, xm1)
    length(x) == length(y) || throw(ArgumentError("`x` and `y` must be the same size"))
    length(x) >= 1 || return
    y[1] = x[1] + xm1
    for n in 2:length(x)
        y[n] = x[n] + x[n-1]
    end

    y
end

"""
    simplp(x, xm1)

Simple lowpass filter. Returns the lowpassed signal.
"""
simplp(x, xm1) = simplp!(similar(x), x, xm1)

You can also make it more compact using range indexing on the arrays rather than a for loop

"""
    simplp!(y, x, xm1)

Simplest lowpass filter, takes `x` as input and puts the output in `y`.
`xm1` is the previous sample of `x` before the first index.

y(n) = x(n) + x(n-1)
"""
function simplp!(y, x, xm1)
    length(x) == length(y) || throw(ArgumentError("`x` and `y` must be the same size"))
    length(x) >= 1 || return
    y[1] = x[1] + xm1
    y[2:end] .= x[2:end] .+ x[1:(end-1)]

    y
end

as a side note - your code will be much easier to read if you wrap it in triple-backticks, which will cause discourse to format it like code.

antoine-levitt · July 11, 2018, 7:50am

And for performance, you just do the same with @inbounds in the loop (after appropriate checking) or @views in the y[2:end] .= ...

Nakul_Tiruviluamala · July 11, 2018, 8:01pm

Thank you for these detailed responses. I’ve got a few questions, which are both about programming and about signal processing.

How do you think of xm1? I am having trouble understanding it. Is the fact that it continually gets rewritten why you call this a function with a mutable argument?
You mentioned “it’s also common to have a non-mutating version without the ! that returns the data rather than filling in a pre-allocated array.”

I take it that’s what this does
simplp(x, xm1) = simplp!(similar(x), x, xm1)

Does Julia manage arrays such that one wouldn’t need to declare the y as a function argument?

How would you execute this code in practice if you’ve got an array of numbers you’d like to lowpass?

JackDevine · July 11, 2018, 9:53pm

How do you think of xm1? I am having trouble understanding it. Is the fact that it continually gets rewritten why you call this a function with a mutable argument?

I don’t think that xm1 does get rewritten. xm1 is just the value that you get if you index one index before the start of the array. So I guess that in Julia it should be called x0, since Julia starts indexing at 1. Quoting from the text:

In this implementation, the first instance of x(n-1) is provided as the procedure argument xm1. That way, both x and y can have the same array bounds ( 0,\dots,M-1 ). For convenience, the value of xm1 appropriate for the next call to simplp is returned as the procedure’s value.

If you wanted to get really Julian, then you could use OffsetArrays to deal with the fact that you have to index outside the range of the array.

How would you execute this code in practice if you’ve got an array of numbers you’d like to lowpass?

julia> x = 1:10 |> float;

julia> xm1 = 0
0

julia> simplp(x,xm1)
10-element Array{Float64,1}:
  1.0
  3.0
  5.0
  7.0
  9.0
 11.0
 13.0
 15.0
 17.0
 19.0

julia> y = similar(x);

julia> simplp!(y,x,xm1)
10-element Array{Float64,1}:
  1.0
  3.0
  5.0
  7.0
  9.0
 11.0
 13.0
 15.0
 17.0
 19.0

Nakul_Tiruviluamala · July 12, 2018, 1:13am

Thanks for the response! This works well.

Can you help me understand what the ‘similar’ function does? I would assume that the y array would always be the same size as the x array so I don’t understand why I would have to have ‘y’ as a procedural argument for the user. I am having trouble grasping it because it seems redundant. As a result, the simplp function (sans exclamation point) makes a lot more sense to me. I think that I am also confused as to what a mutable argument is.

Thanks for the help!

JackDevine · July 12, 2018, 1:53am

If you don’t know what a function does, then you can get documentation on it by typing ? and then the function name at the REPL

help?> similar
search: similar

  similar(array, [element_type=eltype(array)], [dims=size(array)])

  Create an uninitialized mutable array with the given element type and size, based upon the given
  source array. The second and third arguments are both optional, defaulting to the given array's eltype
  and size. The dimensions may be specified either as a single tuple argument or as a series of integer
  arguments.

  Custom AbstractArray subtypes may choose which specific array type is best-suited to return for the
  given element type and dimensionality. If they do not specialize this method, the default is an
  Array{element_type}(dims...).

  For example, similar(1:10, 1, 4) returns an uninitialized Array{Int,2} since ranges are neither
  mutable nor support 2 dimensions:

  julia> similar(1:10, 1, 4)
  1×4 Array{Int64,2}:
   4419743872  4374413872  4419743888  0

  Conversely, similar(trues(10,10), 2) returns an uninitialized BitVector with two elements since
  BitArrays are both mutable and can support 1-dimensional arrays:

  julia> similar(trues(10,10), 2)
  2-element BitArray{1}:
   false
   false

  Since BitArrays can only store elements of type Bool, however, if you request a different element type
  it will create a regular Array instead:

  julia> similar(falses(10), Float64, 2, 4)
  2×4 Array{Float64,2}:
   2.18425e-314  2.18425e-314  2.18425e-314  2.18425e-314
   2.18425e-314  2.18425e-314  2.18425e-314  2.18425e-314

  similar(storagetype, indices)

  Create an uninitialized mutable array analogous to that specified by storagetype, but with indices
  specified by the last argument. storagetype might be a type or a function.

  Examples:

  similar(Array{Int}, indices(A))

  creates an array that "acts like" an Array{Int} (and might indeed be backed by one), but which is
  indexed identically to A. If A has conventional indexing, this will be identical to
  Array{Int}(size(A)), but if A has unconventional indexing then the indices of the result will match A.

  similar(BitArray, (indices(A, 2),))

  would create a 1-dimensional logical array whose indices match those of the columns of A.

  similar(dims->zeros(Int, dims), indices(A))

  would create an array of Int, initialized to zero, matching the indices of A.

It is true that it is annoying to make the user have to provide y in order for them to use your function. The reason that you would do so is if you are going to call the simplp function lots of times. If you use the simplp(x, xm1) method, then Julia will have to allocate a y vector every time that you call the function. Therefore, Julia will have to allocate memory every time that you call simplp(x, xm1). This can get very expensive if you are filtering a lot of large vectors. On the other hand, if you use simplp(y, x, xm1), then you can allocate y once and every time that you call the function, Julia will use the same piece of memory to write to y. Julia also gives you some tools for introspecting these sorts of things, for example:

julia> x = 1:10 |> float;

julia> simplp(x, 0)  # Warm up.
10-element Array{Float64,1}:
  1.0
  3.0
  5.0
  7.0
  9.0
 11.0
 13.0
 15.0
 17.0
 19.0

julia> y = similar(x);

julia> simplp!(y, x, 0)  # Warm up.
10-element Array{Float64,1}:
  1.0
  3.0
  5.0
  7.0
  9.0
 11.0
 13.0
 15.0
 17.0
 19.0

julia> @time simplp(x, 0)
  0.000005 seconds (84 allocations: 6.248 KiB)
10-element Array{Float64,1}:
  1.0
  3.0
  5.0
  7.0
  9.0
 11.0
 13.0
 15.0
 17.0
 19.0

julia> @time simplp!(y, x, 0)
  0.000005 seconds (4 allocations: 160 bytes)
10-element Array{Float64,1}:
  1.0
  3.0
  5.0
  7.0
  9.0
 11.0
 13.0
 15.0
 17.0
 19.0

You can see a significant difference in the number of allocations between the two methods. If you decide to get a bit more serious about measuring the performance of the functions that you write, then I would suggest BenchmarkTools.

Nakul_Tiruviluamala · July 13, 2018, 3:32am

I see. That makes a lot more sense. I can now see how performance is increased by declaring the y as well.

Thank you!