Frequency domain sweep synthesis #199
Conversation
There was a problem hiding this comment.
I'd suggest a different structure by splitting up the functionality:
- add domain option (
'time'or'frequency') tolinear_sweepandexponential_sweep, explanation with ripples in respective other domain - add
perfect_sweepfunction, this would avoid to have a lot of unused parameters - add function to generate sweep with given magnitude spectrum
I think this structure is a little more user-friendly. Of course the functions could share private methods needed for the sweep generation in the background.
|
Yes, this would have the benefit of making the functions more explicit. The disadvantage is that time and frequency domain synthesis require different parameters, thus making the functions less explicit:
That's maybe something we can discuss during the next weekly. I can image both... |
pyfar/signals/deterministic.py
Outdated
| return signal, group_delay | ||
|
|
||
|
|
||
| def magnitude_weighted_sweep( |
There was a problem hiding this comment.
The more I think about this, I'd prefer magnitude_spectrum_weighted_sweep:
| def magnitude_weighted_sweep( | |
| def magnitude_spectrum_weighted_sweep( |
pyfar/signals/deterministic.py
Outdated
| magnitude : Signal, string | ||
| Specify the magnitude response of the sweep. | ||
|
|
||
| signal | ||
| The magnitude response as :py:class:`~pyfar.classes.audio.Signal` | ||
| object. If ``signal.n_samples`` is smaller than `n_samples`, zeros | ||
| are padded to the end of `signal`. Note that `frequency_range` is | ||
| not required in this case. | ||
| ``'linear'`` | ||
| Design a sweep with linearly increasing frequency and a constant | ||
| magnitude spectrum. | ||
| ``'exponential'`` | ||
| Design a sweep with exponentially increasing frequency. The | ||
| magnitude decreases by 3 dB per frequency doubling and has constant | ||
| energy in fiters of relative constant bandwidth (e.g. octaves). | ||
| ``'perfect'`` | ||
| Perfect sweep according to [#]_. Note that the parameters | ||
| `start_margin`, `stop_margin`, `frequency_range` and `double` are | ||
| not required in this case. |
There was a problem hiding this comment.
Calling this magnitude is a bit misleading. The effect on the magnitude is only a secondary effect.
There was a problem hiding this comment.
Renamed it to sweep_type. But it should not be critical any more since it was moved to a private function
Co-authored-by: Marco Berzborn <mberz@users.noreply.github.com>
…/pyfar/pyfar into feature/spectral_sweep_synthesis
pyfar/signals/deterministic.py
Outdated
| The linear sweep can also be generated in the frequency domain (see | ||
| :py:func:`~linear_sweep_freq`). Time domain synthesis exhibits a constant | ||
| temporal envelope in trade of slight ripples in the magnitude response. | ||
| Frequency domain synthesis exhibits smooth magnitude spectra and in trade |
There was a problem hiding this comment.
| Frequency domain synthesis exhibits smooth magnitude spectra and in trade | |
| Frequency domain synthesis exhibits smooth magnitude spectra in trade |
| fade_in : int | ||
| Duration of a squared sine fade-in in samples. The fade starts at the | ||
| first sample of the sweep that is closer than 60 dB to the absolute | ||
| maximum of the sweep time signal. The default ``0`` does not apply |
| fade_out : int | ||
| Duration of a squared cosine fade-out in samples. The fade ends at the | ||
| last sample of the sweep that is closer than 60 dB to the absolute | ||
| maximum of the sweep time signal. The default ``0`` does not apply |
pyfar/signals/deterministic.py
Outdated
| The order of the Butterworth filters that are applied to limit the | ||
| frequency range by a high-pass if ``frequency_range[0] > 0`` and/or by | ||
| a low-pass if ``frequency_range[1] < sampling_rate / 2``. The default |
There was a problem hiding this comment.
I'd suggest to move this information to frequency_range, where the Butterworth filter is already mentioned.
There was a problem hiding this comment.
good idea. Done as suggested
pyfar/signals/deterministic.py
Outdated
| group_delay_sweep : FrequencyData | ||
| The group delay of the sweep in seconds as a single sided spectrum. |
| >>> magnitude = pf.dsp.filter.low_shelve( | ||
| ... pf.signals.impulse(2**16), 500, 20, 2) | ||
| >>> magnitude = pf.dsp.filter.butterworth(magnitude, 8, 50, 'highpass') |
There was a problem hiding this comment.
What about just using a simple bandpass instead? Then no one would need to think about how the low shelve filter looks like (as I did :D)
pyfar/signals/deterministic.py
Outdated
| ------- | ||
| sweep : Signal | ||
| The sweep signal. The Signal is in the time domain, has a maximum | ||
| absolute amplitude of 1 and the ``rms`` FFT normalization |
There was a problem hiding this comment.
I think rms normalization is not correct here
pyfar/signals/deterministic.py
Outdated
| ------- | ||
| sweep : Signal | ||
| The sweep signal. The Signal is in the time domain, has a maximum | ||
| absolute amplitude of 1 and the ``rms`` FFT normalization |
| >>> auto_correlation = np.empty(2**8) | ||
| >>> for idx, shift in enumerate(range(-2**7, 2**7)): | ||
| >>> auto_correlation[idx] = np.dot( | ||
| ... sweep.time.flatten(), | ||
| ... np.roll(sweep.time.flatten(), shift)) | ||
| >>> | ||
| >>> auto_correlation /= pf.dsp.energy(sweep) |
There was a problem hiding this comment.
If I remember correctly, alternative methods to calculate the autocorrelation are the convolution with the time-flipped (real-valued) signal, or via the inverse fourier transform of the squared magnitude spectrum. Maybe these approaches lead to shorter code using pyfar.dsp functionality?
There was a problem hiding this comment.
Indeed, this can be done in a more compact way and especially without a loop. For example using https://docs.scipy.org/doc/scipy/reference/generated/scipy.signal.correlate.html
There was a problem hiding this comment.
I tried this, but its not cyclic (regardless of mode and method) and causes ringing in the auto-correlation:
import pyfar as pf
import numpy as np
import scipy.signal as sgn
import matplotlib.pyplot as plt
%matplotlib inline
sweep = pf.signals.linear_perfect_sweep(2**8)
# compute auto-correlation
auto_correlation = sgn.correlate(
sweep.time.flatten(), sweep.time.flatten())
auto_correlation /= pf.dsp.energy(sweep)
# plot auto-correlation
with pf.plot.context():
plt.plot(np.arange(-2**7, 2**7), auto_correlation)
plt.gca().set_xlim(-2**7, 2**7)
plt.gca().set_xlabel('time lag in samples')
plt.gca().set_ylabel('auto correlation')I'm this leaving this as it was...
There was a problem hiding this comment.
Mh, this is a peculiar issue with the integral approximation in SciPy I'd assume...
pyfar/signals/deterministic.py
Outdated
| n_samples, frequency_range, start_margin, stop_margin, fade_in=None, | ||
| fade_out=None, butterworth_order=8, sampling_rate=44100): |
There was a problem hiding this comment.
Why is the default for the fade in and out None, when the docstring states that the default is 0 samples?
There was a problem hiding this comment.
You might have reviewed an outdated version. The docstring mentions that the default is None after someone found this in an earlier round.
pyfar/signals/deterministic.py
Outdated
| butterworth_order : int, None | ||
| The order of the Butterworth filters that are applied to limit the | ||
| frequency range by a high-pass if ``frequency_range[0] > 0`` and/or by | ||
| a low-pass if ``frequency_range[1] < sampling_rate / 2``. The default | ||
| is ``8``. |
There was a problem hiding this comment.
Maybe it's a bit more intuitive to call this argument bandpass order or similar?
There was a problem hiding this comment.
Good idea, renamed to bandpass_orderbut left the description as it was.
pyfar/signals/deterministic.py
Outdated
| S. Müller, P. Massarani. 'Transfer Function Measurement with Sweeps. | ||
| Directors Cut Including Previously Unreleased Material and some | ||
| Corrections. J. Audio Eng. Soc. 2001, 49 (6), 443-471 (all equations | ||
| referenced in the code refer to this). |
| >>> auto_correlation = np.empty(2**8) | ||
| >>> for idx, shift in enumerate(range(-2**7, 2**7)): | ||
| >>> auto_correlation[idx] = np.dot( | ||
| ... sweep.time.flatten(), | ||
| ... np.roll(sweep.time.flatten(), shift)) | ||
| >>> | ||
| >>> auto_correlation /= pf.dsp.energy(sweep) |
There was a problem hiding this comment.
Indeed, this can be done in a more compact way and especially without a loop. For example using https://docs.scipy.org/doc/scipy/reference/generated/scipy.signal.correlate.html
pyfar/signals/deterministic.py
Outdated
|
|
||
| The exponential sweep can also be generated in the time domain | ||
| (:py:func:`~exponential_sweep_time`). Frequency domain synthesis | ||
| exhibits smooth magnitude spectra in trade of a slightly irregular |
There was a problem hiding this comment.
| exhibits smooth magnitude spectra in trade of a slightly irregular | |
| exhibits smooth magnitude spectra at the cost of a slightly irregular |
There was a problem hiding this comment.
changed the wording everywhere
pyfar/signals/deterministic.py
Outdated
| The length of the sweep in samples. | ||
| sampling_rate : int, optional | ||
| The sampling rate in Hz. The default is ``44100``. | ||
| group_delay : boolean, optional |
There was a problem hiding this comment.
| group_delay : boolean, optional | |
| return_group_delay : boolean, optional |
pyfar/signals/deterministic.py
Outdated
| fade_in=None, | ||
| fade_out=None, | ||
| sampling_rate=sampling_rate, | ||
| group_delay=group_delay) |
There was a problem hiding this comment.
| group_delay=group_delay) | |
| return_group_delay=reuturn_group_delay) |
pyfar/signals/deterministic.py
Outdated
|
|
||
|
|
||
| def _frequency_domain_sweep( | ||
| n_samples, sweep_type, frequency_range, butterworth_order, |
pyfar/signals/deterministic.py
Outdated
| def _frequency_domain_sweep( | ||
| n_samples, sweep_type, frequency_range, butterworth_order, | ||
| start_margin, stop_margin, fade_in, fade_out, sampling_rate, | ||
| group_delay): |
There was a problem hiding this comment.
| group_delay): | |
| return_group_delay): |
pyfar/signals/deterministic.py
Outdated
|
|
||
| # normalize to time domain amplitude of almost one | ||
| # (to avoid clipping if written to fixed point wav file) | ||
| sweep = pyfar.dsp.normalize(sweep) * (1 - 2**-15) |
There was a problem hiding this comment.
| sweep = pyfar.dsp.normalize(sweep) * (1 - 2**-15) | |
| sweep = pyfar.dsp.normalize(sweep) |
That is not something this function should take care of. That's the problem of a wav writer!
pyfar/signals/deterministic.py
Outdated
| n_samples, frequency_range, start_margin, stop_margin, fade_in=None, | ||
| fade_out=None, butterworth_order=8, sampling_rate=44100): |
There was a problem hiding this comment.
You might have reviewed an outdated version. The docstring mentions that the default is None after someone found this in an earlier round.
pyfar/signals/deterministic.py
Outdated
| butterworth_order : int, None | ||
| The order of the Butterworth filters that are applied to limit the | ||
| frequency range by a high-pass if ``frequency_range[0] > 0`` and/or by | ||
| a low-pass if ``frequency_range[1] < sampling_rate / 2``. The default | ||
| is ``8``. |
There was a problem hiding this comment.
Good idea, renamed to bandpass_orderbut left the description as it was.
pyfar/signals/deterministic.py
Outdated
| S. Müller, P. Massarani. 'Transfer Function Measurement with Sweeps. | ||
| Directors Cut Including Previously Unreleased Material and some | ||
| Corrections. J. Audio Eng. Soc. 2001, 49 (6), 443-471 (all equations | ||
| referenced in the code refer to this). |
pyfar/signals/deterministic.py
Outdated
| def exponential_sweep_freq( | ||
| n_samples, frequency_range, start_margin, stop_margin, fade_in=None, | ||
| fade_out=None, butterworth_order=8, sampling_rate=44100, | ||
| group_delay=False): |
There was a problem hiding this comment.
good point, changed everywhere
pyfar/signals/deterministic.py
Outdated
| group_delay : boolean, optional | ||
| Return the analytical group delay of the sweep. This can be used to | ||
| compute the times at which distortion products appear. The default is | ||
| ``False``. |
| >>> import pyfar as pf | ||
| >>> magnitude = pf.dsp.filter.low_shelve( | ||
| ... pf.signals.impulse(2**16), 500, 20, 2) | ||
| >>> magnitude = pf.dsp.filter.butterworth(magnitude, 8, 50, 'highpass') | ||
| >>> sweep = pf.signals.magnitude_spectrum_weighted_sweep( | ||
| ... 2**16, magnitude, 5000, 1000) | ||
| >>> pf.plot.time_freq(sweep) |
There was a problem hiding this comment.
Id like to keep it simple, here - especially while our interpolation function still overshoots :(
pyfar/signals/deterministic.py
Outdated
| system identification. It is orthogonal to delayed versions of itself and | ||
| can be played back in a loop. |
There was a problem hiding this comment.
removed it - although infinitely long signals are hard to play in a loop
| >>> auto_correlation = np.empty(2**8) | ||
| >>> for idx, shift in enumerate(range(-2**7, 2**7)): | ||
| >>> auto_correlation[idx] = np.dot( | ||
| ... sweep.time.flatten(), | ||
| ... np.roll(sweep.time.flatten(), shift)) | ||
| >>> | ||
| >>> auto_correlation /= pf.dsp.energy(sweep) |
There was a problem hiding this comment.
I tried this, but its not cyclic (regardless of mode and method) and causes ringing in the auto-correlation:
import pyfar as pf
import numpy as np
import scipy.signal as sgn
import matplotlib.pyplot as plt
%matplotlib inline
sweep = pf.signals.linear_perfect_sweep(2**8)
# compute auto-correlation
auto_correlation = sgn.correlate(
sweep.time.flatten(), sweep.time.flatten())
auto_correlation /= pf.dsp.energy(sweep)
# plot auto-correlation
with pf.plot.context():
plt.plot(np.arange(-2**7, 2**7), auto_correlation)
plt.gca().set_xlim(-2**7, 2**7)
plt.gca().set_xlabel('time lag in samples')
plt.gca().set_ylabel('auto correlation')I'm this leaving this as it was...
pyfar/signals/deterministic.py
Outdated
|
|
||
| # normalize to time domain amplitude of almost one | ||
| # (to avoid clipping if written to fixed point wav file) | ||
| sweep = pyfar.dsp.normalize(sweep) * (1 - 2**-15) |
pyfar/signals/deterministic.py
Outdated
| n_samples, frequency_range, start_margin, stop_margin, fade_in=None, | ||
| fade_out=None, bandpass_order=8, sampling_rate=44100, |
There was a problem hiding this comment.
| n_samples, frequency_range, start_margin, stop_margin, fade_in=None, | |
| fade_out=None, bandpass_order=8, sampling_rate=44100, | |
| n_samples, frequency_range, start_margin, stop_margin, fade_in=0, | |
| fade_out=0, bandpass_order=8, sampling_rate=44100, |
Is this type conversion required here? The docstring states that the type is 'int'
Also, this parameter is inconsistent with the time domain sweep generation function, which uses n_fade_in.
There was a problem hiding this comment.
Thats true - renamed to n_fade_* and used 0 as the default
pyfar/signals/deterministic.py
Outdated
| n_samples, frequency_range, start_margin, stop_margin, fade_in=None, | ||
| fade_out=None, bandpass_order=8, sampling_rate=44100, |
There was a problem hiding this comment.
| n_samples, frequency_range, start_margin, stop_margin, fade_in=None, | |
| fade_out=None, bandpass_order=8, sampling_rate=44100, | |
| n_samples, frequency_range, start_margin, stop_margin, fade_in=0, | |
| fade_out=0, bandpass_order=8, sampling_rate=44100, |
see above
pyfar/signals/deterministic.py
Outdated
| n_samples, magnitude_spectrum, start_margin, stop_margin, | ||
| fade_in=None, fade_out=None, sampling_rate=44100, |
| >>> auto_correlation = np.empty(2**8) | ||
| >>> for idx, shift in enumerate(range(-2**7, 2**7)): | ||
| >>> auto_correlation[idx] = np.dot( | ||
| ... sweep.time.flatten(), | ||
| ... np.roll(sweep.time.flatten(), shift)) | ||
| >>> | ||
| >>> auto_correlation /= pf.dsp.energy(sweep) |
There was a problem hiding this comment.
Mh, this is a peculiar issue with the integral approximation in SciPy I'd assume...
| .. [#] S. Müller, P. Massarani. 'Transfer Function Measurement with Sweeps. | ||
| Directors Cut Including Previously Unreleased Material and some | ||
| Corrections.' J. Audio Eng. Soc. 2001, 49 (6), 443-471. (all | ||
| equations referenced in the code refer to this). |
There was a problem hiding this comment.
This citation is only suited for non-perfect sweeps, right?
There was a problem hiding this comment.
yes and no I'd say. It does not mention perfect sweeps if I remember correctly but the method generally also holds for the perfect linear sweep at least.
Changes proposed in this pull request:
Add and test frequency domain sweep synthesis
make group delay an optional return parameter
apply 'none' fft normalisation