Skip to content

generate coefficients for WFS FIR prefilters #29

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
wants to merge 5 commits into
base: master
Choose a base branch
from
Open

generate coefficients for WFS FIR prefilters #29

Changes from 2 commits
Commits
Show all changes
5 commits
Select commit Hold shift + click to select a range
057c689
Add wfs_prefilter
Dec 11, 2016
763ec97
Fix wfs_prefilter: normalize to unit energy.
Dec 11, 2016
5e17e24
Clean-up wfs_prefilter.
Dec 12, 2016
47e3cc1
Add wrapper function for WFS prefiltering.
Jan 24, 2017
f61f17c
Use DelayedSignal in wfs_25d_fir_prefilter
Feb 11, 2017
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
62 changes: 62 additions & 0 deletions sfs/time/drivingfunction.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -197,3 +197,65 @@ def apply_delays(signal, delays, fs=None):
for channel, cdelay in enumerate(delays_samples):
out[cdelay:cdelay + len(signal), channel] = signal
return out, offset_samples / fs


def wfs_prefilter(dim='2.5D', N=128, fl=50, fu=2000, fs=None, c=None):
"""Get pre-equalization filter for WFS.

Rising slope with 3dB/oct ('2.5D') or 6dB/oct ('2D' and '3D').
Constant magnitude below fl and above fu.
Type 1 linear phase FIR filter of order N.
Simple design via "frequency sampling method".

Parameters
----------
N : int, optional
Filter order, shall be even. For odd N, N+1 is used.
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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

dim : str, optional
Dimensionality, must be '2D', '2.5D' or '3D'.
fl : int, optional
Lower corner frequency in Hertz.
fu : int, optional
Upper corner frequency in Hertz.
(Should be around spatial aliasing limit.)
fs : int, optional
Sampling frequency in Hertz.
c : float, optional
Speed of sound.

Returns
-------
h : (N+1,) numpy.ndarray
Filter taps.
delay : float
Pre-delay in seconds.

"""
N = 2*(int(N + 1)//2) # for odd N, use N+1 instead
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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?

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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?

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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

if fs is None:
fs = defs.fs
if c is None:
c = defs.c

numbins = int(N/2 + 1)
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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?

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

but I think it's not Pythonic.

You are probably right. bins it shall be.

delta_f = fs / (2*numbins - 1)
f = np.arange(numbins) * delta_f
if dim == '2D' or dim == '3D':
alpha = 1
elif dim == '2.5D':
alpha = 0.5

desired = np.power(2 * np.pi * f / c, alpha)
low_shelf = np.power(2 * np.pi * fl / c, alpha)
high_shelf = np.power(2 * np.pi * fu / c, alpha)

l_index = int(np.ceil(fl / delta_f))
u_index = int(min(np.ceil(fu / delta_f), numbins - 1))
desired[:l_index] = low_shelf
desired[u_index:] = min(high_shelf, desired[u_index])
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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?

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Thanks! It is a lot cleaner using clip.


h = np.fft.ifft(np.concatenate((desired, desired[-1:0:-1])))
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Can't you use irfft() here?

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

done.

h = np.roll(np.real(h), numbins - 1)
h = h / np.sqrt(np.sum(abs(h)**2)) # normalize energy
delay = (numbins - 1) / fs
return h, delay