I have a hard time finding a good "upfirdn"

Specific Domains Signal and Image Processing
1

Hi,

I am translating some code from Matlab to Julia and can’t find an equvalent to upfirdn. The only working option I have is to use PyCall and use the Python (SciPy) version. Like this:

scipy_signal = pyimport("scipy.signal")
s = scipy_signal.upfirdn(filter, s, p, q)

I also tried to make my own version using DSP.resample and DSP.Filters.filt. But it doesn’t work…yet :slight_smile: .

Have I missed some great library?

1 Like
3

It seems like you should be able to implement this easily with resample and filter. What is the issue you are having?

4

I had problems with padding. It works now, but an explanation why it works would be appreciated :slight_smile:

function upfirdn(s, filter, p, q)
    pad = length(s)*p+length(filter)-p
    if p != 1
        s_padded = [s; zeros(ComplexF64, length(filter))] # need padding for the filter
        s = DSP.resample(s_padded, p, [1]) # [1] makes it a zero insertion upsampling
    end

    s = DSP.Filters.filt(FIRFilter(filter), s)
    s = s[1:pad] # Take away padding

    if q != 1
        #result_length = Int(p * length(s))
        s = DSP.resample(s, 1/q, [1])  # downsample
        #s = [s; zeros(ComplexF64, result_length - length(s))] # zero padding to get correct length
    end

    #=
    @show size(s), p, q, size(filter)
    scipy_signal = pyimport("scipy.signal")
    s = scipy_signal.upfirdn(filter, s, p, q)
    @show size(s)
    =#
    return s

end

It works but it is slightly slower than using PyCall. So I need to find a way to speed it up.

5

There are a number of things that you could try, but the first thing would be to

using BenchmarkTools

@benchmark upfirdn(...)

to figure out how much time it’s taking and how much GC it’s doing.

Your best bets are going to be to reduce allocations. perhaps create a buffer that you reuse by calling filt! and then use a @view to strip padding etc.

2 Likes
6

Hi @psteffensen,I just needed to point something out here.
What you have implemented is exactly correct for getting the result, but it is a very slow and inefficient algorithm. You should use the FIRFilter stateful filter object with a rational rate, which uses an efficient polyphase resampler.

FIRFilter(h[, ratio])

Construct a stateful FIRFilter object from the vector of filter taps h . ratio is an optional rational integer which specifies the input to output sample rate relationship (e.g. 147//160 for converting recorded audio from 48 KHz to 44.1 KHz).

https://github.com/JuliaDSP/DSP.jl/blob/7beb55d71231ec1e66abcfe23791553a549a2f44/src/Filters/stream_filt.jl#L148-L155

Best regards,
Tom

2 Likes
7

Thank you Tom (sirtom67) for pointing that out,

I actually ended up using “DSP.Filters.filt!” in my code. I guess it is the same as FIRfilter. I used it because it was much faster, like you say, but I thought I was cheating :blush:

My code only needed “upfir” and not the “down” part, so it ended up looking like this:

function upfirdn(s, filter, p, q)
    ## FIR (Filter)
    result = zeros(ComplexF64, length(s)*p+length(filter))
    DSP.Filters.filt!(result, FIRFilter(filter, p), s)

    return result
end

Then for the “down” part i guess you can you resample to decimate, like this:

    if q != 1
        s = DSP.resample(s, 1/q, [1])
    end
8

Looks good.

The “down” part here suggests part of the problem - for downsampling with factor q > 1, it is inefficient to call filt first, and then throw away (q-1) out of q samples - you’ve done “q”-times as much work as required! So - if it matters (and it may very well not matter if the code runs fast enough due to small signal and / or filter lengths) - it is better to include the “q” downsample factor in the FIRFilter.

1 Like