Skip to content
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.

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

add deconv fuction to pyfar/dsp/dsp #18

Draft
wants to merge 10 commits into
base: develop
Choose a base branch
from

Conversation

al-de
Copy link

@al-de al-de commented Apr 5, 2021

This code implements a function named 'deconv' in the DSP module of pyfar.
It returns the deconvolution of a measurement signal and an excitation signal and makes use of the regularized spectrum inversion.
It can be used to calculate an impulse response for example from two sweeps recorded before and after passing a device under test.

@al-de al-de added enhancement New feature or request dsp labels Apr 5, 2021
@al-de al-de requested a review from f-brinkmann April 5, 2021 14:15
@f-brinkmann f-brinkmann requested a review from mberz April 5, 2021 15:45
Copy link

@f-brinkmann f-brinkmann left a comment

Choose a reason for hiding this comment

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

This already looks good to me. Most comments are regarding code style and documentation only.
After @mberz gave feedback we could think about testing.

pyfar/dsp/dsp.py Outdated
Comment on lines 265 to 271
excitation = Signal(np.concatenate((excitation.time,
np.zeros(measurement.n_samples -
excitation.n_samples))),
excitation.sampling_rate,
fft_norm=excitation.fft_norm,
dtype=excitation.dtype,
comment=excitation.comment)

Choose a reason for hiding this comment

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

  • Can be shorter by directly using excitation.time = np.concatenate(...).
  • Does this work for multichannel signals or do we need np.zeros(excitation.cshape + (exc.n_samples - meas.n_samples, )) for that?

Copy link
Author

@al-de al-de Apr 7, 2021

Choose a reason for hiding this comment

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

  1. That's a good point, I implemented the shorter version in my latest commit.
  2. At the moment I think this code would only work for single channel signals, I will use your suggestion to make it work for multichannel signals, too. Thank you. :)

pyfar/dsp/dsp.py Outdated
Comment on lines 274 to 281
# Add Zeros to measurement
measurement = Signal(np.concatenate((measurement.time,
np.zeros(excitation.n_samples -
measurement.n_samples))),
measurement.sampling_rate,
fft_norm=measurement.fft_norm,
dtype=measurement.dtype,
comment=measurement.comment)

Choose a reason for hiding this comment

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

see above :)

Copy link
Author

Choose a reason for hiding this comment

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

check

pyfar/dsp/dsp.py Outdated

Returns
-------
pyfar.Signal

Choose a reason for hiding this comment

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

This will reference the signal class in the Sphinx documentation :)

Suggested change
pyfar.Signal
signal : Signal

You should as well replace pyfar.Signal with Signal above...

Copy link
Author

Choose a reason for hiding this comment

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

check

pyfar/dsp/dsp.py Outdated

# Transform back to time domain and return the impulse resonse
result.fft_norm = None
result.domain = 'time'

Choose a reason for hiding this comment

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

I think this could also be omitted. In case the user needs frequency data, we save two FFTs without this line.

Suggested change
result.domain = 'time'

Copy link
Author

Choose a reason for hiding this comment

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

check

@@ -223,3 +224,75 @@ def spectrogram(signal, dB=True, log_prefix=20, log_reference=1,
times -= times[0]

return frequencies, times, spectrogram


def deconv(measurement, excitation, **kwargs):

Choose a reason for hiding this comment

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

What would be the best name for this?

  • devonv
  • deconvolution
  • spectral_deconvolution

pyfar/dsp/dsp.py Outdated


def deconv(measurement, excitation, **kwargs):
"""Function to calculate the deconvolution of two signals.

Choose a reason for hiding this comment

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

This comment is a bit picky - I think the one line docstring guideline is to start with a verb :)

Suggested change
"""Function to calculate the deconvolution of two signals.
"""Calculate transfer functions by spectral deconvolution of two signals.

Copy link
Author

Choose a reason for hiding this comment

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

check. :)


def deconv(measurement, excitation, **kwargs):
"""Function to calculate the deconvolution of two signals.

Choose a reason for hiding this comment

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

We could add information to make things really obvious for novices:

The transfer function :math:`H(\\omega)` is calculated by spectral
deconvolution (spectral division)

.. math::
   \\H(\\omega) = \\frac{Y(\\omega)}{X(\\omega)}, 

where :math:`X(\\omega)` is the excitation signal and :math:`Y(\\omega)`
the measured signal. Regulated inversion is used to avoid numerical issues
in calculating :math:`\\hat{X(\\omega)} = 1/X(\\omega)` for small values of
:math:`X(\\omega)` (see :py:func:`pyfar.dsp.regulated_spectrum_inversion`).
The transfer function is thus calculated as

.. mat::
    \\H(\\omega) = Y(\\omega)\\hat{X(\\omega)}.

For more information, refer to [#]_

and add below

References
-----------
.. [#] S. Mueller and P. Masserani "Transfer function measurement with sweeps. Directors cut.
       J. Audio Eng. Soc. 49(6):443-471, (2001, June).

Copy link
Author

Choose a reason for hiding this comment

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

Thank you for that suggestion, I added this with my latest commit.

Returns
-------
pyfar.Signal
The resulting signal after deconvolution.

Choose a reason for hiding this comment

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

Mention that you set the fft_norm to 'none'

Copy link
Author

Choose a reason for hiding this comment

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

check

pyfar/dsp/dsp.py Outdated
Comment on lines 289 to 293
if not (measurement.comment is None and excitation.comment is None):
result.comment = (f"IR calculated with deconvolution: [1]"
f"{measurement.comment}; [2]{excitation.comment}")
else:
result.comment = "IR calculated with deconvolution"

Choose a reason for hiding this comment

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

I like the idea with the comment, but would suggest the following to make it a bit more obvious:

result.comment = "Calculated with pyfar.dsp.deconv."
if measurement.comment != 'none':
    result.comment += f" Measured signal: {measurement.comment}."
if excitation.comment != 'none':
    result.comment += f" Excitation signal: {excitation.comment}."

Copy link
Author

Choose a reason for hiding this comment

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

Good idea, this looks much cleaner than in my draft.

----------
measurement : pyfar.Signal
The measurement signal, recorded after passing the device under test.
excitation : pyfar.Signal

Choose a reason for hiding this comment

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

You should add that measurement is zero padded if it is shorter than excitation. (And add the opposite below to excitation)

Copy link
Author

Choose a reason for hiding this comment

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

check

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
dsp enhancement New feature or request
Projects
None yet
Development

Successfully merging this pull request may close these issues.

2 participants