generate coefficients for WFS FIR prefilters

[trettberg]

A simple FIR prefilter for WFS, as in the Matlab version.

TODO

• return a signal if Add signal class and helper function #24 gets accepted
• decide on a function name
• think again about default args

[mgeier]

This is only use once below, where it is divided by two and one is added. That's a lot of dividing by 2 and adding one, isn't it?

[trettberg]

It's ugly, but intentional.
The body uses only numbins, but that would be a rather uncommon parameter.

[mgeier]

I wasn't talking about that. But I admit that I wasn't quite clear ...

I was talking about this:

bins = int(2*(int(N + 1)//2)/2 + 1)

This makes me dizzy.

Is this really the canonical way to write this?

[trettberg]

Got rid of it by raising an error for odd N, as you suggest below.

[mgeier]

I don't like "num" prefixes. What about just calling it bins?

[trettberg]

Hmm. bins sounds like a list of bins to me.
I suppose you would not like nbins either?

[mgeier]

Indeed, nbins is as bad as numbins. This is of course very much a matter of taste, but I think it's not Pythonic.
Do you know any Python code of the standard library (or in code of any of the Python core contributors) that uses this naming convention?

I understand that bins sounds like a list of bins, but it should be immediately clear from its usage that it's not.
I personally wouldn't see this as a problem.

If you want to make it abundantly clear that it is a quantity rather than a list of things, you should call it something like filtersize. Or simply size, since it's the only "size" thats relevant in this context.
What do you think about this?

[trettberg]

but I think it's not Pythonic.

You are probably right. bins it shall be.

[mgeier]

Can't you use irfft() here?

[trettberg]

done.

[mgeier]

Can desired even be smaller here? Isn't min() superfluous?

Also, do you really need to calculate both *_shelf and *_index?

Couldn't you just clip desired to low_shelf and high_shelf (without needing the indices)?

Or, alternatively, use desired[l_index] and desired[u_index] without calculating low_shelf and high_shelf?

[trettberg]

Thanks! It is a lot cleaner using clip.

[mgeier]

Why don't you just raise a ValueError if the user provides wrong input?

[mgeier]

Thanks, but this isn't a filter!
It is a function that creates filter coefficients, therefore the name doesn't seem to fit.

Anyway, it would be interesting to see how this is supposed to be used.
Without usage, I cannot tell if the API (or the name, for that matter) makes sense.

[trettberg]

Do you mean a function called "filter" should actually apply a filter?

For the moment, it would be applied manually by the user.

[chris-hld]

In general, I would expect any function filter to apply a filter, too. Is there a benefit in returning just coefficients?
And I think the prefiltering will be called a lot, it's probably easier to get back the filtered signals or otherwise provide a function that does the actual filtering (maybe for all driving signals at once).
Something like d , t_filterdelay = wfs_prefilter(d, ...) would be perfect.

[trettberg]

Then I guess get_wfs_fir_prefilter might be a better name?

@mgeier : Could you point me to a better way to convolve N (or 1) signals with N (or 1) filters?
I could only come up with something like this:

def fftconvolve_1d(a, v, axis=1):
    """1-d Convolution along axis. 
    """
    if a.ndim > 2 or v.ndim > 2:
        raise ValueError('dim mismatch');
    output_length = a.shape[axis] + v.shape[axis] -1
    nfft = int(2**np.ceil(np.log2(output_length)))
    A = np.fft.rfft(a, nfft, axis)
    V = np.fft.rfft(v, nfft, axis)
    out = np.fft.irfft(A * V, nfft, axis)
    return np.delete(out, np.s_[output_length:], axis)

[spors]

Since the filter is specific to 2.5-dimensional synthesis I would expect 25D in the function name.

[trettberg]

Since the filter is specific to 2.5-dimensional synthesis I would expect 25D in the function name.

@spors : The docs use the term "pre-equalization filter" also for the 6dB high pass in the 2D and 3D case. Is this wrong?

[spors]

get_wfs_25d_fir_prefilter would be in line with the naming scheme of the driving functions we used so far.

[mgeier]

Yes, IMHO a function called "filter" should either be a filter or create and return a filter.
A set of FIR coefficients is not a filter.

In a Pythonic naming scheme, the word "get" should be avoided as much as possible. Most of the times, it's utterly meaningless. Unless the "getting" is really one of the main tasks of a function, e.g. getattr().

re FFT convolution:
This is a bit of a sad story in the scientific Python universe. I've made a few notes about this in the communication acoustics exercises.

As long as it stays 1D, scipy.signal.fftconvolve() should be the weapon of choice. However, for 2D arrays, this does a two-dimensional convolution, which is normally not what we want.
It would be nice to have this as a column-wise convolution (similar to fftfilt() in GNU Octave), but this doesn't exist yet in NumPy/SciPy. There is an open issue: scipy/scipy#1364 (last activity summer 2014).

[trettberg]

Then what about wfs_prefilter_fir or wfs_fir_prefilter_coefficients?
We can use wrappers if dimensionality should be part of the function name.

Since this only about filter coefficients, no assumptions about usage are necessary.
I've changed the PR title accordingly.

@mgeier: Thanks for clarifying the FFT convolution situation!
I agree scipy.signal.fftconvolve is fine for the job here.
However this is will be an issue e.g. for fraction-of-sample interpolation.

[mgeier]

@trettberg

Since this only about filter coefficients, no assumptions about usage are necessary.

This is not true.
You are already assuming that the user will call a function to get FIR coefficients, right?
But is this really necessary?
Is it really the best user interface?
Or should the function actually return a filter function?
Or something completely different?
I think those questions can be best discussed based on the suggested usage ...

[mgeier]

@trettberg

However this is will be an issue e.g. for fraction-of-sample interpolation.

I don't understand, can you please elaborate?

[trettberg]

You are already assuming that the user will call a function to get FIR coefficients, right?

No. It is merely assumed that filter coefficients are needed for FIR filtering.

Personally, I'd use them directly:

signal , _ = wfs_prefiter_fir()
d = driving_signals(..., signal, ...)

But I do not suggest this or any particular usage.

However this is will be an issue e.g. for fraction-of-sample interpolation.

I don't understand, can you please elaborate?

Sorry, I should not have brought it up as it does not belong here.
For fractional delay interpolation with FIR filters, some kind of "columnwise convolution" is necessary.
But we don't need to discuss this now.

[mgeier]

No. It is merely assumed that filter coefficients are needed for FIR filtering.

Does that mean that wfs_prefilter_fir() is supposed to be an "internal" function that should not be called by users?

Your usage above looks like a "shortcut".
In the end we want to be able to also plot sound fields generated by arbitrary signals, not only impulse responses, right?

I think we actually should suggest a certain usage, and this usage should include an explicit filtering step.

Initially, I thought we could just use generic filter functions everywhere, which could be FIR or IIR filters or whatever nonlinear stuff.
To get an impulse response, I imagined to just feed those functions with the trivial signal [1] (assuming a linear system).
Sadly, this doesn't work, because the impulse response could be infinite ...
So we either have to cut the output signal to the same length as the input (which would mean that we have to use the much less convenient [1, 0, 0, 0, 0, ...] to get a meaningful impulse response), or we somehow have to specify how much longer the output should be (which wouldn't be very convenient either).

Or should we limit ourselves to FIR filters?

[trettberg]

I think we actually should suggest a certain usage, and this usage should include an explicit filtering step.

Initially, I thought we could just use generic filter functions everywhere, which could be FIR or IIR filters or whatever nonlinear stuff.

By "generic filter function" you mean a function that takes a signal and yields a signal?

Then the FIR filter might be used like this:

def wfs_fir_prefilter(sig, ...):
    h = _wfs_prefilter_fir(...)
    return fftconvolve(h, sig)

re IR length:
I don't see why this would be different from the source functions:
The filter function decides how many samples it returns.

[mgeier]

By "generic filter function" you mean a function that takes a signal and yields a signal?

Exactement!

And the question about the API would be if we should have the impulse-reponse-generating function as a user-callable function of just as an implementation detail (or not at all).

I don't see why this would be different from the source functions:

AFAICT, the source functions don't have this kind of freedom. They get a grid and have to use exactly those points.

What's the similarity there?

The filter function decides how many samples it returns.

That sounds reasonable.
I think it is better than the users having to zero-pad their input signals.
But what would be a meaningful length to return for an IIR filter (in case we are going to have those at some point)?

[trettberg]

And the question about the API would be if we should have the impulse-reponse-generating function as a user-callable function of just as an implementation detail (or not at all).

Implementation detail is fine I think. As long as we have only 2.5D WFS it does not have to be a separate function. But it doesn't hurt either.

the source functions don't have this kind of freedom. [...]
What's the similarity there?

Sorry, that was a mistake. Somehow I forgot that source functions return the field at a single time instant only.

But what would be a meaningful length to return for an IIR filter?

Decay below some threshold maybe?

[trettberg]

I've tried to pick up the open ends:

• Coefficient generation is done internally in _wfs_prefilter_fir. The user calls the API function wfs_25d_fir_prefilter, see example below.

• Since the filter design is pretty much "inspired" by the Matlab version, the same default parameters are used.
  This is of course entirely arbitrary, but IMHO arbitrary default parameters are better than none.

Any objections? Anything else to do?

import numpy as np
import sfs

grid = sfs.util.xyz_grid([-0.5, 3.5], [-2, 2], 0, spacing=0.01)
x0, n0, a0 = sfs.array.linear(21, 0.1)
delays, weights = sfs.time.drivingfunction.wfs_25d_point(x0, n0, [-1, 0, 0])

sig = [1]

# apply pre-equalization & individual delays/weights
sig, t_offset_filter = sfs.time.drivingfunction.wfs_25d_fir_prefilter(sig)
d, t_offset_delays = sfs.time.drivingfunction.driving_signals(delays, weights, sig)

t = 0.01
t -= t_offset_delays + t_offset_filter

p = sfs.time.soundfield.p_array(x0, d, a0, t, grid)
sfs.plot.level(p, grid, vmin=-80, vmax =-40, cmap='YlOrRd')
sfs.plot.loudspeaker_2d(x0, n0)

[chris-hld]

I just tried it, too. Works pretty nice!

As a little remak: I was wondering, if it's worth shifting the call of wfs_25d_fir_prefilter to drivingfunction.wfs_25d_*.
I think this would be more consistent with the theory...

[trettberg]

I was wondering, if it's worth shifting the call of wfs_25d_fir_prefilter to drivingfunction.wfs_25d_*.

We could plug a filter function into drivingfunction.driving_signals().
However, I'd wait until we have more a bit more functionality (e.g. 2D and 3D, possibly a more sophisticated delay-handling). It will be easier to see "what goes where".

[mgeier]

FYI: overlap-add convolution has been implemented in SciPy: scipy/scipy#10869