From 7a35a41cccef7b551f8acc0f1e51b3620e90d814 Mon Sep 17 00:00:00 2001 From: Fabian Brinkmann Date: Wed, 2 Jun 2021 12:59:16 +0200 Subject: [PATCH 01/42] added docstring for discussion --- pyfar/signals/deterministic.py | 130 +++++++++++++++++++++++++++++++++ 1 file changed, 130 insertions(+) diff --git a/pyfar/signals/deterministic.py b/pyfar/signals/deterministic.py index cbb60df63..2d4e94d46 100644 --- a/pyfar/signals/deterministic.py +++ b/pyfar/signals/deterministic.py @@ -143,6 +143,14 @@ def linear_sweep(n_samples, frequency_range, n_fade_out=90, amplitude=1, seconds, and the frequency limits :math:`f_\\mathrm{low}` and :math:`f_\\mathrm{high}`. + .. note:: + The linear sweep can also be generated in the frequency domain + (see :py:func:`~general_sweep_synthesis`). 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 of a slightly irregular temporal + envelope. + Parameters ---------- n_samples : int @@ -200,6 +208,14 @@ def exponential_sweep(n_samples, frequency_range, n_fade_out=90, amplitude=1, seconds, and the frequency limits :math:`f_\\mathrm{low}` and :math:`f_\\mathrm{high}`. + .. note:: + The exponential sweep can also be generated in the frequency domain + (see :py:func:`~general_sweep_synthesis`). 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 of a slightly irregular temporal + envelope. + Parameters ---------- n_samples : int @@ -242,6 +258,120 @@ def exponential_sweep(n_samples, frequency_range, n_fade_out=90, amplitude=1, return signal +def general_sweep_synthesis( + n_samples, spectrum, start_margin=None, stop_margin=None, + requency_range=None, double=True, sampling_rate=44100): + """ + Frequency domain sweep synthesis with arbitrary magnitude response. + + TODO Link to fractional octave smoothing in Notes + + TODO Implement calculation of group delay of impulse responses + + TODO Examples + + Sweep sweep synthesis according to [#]_ + + .. note:: + The linear and exponential sweep can also be generated in the time + domain (see :py:func:`~linear_sweep`, :py:func:`~exponential_sweep`). + Frequency domain synthesis exhibits smooth magnitude spectra and in + trade of a slightly irregular temporal envelope. Time domain synthesis + exhibits a constant temporal envelope in trade of slight ripples in the + magnitude response. + + Parameters + ---------- + n_samples : int + The length of the sweep in samples. + spectrum : 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. + + start_margin : int, optional + The time in samples, at which the sweep starts. The start margin is + required because the frequency domain sweep synthesis has pre-ringing + in the time domain. Not required if `spectrum` is ``'perfect'``. + stop_margin : int, optional + Time in samples, at which the sweep stops. This is relative to + `n_samples`, e.g., a stop margin of 100 samples means that the sweep + ends at sample ``n_samples-10``. This is required, because the + frequency domain sweep synthesis has post-ringing in the time domain. + Not required if `spectrum` is ``'perfect'``. + frequency_range : array like, optional + Frequency range of the sweep given by the lower and upper cut-off + frequency in Hz. The restriction of the frequency range is realized + by appling 8th order Butterworth filters at the specified frequencies. + Not required if `spectrum` is ``'perfect'`` or `signal`. + double : bool, optional + Double `n_samples` during the sweep calculation (recommended). The + default is ``True``. Not required if `spectrum` is ``'perfect'``. + sampling_rate : int, optional + The sampling rate in Hz. The default is ``44100``. + + Returns + ------- + sweep : Signal + The sweep signal. The Signal is in the time domain and has the ``none`` + FFT normalization (see :py:func:`~pyfar.dsp.fft.normalization`). The + sweep parameters are written to `comment`. + group_delay_sweep : FrequencyData + The group delay of the sweep as a single sided spectrum between 0 Hz + and ``sampling_rate/2``. + + TODO add this after implementation is complete: + + This can be used to calculate the group delay of the impulse responses + of linear and harmonic distortion products after deconvoloution (see + :py:func:`~pyfar.dsp...`). + + Notes + ----- + The envelope of the sweep time signal should be constant, appart from + slight overshoots at the beginning and end. If this is not the case, try to + provide a smoother spectrum (if `spectrum` is `signal`) or increase + `n_samples`. + + References + ---------- + .. [#] 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. + .. [#] C. Antweiler, A. Telle, P. Vary, G. Enzner. 'Perfect-Sweep NLMS for + Time-Variant Acoustic System Identification,' IEEE Int. Conf. + Acoustics, Speech and Signal Processing (ICASSP), + Prague, Czech Republic, 2011. doi: 10.1109/ICASSP.2012.6287930. + + Examples + -------- + TODO Example with spectrum=singal + (e.g., Bass emphasis by means of low shelve filter) + + TODO Examples with spectrum="linear" + + TODO Examples with spectrum="exponential" + + TODO Examples with spectrum="perfect" + """ + pass + + def _time_domain_sweep(n_samples, frequency_range, n_fade_out, amplitude, sampling_rate, sweep_type, sweep_rate=None): From 3eb99d32bdd028c2aadd11e6241b3db8031d1b63 Mon Sep 17 00:00:00 2001 From: Fabian Brinkmann Date: Wed, 2 Jun 2021 13:07:48 +0200 Subject: [PATCH 02/42] flake8 fix --- pyfar/signals/deterministic.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/pyfar/signals/deterministic.py b/pyfar/signals/deterministic.py index 2d4e94d46..c12fb1c38 100644 --- a/pyfar/signals/deterministic.py +++ b/pyfar/signals/deterministic.py @@ -259,8 +259,8 @@ def exponential_sweep(n_samples, frequency_range, n_fade_out=90, amplitude=1, def general_sweep_synthesis( - n_samples, spectrum, start_margin=None, stop_margin=None, - requency_range=None, double=True, sampling_rate=44100): + n_samples, spectrum, start_margin=None, stop_margin=None, + requency_range=None, double=True, sampling_rate=44100): """ Frequency domain sweep synthesis with arbitrary magnitude response. From b39f15ac2b0aacbe12ce59a58f9327642dc748ca Mon Sep 17 00:00:00 2001 From: Fabian Brinkmann Date: Thu, 17 Jun 2021 17:12:20 +0200 Subject: [PATCH 03/42] [skip ci] add function definitions for all new sweep synthesis methods --- pyfar/signals/deterministic.py | 56 +++++++++++++++++++++++++++++++--- 1 file changed, 51 insertions(+), 5 deletions(-) diff --git a/pyfar/signals/deterministic.py b/pyfar/signals/deterministic.py index 9ab4d9ad5..c55f0f19d 100644 --- a/pyfar/signals/deterministic.py +++ b/pyfar/signals/deterministic.py @@ -204,6 +204,17 @@ def linear_sweep_time(n_samples, frequency_range, n_fade_out=90, amplitude=1, return signal +def linear_sweep_freq( + n_samples, start_margin=None, stop_margin=None, frequency_range=None, + butterworth_order=8, double=True, sampling_rate=44100): + + signal, group_delay = _sweep_synthesis_freq( + n_samples, "linear", start_margin, stop_margin, + frequency_range, butterworth_order, double, sampling_rate) + + return signal, group_delay + + def exponential_sweep(n_samples, frequency_range, n_fade_out=90, amplitude=1, sweep_rate=None, sampling_rate=44100): """ @@ -287,9 +298,41 @@ def exponential_sweep_time(n_samples, frequency_range, n_fade_out=90, return signal -def general_sweep_synthesis( - n_samples, spectrum, start_margin=None, stop_margin=None, - requency_range=None, double=True, sampling_rate=44100): +def exponential_sweep_freq( + n_samples, start_margin=None, stop_margin=None, frequency_range=None, + butterworth_order=8, double=True, sampling_rate=44100): + + signal, group_delay = _sweep_synthesis_freq( + n_samples, "exponential", start_margin, stop_margin, + frequency_range, butterworth_order, double, sampling_rate) + + return signal, group_delay + + +def magnitude_weighted_sweep( + n_samples, magnitude, start_margin=None, stop_margin=None, + double=True, sampling_rate=44100): + + signal, group_delay = _sweep_synthesis_freq( + n_samples, magnitude, start_margin, stop_margin, + None, None, double, sampling_rate) + + return signal, group_delay + + +def perfect_sweep( + n_samples, sampling_rate=44100): + + signal, group_delay = _sweep_synthesis_freq( + n_samples, "perfect", None, None, None, None, False, sampling_rate) + + return signal, group_delay + + +def _sweep_synthesis_freq( + n_samples, magnitude, start_margin=None, stop_margin=None, + frequency_range=None, buterworth_order=8, double=True, + sampling_rate=44100): """ Frequency domain sweep synthesis with arbitrary magnitude response. @@ -313,7 +356,7 @@ def general_sweep_synthesis( ---------- n_samples : int The length of the sweep in samples. - spectrum : Signal, string + magnitude : Signal, string Specify the magnitude response of the sweep. signal @@ -346,8 +389,11 @@ def general_sweep_synthesis( frequency_range : array like, optional Frequency range of the sweep given by the lower and upper cut-off frequency in Hz. The restriction of the frequency range is realized - by appling 8th order Butterworth filters at the specified frequencies. + by appling a Butterworth band-pass with the specified frequencies. Not required if `spectrum` is ``'perfect'`` or `signal`. + butterworth_order : int, optional + The order of the Butterworth filters that are applied to limit the + frequency range. The default is ``8``. double : bool, optional Double `n_samples` during the sweep calculation (recommended). The default is ``True``. Not required if `spectrum` is ``'perfect'``. From 7a4647eeab4a94068aa8dca2fd6bf561f27db8b8 Mon Sep 17 00:00:00 2001 From: Fabian Brinkmann Date: Fri, 18 Jun 2021 12:43:02 +0200 Subject: [PATCH 04/42] [skip ci] suggestion from mberz --- pyfar/signals/deterministic.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pyfar/signals/deterministic.py b/pyfar/signals/deterministic.py index c55f0f19d..0552cd141 100644 --- a/pyfar/signals/deterministic.py +++ b/pyfar/signals/deterministic.py @@ -309,7 +309,7 @@ def exponential_sweep_freq( return signal, group_delay -def magnitude_weighted_sweep( +def magnitude_spectrum_weighted_sweep( n_samples, magnitude, start_margin=None, stop_margin=None, double=True, sampling_rate=44100): From b048bb8927e2be8fa1108539b7f1c5276d890d06 Mon Sep 17 00:00:00 2001 From: Fabian Brinkmann Date: Fri, 18 Jun 2021 13:36:36 +0200 Subject: [PATCH 05/42] unique naming Co-authored-by: Marco Berzborn --- pyfar/signals/deterministic.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pyfar/signals/deterministic.py b/pyfar/signals/deterministic.py index 862b80611..28300b2d9 100644 --- a/pyfar/signals/deterministic.py +++ b/pyfar/signals/deterministic.py @@ -322,7 +322,7 @@ def magnitude_spectrum_weighted_sweep( return signal, group_delay -def perfect_sweep( +def perfect_linear_sweep( n_samples, sampling_rate=44100): signal, group_delay = _sweep_synthesis_freq( From 89573c6420ff633aa58eff59daa9256a80acad5e Mon Sep 17 00:00:00 2001 From: Fabian Brinkmann Date: Mon, 28 Jun 2021 10:57:16 +0200 Subject: [PATCH 06/42] [skip ci] changed docstring --- pyfar/signals/deterministic.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pyfar/signals/deterministic.py b/pyfar/signals/deterministic.py index 28300b2d9..c7b35ede2 100644 --- a/pyfar/signals/deterministic.py +++ b/pyfar/signals/deterministic.py @@ -422,7 +422,7 @@ def _sweep_synthesis_freq( ----- The envelope of the sweep time signal should be constant, appart from slight overshoots at the beginning and end. If this is not the case, try to - provide a smoother spectrum (if `spectrum` is `signal`) or increase + provide a smoother spectrum (if `magnitude` is `signal`) or increase `n_samples`. References From e7199462d321d148e2da8899d8c8d0d5a0022e45 Mon Sep 17 00:00:00 2001 From: Alexander Deutsch Date: Fri, 21 Jan 2022 14:07:58 +0100 Subject: [PATCH 07/42] [skip ci] started to port implementation from AKtools, not ready --- pyfar/signals/deterministic.py | 113 ++++++++++++++++++++++++++++++++- 1 file changed, 112 insertions(+), 1 deletion(-) diff --git a/pyfar/signals/deterministic.py b/pyfar/signals/deterministic.py index c7b35ede2..719556b2a 100644 --- a/pyfar/signals/deterministic.py +++ b/pyfar/signals/deterministic.py @@ -446,7 +446,118 @@ def _sweep_synthesis_freq( TODO Examples with spectrum="perfect" """ - pass + + # calculte min and max frequency range limits for default frequency_range + delta_f = sampling_rate / n_samples + nyq = int(n_samples / 2) + 1 + t_end = n_samples / sampling_rate + if frequency_range is None: + f_max = sampling_rate / 2 - delta_f + frequency_range = [delta_f, f_max] + + # check input + if not isinstance(magnitude, (pyfar.Signal, str)): + raise TypeError("magnitude must be type Signal or str.") + if isinstance(magnitude, str) and magnitude not in ['linear', + 'exponential', + 'perfect']: + raise ValueError("magnitude must be 'linear', 'exponential' or", + "'perfect', when its type is str.") + if np.atleast_1d(frequency_range).size != 2: + raise ValueError( + "frequency_range must be an array like with to elements.") + if frequency_range[1] > sampling_rate/2: + raise ValueError( + "Upper frequency limit is larger than half the sampling rate.") + if frequency_range[0] == 0 and magnitude == "exponential": + raise ValueError("The exponential sweep can not start at 0 Hz.") + + # define start_margin and stop_margin if not given + if start_margin is None: + start_margin = 0 + if stop_margin is None: + stop_margin = 0 + + # double n_samples + if double and magnitude != 'perfect': + n_samples *= 2 + + # zero pad magnitude Signal, if needed + if isinstance(magnitude, pyfar.Signal): + if n_samples > magnitude.n_samples: + magnitude = pyfar.dsp.pad_zeros(magnitude, + n_samples - magnitude.n_samples) + + if magnitude == 'linear': + h_sweep = np.ones(int(n_samples)) + elif magnitude == 'exponential': + h_sweep = np.zeros(int(n_samples)) + h_sweep[1:int(n_samples/2)] = 1 / np.sqrt(2 * np.pi * np.linspace(frequency_range[0], frequency_range[1], len(h_sweep[1:int(n_samples/2)]))) + elif magnitude == 'perfect': + pass + + # mirror halfvector up to fs-df and get absolute values + sweep_abs = np.abs(h_sweep) + sweep_abs = np.concatenate((sweep_abs, + np.zeros(1), + sweep_abs[1:][::-1]), axis=None) + + # initialize group delay + sweep_tg = np.zeros(nyq + 1) + + # groupdelay at nyquist frequency + if double: + sweep_tg[nyq] = t_end / 2 - stop_margin / sampling_rate + else: + sweep_tg[nyq] = t_end - stop_margin / sampling_rate + + # groupdelay at DC + sweep_tg[0] = 0 + # groupdelay for first frequency bin + sweep_tg[1] = start_margin / sampling_rate + + # FORMULA (11, p.40 ) + sweep_power = np.sum(np.abs(sweep_abs[2:nyq] ** 2)) + C = (sweep_tg[0] - sweep_tg[1]) / sweep_power + + # FORMULA (10, p.40 ) + for k in range(2, nyq): + sweep_tg[k] = sweep_tg[k+1] + C * np.abs(sweep_abs[k]) ** 2 + + # calculate phase from group delay + sweep_phase = -1 * np.cumsum(sweep_tg) * 2 * np.pi * delta_f + + # mirror phase up to sampling_rate - delta_f + sweep_phase = np.concatenate((sweep_phase[:nyq], + (-1 * sweep_phase[1:nyq][::-1])), axis=None) + sweep_phase[0] = np.pi + + # calculate the complex spectrum + sweep_freq = sweep_abs * np.exp(1j * sweep_phase) + sweep_freq[0] = np.abs(sweep_freq[0]) + + sweep = pyfar.Signal(sweep_freq, sampling_rate, domain='freq') + + # transform into time domain + sweep.domain = 'time' + + # put group delay into FrequencyData + group_delay_sweep = pyfar.FrequencyData(sweep_tg[:nyq], + np.arange(0, sampling_rate/2, + step=delta_f)) + + # normalize to avoid clipping when written to wav with 16bit + # (worst case; LSB = 2^-15) + sweep = sweep / np.max(np.abs(sweep.time)) * (1 - 2**(-15)) + + # cut the second half of sweep, if n_samples was doubled while synthesis + if double: + sweep = pyfar.Signal(sweep.time[:int(sweep.n_samples/2)], + sampling_rate) + + # add comment + + return sweep, group_delay_sweep def _time_domain_sweep(n_samples, frequency_range, n_fade_out, amplitude, From 299acc62a4aa27f3c09e4b3680a0495bb212d40c Mon Sep 17 00:00:00 2001 From: Alexander Deutsch Date: Fri, 10 Jun 2022 12:04:10 +0200 Subject: [PATCH 08/42] [skip ci] add current status of _sweep_synthesis_freq, not working correctly --- pyfar/signals/deterministic.py | 107 +++++++++++++++------------------ 1 file changed, 48 insertions(+), 59 deletions(-) diff --git a/pyfar/signals/deterministic.py b/pyfar/signals/deterministic.py index 719556b2a..a69be8536 100644 --- a/pyfar/signals/deterministic.py +++ b/pyfar/signals/deterministic.py @@ -450,22 +450,21 @@ def _sweep_synthesis_freq( # calculte min and max frequency range limits for default frequency_range delta_f = sampling_rate / n_samples nyq = int(n_samples / 2) + 1 - t_end = n_samples / sampling_rate if frequency_range is None: f_max = sampling_rate / 2 - delta_f frequency_range = [delta_f, f_max] # check input if not isinstance(magnitude, (pyfar.Signal, str)): - raise TypeError("magnitude must be type Signal or str.") + raise TypeError("Magnitude must be type Signal or str.") if isinstance(magnitude, str) and magnitude not in ['linear', 'exponential', 'perfect']: - raise ValueError("magnitude must be 'linear', 'exponential' or", + raise ValueError("Magnitude must be 'linear', 'exponential' or", "'perfect', when its type is str.") if np.atleast_1d(frequency_range).size != 2: raise ValueError( - "frequency_range must be an array like with to elements.") + "Frequency_range must be an array like with two elements.") if frequency_range[1] > sampling_rate/2: raise ValueError( "Upper frequency limit is larger than half the sampling rate.") @@ -488,76 +487,66 @@ def _sweep_synthesis_freq( magnitude = pyfar.dsp.pad_zeros(magnitude, n_samples - magnitude.n_samples) - if magnitude == 'linear': - h_sweep = np.ones(int(n_samples)) - elif magnitude == 'exponential': - h_sweep = np.zeros(int(n_samples)) - h_sweep[1:int(n_samples/2)] = 1 / np.sqrt(2 * np.pi * np.linspace(frequency_range[0], frequency_range[1], len(h_sweep[1:int(n_samples/2)]))) - elif magnitude == 'perfect': - pass - - # mirror halfvector up to fs-df and get absolute values - sweep_abs = np.abs(h_sweep) - sweep_abs = np.concatenate((sweep_abs, - np.zeros(1), - sweep_abs[1:][::-1]), axis=None) + # spacing between frequency bins of FFT + df = sampling_rate / n_samples + # index of nyquist frequency + nyq = int(n_samples / 2) - # initialize group delay - sweep_tg = np.zeros(nyq + 1) + # magnitude coloration + if magnitude == 'lin': + h_sweep = np.ones(int(n_samples / 2)) + elif magnitude == 'exp': + h_sweep = np.zeros(int(n_samples / 2)) + f_min = df + f_max = sampling_rate / 2 + h_sweep[1:] = 1 / np.sqrt(2 * np.pi * np.arange(f_min, f_max, df)) - # groupdelay at nyquist frequency - if double: - sweep_tg[nyq] = t_end / 2 - stop_margin / sampling_rate - else: - sweep_tg[nyq] = t_end - stop_margin / sampling_rate + # TODO: generate filter for frequency_range and multiply with h_sweep - # groupdelay at DC - sweep_tg[0] = 0 - # groupdelay for first frequency bin - sweep_tg[1] = start_margin / sampling_rate + # initialize group delay + tg = np.zeros(nyq) + tg[1] = start_margin + tg[-1] = stop_margin # FORMULA (11, p.40 ) - sweep_power = np.sum(np.abs(sweep_abs[2:nyq] ** 2)) - C = (sweep_tg[0] - sweep_tg[1]) / sweep_power + sweep_power = np.sum(np.abs(h_sweep**2)) + C = (tg[-1] - tg[1]) / sweep_power # FORMULA (10, p.40 ) - for k in range(2, nyq): - sweep_tg[k] = sweep_tg[k+1] + C * np.abs(sweep_abs[k]) ** 2 - - # calculate phase from group delay - sweep_phase = -1 * np.cumsum(sweep_tg) * 2 * np.pi * delta_f + for k in range(2, nyq): # index 2 to nyq + tg[k] = tg[k-1] + C * np.abs(h_sweep[k])**2 - # mirror phase up to sampling_rate - delta_f - sweep_phase = np.concatenate((sweep_phase[:nyq], - (-1 * sweep_phase[1:nyq][::-1])), axis=None) - sweep_phase[0] = np.pi + sweep_ang = -1 * np.cumsum(tg) * 2 * np.pi * df - # calculate the complex spectrum - sweep_freq = sweep_abs * np.exp(1j * sweep_phase) - sweep_freq[0] = np.abs(sweep_freq[0]) - - sweep = pyfar.Signal(sweep_freq, sampling_rate, domain='freq') + # TODO: add correct phase similar to AKtools AKsweepFD + """ + # correct phase + if sweep_ang[-1] != 0: - # transform into time domain - sweep.domain = 'time' + length = len(sweep_ang) + phi_end = sweep_ang[-1] + df_phase = sampling_rate / (length) - # put group delay into FrequencyData - group_delay_sweep = pyfar.FrequencyData(sweep_tg[:nyq], - np.arange(0, sampling_rate/2, - step=delta_f)) + temp = np.ones(length) + temp = np.cumsum(temp) - 1 + offset = df_phase * phi_end / (sampling_rate / 2) + correct = temp * offset - # normalize to avoid clipping when written to wav with 16bit - # (worst case; LSB = 2^-15) - sweep = sweep / np.max(np.abs(sweep.time)) * (1 - 2**(-15)) + sweep_ang = sweep_ang - correct + """ - # cut the second half of sweep, if n_samples was doubled while synthesis - if double: - sweep = pyfar.Signal(sweep.time[:int(sweep.n_samples/2)], - sampling_rate) + # combine magnitude and phase of sweep + SWEEP = h_sweep * np.exp(1j * sweep_ang) + SWEEP[0] = np.abs(SWEEP[0]) - # add comment + # put sweep in pyfar.Signal an transform to time domain + sweep = pyfar.Signal(SWEEP, + sampling_rate, + n_samples=n_samples, + domain='freq') + sweep.domain = 'time' - return sweep, group_delay_sweep + return sweep, pyfar.dsp.group_delay(sweep) def _time_domain_sweep(n_samples, frequency_range, n_fade_out, amplitude, From 9e49d03efdc0e1888363123c2178dae60cb46c80 Mon Sep 17 00:00:00 2001 From: Fabian Brinkmann Date: Fri, 17 Jun 2022 14:05:19 +0200 Subject: [PATCH 09/42] [skip ci] update frequency domain sweep --- pyfar/signals/deterministic.py | 37 +++++++++++++++++----------------- 1 file changed, 18 insertions(+), 19 deletions(-) diff --git a/pyfar/signals/deterministic.py b/pyfar/signals/deterministic.py index a69be8536..77ac3fb6b 100644 --- a/pyfar/signals/deterministic.py +++ b/pyfar/signals/deterministic.py @@ -447,12 +447,13 @@ def _sweep_synthesis_freq( TODO Examples with spectrum="perfect" """ - # calculte min and max frequency range limits for default frequency_range - delta_f = sampling_rate / n_samples - nyq = int(n_samples / 2) + 1 + # spacing between frequency bins of FFT + df = sampling_rate / n_samples + + # calculate min and max frequency range limits for default frequency_range if frequency_range is None: - f_max = sampling_rate / 2 - delta_f - frequency_range = [delta_f, f_max] + f_max = sampling_rate / 2 - df + frequency_range = [df, f_max] # check input if not isinstance(magnitude, (pyfar.Signal, str)): @@ -481,22 +482,20 @@ def _sweep_synthesis_freq( if double and magnitude != 'perfect': n_samples *= 2 + # get number of bins (works for even and odd n_samples) + n_bins = n_samples // 2 + 1 + # zero pad magnitude Signal, if needed if isinstance(magnitude, pyfar.Signal): if n_samples > magnitude.n_samples: magnitude = pyfar.dsp.pad_zeros(magnitude, n_samples - magnitude.n_samples) - # spacing between frequency bins of FFT - df = sampling_rate / n_samples - # index of nyquist frequency - nyq = int(n_samples / 2) - # magnitude coloration - if magnitude == 'lin': - h_sweep = np.ones(int(n_samples / 2)) - elif magnitude == 'exp': - h_sweep = np.zeros(int(n_samples / 2)) + if magnitude == 'linear': + h_sweep = np.ones(n_bins) + elif magnitude == 'exponential': + h_sweep = np.zeros(n_bins) f_min = df f_max = sampling_rate / 2 h_sweep[1:] = 1 / np.sqrt(2 * np.pi * np.arange(f_min, f_max, df)) @@ -504,16 +503,16 @@ def _sweep_synthesis_freq( # TODO: generate filter for frequency_range and multiply with h_sweep # initialize group delay - tg = np.zeros(nyq) - tg[1] = start_margin - tg[-1] = stop_margin + tg = np.zeros(n_bins) + tg[1] = start_margin / sampling_rate + tg[-1] = (n_samples - stop_margin) / sampling_rate # FORMULA (11, p.40 ) sweep_power = np.sum(np.abs(h_sweep**2)) C = (tg[-1] - tg[1]) / sweep_power # FORMULA (10, p.40 ) - for k in range(2, nyq): # index 2 to nyq + for k in range(2, n_bins): # index 2 to nyq tg[k] = tg[k-1] + C * np.abs(h_sweep[k])**2 sweep_ang = -1 * np.cumsum(tg) * 2 * np.pi * df @@ -546,7 +545,7 @@ def _sweep_synthesis_freq( domain='freq') sweep.domain = 'time' - return sweep, pyfar.dsp.group_delay(sweep) + return sweep, tg def _time_domain_sweep(n_samples, frequency_range, n_fade_out, amplitude, From e51f2b2272074e28fb662eb71fbce49c0f916432 Mon Sep 17 00:00:00 2001 From: Fabian Brinkmann Date: Wed, 18 Oct 2023 15:09:37 +0200 Subject: [PATCH 10/42] remove defaults from private function to be more explicit --- pyfar/signals/deterministic.py | 57 ++++++++++++++-------------------- 1 file changed, 24 insertions(+), 33 deletions(-) diff --git a/pyfar/signals/deterministic.py b/pyfar/signals/deterministic.py index f4e45fc67..f74081fea 100644 --- a/pyfar/signals/deterministic.py +++ b/pyfar/signals/deterministic.py @@ -308,15 +308,14 @@ def perfect_linear_sweep( n_samples, sampling_rate=44100): signal, group_delay = _sweep_synthesis_freq( - n_samples, "perfect", None, None, None, None, False, sampling_rate) + n_samples, "perfect_linear", None, None, None, None, False, sampling_rate) return signal, group_delay def _sweep_synthesis_freq( - n_samples, magnitude, start_margin=None, stop_margin=None, - frequency_range=None, buterworth_order=8, double=True, - sampling_rate=44100): + n_samples, magnitude, start_margin, stop_margin, frequency_range, + buterworth_order, double, sampling_rate): """ Frequency domain sweep synthesis with arbitrary magnitude response. @@ -355,34 +354,34 @@ def _sweep_synthesis_freq( 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 + ``'perfect_linear'`` + Perfect linear sweep according to [#]_. Note that the parameters `start_margin`, `stop_margin`, `frequency_range` and `double` are not required in this case. start_margin : int, optional The time in samples, at which the sweep starts. The start margin is required because the frequency domain sweep synthesis has pre-ringing - in the time domain. Not required if `spectrum` is ``'perfect'``. + in the time domain. Set to ``0`` if `spectrum` is ``'perfect_linear'``. stop_margin : int, optional Time in samples, at which the sweep stops. This is relative to `n_samples`, e.g., a stop margin of 100 samples means that the sweep ends at sample ``n_samples-10``. This is required, because the frequency domain sweep synthesis has post-ringing in the time domain. - Not required if `spectrum` is ``'perfect'``. + Set to ``0`` if `spectrum` is ``'perfect_linear'``. frequency_range : array like, optional Frequency range of the sweep given by the lower and upper cut-off frequency in Hz. The restriction of the frequency range is realized by appling a Butterworth band-pass with the specified frequencies. - Not required if `spectrum` is ``'perfect'`` or `signal`. + Ignored if `spectrum` is ``'perfect_linear'`` or `signal`. butterworth_order : int, optional The order of the Butterworth filters that are applied to limit the - frequency range. The default is ``8``. + frequency range. double : bool, optional - Double `n_samples` during the sweep calculation (recommended). The - default is ``True``. Not required if `spectrum` is ``'perfect'``. + Double `n_samples` during the sweep calculation (recommended). Set to + ``False`` if `spectrum` is ``'perfect_linear'``. sampling_rate : int, optional - The sampling rate in Hz. The default is ``44100``. + The sampling rate in Hz. Returns ------- @@ -426,25 +425,17 @@ def _sweep_synthesis_freq( TODO Examples with spectrum="exponential" - TODO Examples with spectrum="perfect" + TODO Examples with spectrum="perfect_linear" """ - # spacing between frequency bins of FFT - df = sampling_rate / n_samples - - # calculate min and max frequency range limits for default frequency_range - if frequency_range is None: - f_max = sampling_rate / 2 - df - frequency_range = [df, f_max] - # check input if not isinstance(magnitude, (pyfar.Signal, str)): raise TypeError("Magnitude must be type Signal or str.") if isinstance(magnitude, str) and magnitude not in ['linear', 'exponential', - 'perfect']: + 'perfect_linear']: raise ValueError("Magnitude must be 'linear', 'exponential' or", - "'perfect', when its type is str.") + "'perfect_linear', when its type is str.") if np.atleast_1d(frequency_range).size != 2: raise ValueError( "Frequency_range must be an array like with two elements.") @@ -454,24 +445,24 @@ def _sweep_synthesis_freq( if frequency_range[0] == 0 and magnitude == "exponential": raise ValueError("The exponential sweep can not start at 0 Hz.") - # define start_margin and stop_margin if not given - if start_margin is None: - start_margin = 0 - if stop_margin is None: - stop_margin = 0 - # double n_samples - if double and magnitude != 'perfect': + if double and magnitude != 'perfect_linear': n_samples *= 2 + # spacing between frequency bins of FFT + df = sampling_rate / n_samples + # get number of bins (works for even and odd n_samples) n_bins = n_samples // 2 + 1 # zero pad magnitude Signal, if needed if isinstance(magnitude, pyfar.Signal): if n_samples > magnitude.n_samples: - magnitude = pyfar.dsp.pad_zeros(magnitude, - n_samples - magnitude.n_samples) + magnitude = pyfar.dsp.pad_zeros( + magnitude, n_samples - magnitude.n_samples) + elif magnitude.n_samples > n_samples: + raise ValueError((f'magnitue can has {magnitude.n_samples} samples' + f' but must not be longer than {n_samples}')) # magnitude coloration if magnitude == 'linear': From 868564977fc942e732ca9a887c30cb0616773cb3 Mon Sep 17 00:00:00 2001 From: Fabian Brinkmann Date: Wed, 18 Oct 2023 15:35:35 +0200 Subject: [PATCH 11/42] progress --- pyfar/signals/deterministic.py | 55 ++++++++++++++++++++++------------ 1 file changed, 36 insertions(+), 19 deletions(-) diff --git a/pyfar/signals/deterministic.py b/pyfar/signals/deterministic.py index f74081fea..dd33d888c 100644 --- a/pyfar/signals/deterministic.py +++ b/pyfar/signals/deterministic.py @@ -308,14 +308,15 @@ def perfect_linear_sweep( n_samples, sampling_rate=44100): signal, group_delay = _sweep_synthesis_freq( - n_samples, "perfect_linear", None, None, None, None, False, sampling_rate) + n_samples, "perfect_linear", None, None, None, None, False, + sampling_rate) return signal, group_delay def _sweep_synthesis_freq( n_samples, magnitude, start_margin, stop_margin, frequency_range, - buterworth_order, double, sampling_rate): + butterworth_order, double, sampling_rate): """ Frequency domain sweep synthesis with arbitrary magnitude response. @@ -362,24 +363,25 @@ def _sweep_synthesis_freq( start_margin : int, optional The time in samples, at which the sweep starts. The start margin is required because the frequency domain sweep synthesis has pre-ringing - in the time domain. Set to ``0`` if `spectrum` is ``'perfect_linear'``. + in the time domain. Set to ``0`` if `magnitude` is ``'perfect_linear'``. stop_margin : int, optional Time in samples, at which the sweep stops. This is relative to `n_samples`, e.g., a stop margin of 100 samples means that the sweep ends at sample ``n_samples-10``. This is required, because the frequency domain sweep synthesis has post-ringing in the time domain. - Set to ``0`` if `spectrum` is ``'perfect_linear'``. + Set to ``0`` if `magnitude` is ``'perfect_linear'``. frequency_range : array like, optional Frequency range of the sweep given by the lower and upper cut-off frequency in Hz. The restriction of the frequency range is realized by appling a Butterworth band-pass with the specified frequencies. - Ignored if `spectrum` is ``'perfect_linear'`` or `signal`. + Ignored if `magnitude` is ``'perfect_linear'`` or `signal`. butterworth_order : int, optional The order of the Butterworth filters that are applied to limit the - frequency range. + frequency range. Must be ``None`` if `magnitude` is + ``'perfect_linear'``. double : bool, optional Double `n_samples` during the sweep calculation (recommended). Set to - ``False`` if `spectrum` is ``'perfect_linear'``. + ``False`` if `magnitude` is ``'perfect_linear'``. sampling_rate : int, optional The sampling rate in Hz. @@ -418,14 +420,18 @@ def _sweep_synthesis_freq( Examples -------- - TODO Example with spectrum=singal + TODO Example with magnitude=singal (e.g., Bass emphasis by means of low shelve filter) - TODO Examples with spectrum="linear" + TODO Examples with magnitude="linear" - TODO Examples with spectrum="exponential" + TODO Examples with magnitude="exponential" - TODO Examples with spectrum="perfect_linear" + TODO Examples with magnitude="perfect_linear" + + TODO rename variable h_sweep to sweep_abs + + TODO rename SWEEP to sweep """ # check input @@ -455,27 +461,38 @@ def _sweep_synthesis_freq( # get number of bins (works for even and odd n_samples) n_bins = n_samples // 2 + 1 - # zero pad magnitude Signal, if needed + # handle magnitude parameter if isinstance(magnitude, pyfar.Signal): + # zero pad magnitude Signal or raise error if needed if n_samples > magnitude.n_samples: magnitude = pyfar.dsp.pad_zeros( magnitude, n_samples - magnitude.n_samples) elif magnitude.n_samples > n_samples: raise ValueError((f'magnitue can has {magnitude.n_samples} samples' f' but must not be longer than {n_samples}')) - - # magnitude coloration - if magnitude == 'linear': + h_sweep = np.abs(magnitude.freq_raw) + elif magnitude == 'linear': h_sweep = np.ones(n_bins) elif magnitude == 'exponential': h_sweep = np.zeros(n_bins) f_min = df f_max = sampling_rate / 2 h_sweep[1:] = 1 / np.sqrt(2 * np.pi * np.arange(f_min, f_max, df)) - - # TODO: generate filter for frequency_range and multiply with h_sweep - - # initialize group delay + elif magnitude == 'perfect_linear': + if start_margin != 0 or stop_margin != 0 or double or \ + butterworth_order is not None or \ + np.array(frequency_range) != np.array([0, sampling_rate / 2]): + # internal warning. Users will not call this function directly + # and can not cause this error. + raise ValueError(('Found conflicting parameters')) + + # apply band limit to magnitude + if magnitude in ['linear', 'exponential']: + bandpass = pyfar.dsp.filter.butterworth( + pyfar.signals.impulse(n_samples, sampling_rate=sampling_rate)) + h_sweep *= np.abs(bandpass.freq.flatten()) + + # initialize group delay in seconds at 0 Hz, df and Nyquist tg = np.zeros(n_bins) tg[1] = start_margin / sampling_rate tg[-1] = (n_samples - stop_margin) / sampling_rate From 382cf149dd70bda5807281093a350d477d4ed1b7 Mon Sep 17 00:00:00 2001 From: Fabian Brinkmann Date: Wed, 18 Oct 2023 19:00:09 +0200 Subject: [PATCH 12/42] first working version --- pyfar/signals/deterministic.py | 106 +++++++++++++++++++-------------- 1 file changed, 61 insertions(+), 45 deletions(-) diff --git a/pyfar/signals/deterministic.py b/pyfar/signals/deterministic.py index dd33d888c..797d1e4d4 100644 --- a/pyfar/signals/deterministic.py +++ b/pyfar/signals/deterministic.py @@ -360,29 +360,29 @@ def _sweep_synthesis_freq( `start_margin`, `stop_margin`, `frequency_range` and `double` are not required in this case. - start_margin : int, optional + start_margin : int, float The time in samples, at which the sweep starts. The start margin is required because the frequency domain sweep synthesis has pre-ringing in the time domain. Set to ``0`` if `magnitude` is ``'perfect_linear'``. - stop_margin : int, optional + stop_margin : int, float Time in samples, at which the sweep stops. This is relative to `n_samples`, e.g., a stop margin of 100 samples means that the sweep ends at sample ``n_samples-10``. This is required, because the frequency domain sweep synthesis has post-ringing in the time domain. Set to ``0`` if `magnitude` is ``'perfect_linear'``. - frequency_range : array like, optional + frequency_range : array like Frequency range of the sweep given by the lower and upper cut-off frequency in Hz. The restriction of the frequency range is realized by appling a Butterworth band-pass with the specified frequencies. Ignored if `magnitude` is ``'perfect_linear'`` or `signal`. - butterworth_order : int, optional + butterworth_order : int, None The order of the Butterworth filters that are applied to limit the - frequency range. Must be ``None`` if `magnitude` is - ``'perfect_linear'``. - double : bool, optional + frequency range by a high-pass if ``frequency_range[0] > 0`` and/or by + a low-pass if ``frequency_range[1] < sampling_rate / 2``. + double : bool Double `n_samples` during the sweep calculation (recommended). Set to ``False`` if `magnitude` is ``'perfect_linear'``. - sampling_rate : int, optional + sampling_rate : int The sampling rate in Hz. Returns @@ -431,7 +431,7 @@ def _sweep_synthesis_freq( TODO rename variable h_sweep to sweep_abs - TODO rename SWEEP to sweep + TODO make group delay a FrequencyData object """ # check input @@ -449,10 +449,20 @@ def _sweep_synthesis_freq( raise ValueError( "Upper frequency limit is larger than half the sampling rate.") if frequency_range[0] == 0 and magnitude == "exponential": - raise ValueError("The exponential sweep can not start at 0 Hz.") + Warning(( + "The exponential sweep has a 1/frequency magnitude spectrum. " + "The magnitude is set to 0 at 0 Hz to avoid division by zero.")) + if magnitude == 'perfect_linear' and \ + (start_margin != 0 or stop_margin != 0 or double or + frequency_range[0] != 0 or + frequency_range[1] != sampling_rate / 2): + # internal warning. Users will not call this function directly + # and can not cause this error. + raise ValueError(('Found conflicting parameters')) # double n_samples if double and magnitude != 'perfect_linear': + stop_margin += n_samples n_samples *= 2 # spacing between frequency bins of FFT @@ -471,26 +481,31 @@ def _sweep_synthesis_freq( raise ValueError((f'magnitue can has {magnitude.n_samples} samples' f' but must not be longer than {n_samples}')) h_sweep = np.abs(magnitude.freq_raw) - elif magnitude == 'linear': + elif magnitude in ['linear', 'perfect_linear']: + # constant spectrum h_sweep = np.ones(n_bins) elif magnitude == 'exponential': + # 1/f spectrum h_sweep = np.zeros(n_bins) - f_min = df - f_max = sampling_rate / 2 - h_sweep[1:] = 1 / np.sqrt(2 * np.pi * np.arange(f_min, f_max, df)) - elif magnitude == 'perfect_linear': - if start_margin != 0 or stop_margin != 0 or double or \ - butterworth_order is not None or \ - np.array(frequency_range) != np.array([0, sampling_rate / 2]): - # internal warning. Users will not call this function directly - # and can not cause this error. - raise ValueError(('Found conflicting parameters')) + h_sweep[1:] = 1 / np.sqrt(2 * np.pi * np.arange(1, n_bins) * df) # apply band limit to magnitude if magnitude in ['linear', 'exponential']: - bandpass = pyfar.dsp.filter.butterworth( - pyfar.signals.impulse(n_samples, sampling_rate=sampling_rate)) - h_sweep *= np.abs(bandpass.freq.flatten()) + if frequency_range[0] > 0 and frequency_range[1] < sampling_rate / 2: + band_limit = pyfar.dsp.filter.butterworth( + pyfar.signals.impulse(n_samples, sampling_rate=sampling_rate), + butterworth_order, frequency_range, 'bandpass') + elif frequency_range[0] > 0: + band_limit = pyfar.dsp.filter.butterworth( + pyfar.signals.impulse(n_samples, sampling_rate=sampling_rate), + butterworth_order, frequency_range[0], 'highpass') + elif frequency_range[1] < sampling_rate / 2: + band_limit = pyfar.dsp.filter.butterworth( + pyfar.signals.impulse(n_samples, sampling_rate=sampling_rate), + butterworth_order, frequency_range[1], 'lowpass') + else: + band_limit = np.ones_like(h_sweep) + h_sweep *= np.abs(band_limit.freq.flatten()) # initialize group delay in seconds at 0 Hz, df and Nyquist tg = np.zeros(n_bins) @@ -507,33 +522,34 @@ def _sweep_synthesis_freq( sweep_ang = -1 * np.cumsum(tg) * 2 * np.pi * df - # TODO: add correct phase similar to AKtools AKsweepFD - """ - # correct phase - if sweep_ang[-1] != 0: + # wrap and correct phase to be real 0 at Nyquist + sweep_ang = pyfar.dsp.wrap_to_2pi(sweep_ang) + sweep_ang[sweep_ang > np.pi] -= 2*np.pi - length = len(sweep_ang) - phi_end = sweep_ang[-1] - df_phase = sampling_rate / (length) + if sweep_ang[-1] != 0 and not n_samples % 2: - temp = np.ones(length) - temp = np.cumsum(temp) - 1 - offset = df_phase * phi_end / (sampling_rate / 2) - correct = temp * offset - - sweep_ang = sweep_ang - correct - """ + factor = np.cumsum(np.ones_like(sweep_ang)) - 1 + offset = df * sweep_ang[-1] / (sampling_rate / 2) + sweep_ang -= factor * offset + sweep_ang[-1] = np.abs(sweep_ang[-1]) # combine magnitude and phase of sweep - SWEEP = h_sweep * np.exp(1j * sweep_ang) - SWEEP[0] = np.abs(SWEEP[0]) + sweep = h_sweep * np.exp(1j * sweep_ang) # put sweep in pyfar.Signal an transform to time domain - sweep = pyfar.Signal(SWEEP, - sampling_rate, - n_samples=n_samples, - domain='freq') - sweep.domain = 'time' + sweep = pyfar.Signal(sweep, sampling_rate, n_samples, 'freq', 'rms') + + # cut to originally desired length + if double: + n_samples = n_samples // 2 + stop_margin -= n_samples + sweep.time = sweep.time[..., :n_samples] + + # TODO fade-in/out + + # 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) return sweep, tg From 3eec2c079a3d95b417d0c64ae45d37b438d6bebd Mon Sep 17 00:00:00 2001 From: Fabian Brinkmann Date: Thu, 19 Oct 2023 08:40:08 +0200 Subject: [PATCH 13/42] [skip ci] flake8 fix --- pyfar/signals/deterministic.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/pyfar/signals/deterministic.py b/pyfar/signals/deterministic.py index 797d1e4d4..3f5b4991b 100644 --- a/pyfar/signals/deterministic.py +++ b/pyfar/signals/deterministic.py @@ -363,7 +363,8 @@ def _sweep_synthesis_freq( start_margin : int, float The time in samples, at which the sweep starts. The start margin is required because the frequency domain sweep synthesis has pre-ringing - in the time domain. Set to ``0`` if `magnitude` is ``'perfect_linear'``. + in the time domain. Set to ``0`` if `magnitude` is + ``'perfect_linear'``. stop_margin : int, float Time in samples, at which the sweep stops. This is relative to `n_samples`, e.g., a stop margin of 100 samples means that the sweep From f2b9b6dd91615fe8cd54c65c3903cc6721dddc54 Mon Sep 17 00:00:00 2001 From: Fabian Brinkmann Date: Thu, 19 Oct 2023 09:55:55 +0200 Subject: [PATCH 14/42] further progress --- pyfar/signals/deterministic.py | 69 +++++++++++++++++++--------------- 1 file changed, 38 insertions(+), 31 deletions(-) diff --git a/pyfar/signals/deterministic.py b/pyfar/signals/deterministic.py index 3f5b4991b..7fac59919 100644 --- a/pyfar/signals/deterministic.py +++ b/pyfar/signals/deterministic.py @@ -363,8 +363,7 @@ def _sweep_synthesis_freq( start_margin : int, float The time in samples, at which the sweep starts. The start margin is required because the frequency domain sweep synthesis has pre-ringing - in the time domain. Set to ``0`` if `magnitude` is - ``'perfect_linear'``. + in the time domain. Set to ``0`` if `magnitude` is ``'perfect_linear'``. stop_margin : int, float Time in samples, at which the sweep stops. This is relative to `n_samples`, e.g., a stop margin of 100 samples means that the sweep @@ -429,13 +428,9 @@ def _sweep_synthesis_freq( TODO Examples with magnitude="exponential" TODO Examples with magnitude="perfect_linear" - - TODO rename variable h_sweep to sweep_abs - - TODO make group delay a FrequencyData object """ - # check input + # check input ------------------------------------------------------------- if not isinstance(magnitude, (pyfar.Signal, str)): raise TypeError("Magnitude must be type Signal or str.") if isinstance(magnitude, str) and magnitude not in ['linear', @@ -454,13 +449,14 @@ def _sweep_synthesis_freq( "The exponential sweep has a 1/frequency magnitude spectrum. " "The magnitude is set to 0 at 0 Hz to avoid division by zero.")) if magnitude == 'perfect_linear' and \ - (start_margin != 0 or stop_margin != 0 or double or - frequency_range[0] != 0 or + (start_margin != 0 or stop_margin != 0 or double or \ + frequency_range[0] != 0 or \ frequency_range[1] != sampling_rate / 2): # internal warning. Users will not call this function directly # and can not cause this error. raise ValueError(('Found conflicting parameters')) + # initialize basic parameters --------------------------------------------- # double n_samples if double and magnitude != 'perfect_linear': stop_margin += n_samples @@ -472,7 +468,7 @@ def _sweep_synthesis_freq( # get number of bins (works for even and odd n_samples) n_bins = n_samples // 2 + 1 - # handle magnitude parameter + # compute magnitude spectrum ---------------------------------------------- if isinstance(magnitude, pyfar.Signal): # zero pad magnitude Signal or raise error if needed if n_samples > magnitude.n_samples: @@ -481,16 +477,16 @@ def _sweep_synthesis_freq( elif magnitude.n_samples > n_samples: raise ValueError((f'magnitue can has {magnitude.n_samples} samples' f' but must not be longer than {n_samples}')) - h_sweep = np.abs(magnitude.freq_raw) + sweep_abs = np.abs(magnitude.freq_raw) elif magnitude in ['linear', 'perfect_linear']: # constant spectrum - h_sweep = np.ones(n_bins) + sweep_abs = np.ones(n_bins) elif magnitude == 'exponential': # 1/f spectrum - h_sweep = np.zeros(n_bins) - h_sweep[1:] = 1 / np.sqrt(2 * np.pi * np.arange(1, n_bins) * df) + sweep_abs = np.zeros(n_bins) + sweep_abs[1:] = 1 / np.sqrt(2 * np.pi * np.arange(1, n_bins) * df) - # apply band limit to magnitude + # band limit to magnitude spectrum if magnitude in ['linear', 'exponential']: if frequency_range[0] > 0 and frequency_range[1] < sampling_rate / 2: band_limit = pyfar.dsp.filter.butterworth( @@ -505,37 +501,48 @@ def _sweep_synthesis_freq( pyfar.signals.impulse(n_samples, sampling_rate=sampling_rate), butterworth_order, frequency_range[1], 'lowpass') else: - band_limit = np.ones_like(h_sweep) - h_sweep *= np.abs(band_limit.freq.flatten()) - - # initialize group delay in seconds at 0 Hz, df and Nyquist - tg = np.zeros(n_bins) - tg[1] = start_margin / sampling_rate - tg[-1] = (n_samples - stop_margin) / sampling_rate + band_limit = pyfar.Signal(np.ones_like(sweep_abs), sampling_rate, + n_samples, 'freq') + sweep_abs *= np.abs(band_limit.freq.flatten()) + + # compute group delay ----------------------------------------------------- + # group delay at 0 Hz must be 0 + sweep_gd = np.zeros(n_bins) + # group delay at df equals starting time unless it's 0 + if start_margin > 0: + sweep_gd[1] = start_margin / sampling_rate + tg_start = 2 + else: + tg_start = 1 + # group delay at Nyquist equals stopping time + sweep_gd[-1] = (n_samples - stop_margin) / sampling_rate # FORMULA (11, p.40 ) - sweep_power = np.sum(np.abs(h_sweep**2)) - C = (tg[-1] - tg[1]) / sweep_power + sweep_power = np.sum(np.abs(sweep_abs**2)) + C = (sweep_gd[-1] - sweep_gd[1]) / sweep_power # FORMULA (10, p.40 ) - for k in range(2, n_bins): # index 2 to nyq - tg[k] = tg[k-1] + C * np.abs(h_sweep[k])**2 + for k in range(tg_start, n_bins): # index 2 to nyq + sweep_gd[k] = sweep_gd[k-1] + C * np.abs(sweep_abs[k])**2 - sweep_ang = -1 * np.cumsum(tg) * 2 * np.pi * df + # compute phase from group delay ------------------------------------------ + sweep_ang = -1 * np.cumsum(sweep_gd) * 2 * np.pi * df # wrap and correct phase to be real 0 at Nyquist sweep_ang = pyfar.dsp.wrap_to_2pi(sweep_ang) sweep_ang[sweep_ang > np.pi] -= 2*np.pi - if sweep_ang[-1] != 0 and not n_samples % 2: - + if magnitude == 'perfect_linear': + sweep_ang[-1] = 0 + elif sweep_ang[-1] != 0 and not n_samples % 2: factor = np.cumsum(np.ones_like(sweep_ang)) - 1 offset = df * sweep_ang[-1] / (sampling_rate / 2) sweep_ang -= factor * offset sweep_ang[-1] = np.abs(sweep_ang[-1]) + # compute and finalize return data ---------------------------------------- # combine magnitude and phase of sweep - sweep = h_sweep * np.exp(1j * sweep_ang) + sweep = sweep_abs * np.exp(1j * sweep_ang) # put sweep in pyfar.Signal an transform to time domain sweep = pyfar.Signal(sweep, sampling_rate, n_samples, 'freq', 'rms') @@ -552,7 +559,7 @@ def _sweep_synthesis_freq( # (to avoid clipping if written to fixed point wav file) sweep = pyfar.dsp.normalize(sweep) * (1 - 2**-15) - return sweep, tg + return sweep, sweep_gd def _time_domain_sweep(n_samples, frequency_range, n_fade_out, amplitude, From 7d7a73e93e27671a982648da9c1be63b4c5c197a Mon Sep 17 00:00:00 2001 From: Fabian Brinkmann Date: Thu, 19 Oct 2023 10:02:53 +0200 Subject: [PATCH 15/42] return group delay as pyfar Audio object --- pyfar/signals/deterministic.py | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/pyfar/signals/deterministic.py b/pyfar/signals/deterministic.py index 7fac59919..c39ac9e38 100644 --- a/pyfar/signals/deterministic.py +++ b/pyfar/signals/deterministic.py @@ -547,6 +547,10 @@ def _sweep_synthesis_freq( # put sweep in pyfar.Signal an transform to time domain sweep = pyfar.Signal(sweep, sampling_rate, n_samples, 'freq', 'rms') + # put group delay on pyfar FrequencyData + sweep_gd = pyfar.FrequencyData( + sweep_gd, pyfar.dsp.fft.rfftfreq(n_samples, sampling_rate)) + # cut to originally desired length if double: n_samples = n_samples // 2 From be43e4e5f843d4ec84d11f36e3df6903ab138d2a Mon Sep 17 00:00:00 2001 From: Fabian Brinkmann Date: Thu, 19 Oct 2023 10:38:14 +0200 Subject: [PATCH 16/42] more flake fixes --- pyfar/signals/deterministic.py | 9 +++++---- 1 file changed, 5 insertions(+), 4 deletions(-) diff --git a/pyfar/signals/deterministic.py b/pyfar/signals/deterministic.py index c39ac9e38..abd35acf8 100644 --- a/pyfar/signals/deterministic.py +++ b/pyfar/signals/deterministic.py @@ -363,7 +363,8 @@ def _sweep_synthesis_freq( start_margin : int, float The time in samples, at which the sweep starts. The start margin is required because the frequency domain sweep synthesis has pre-ringing - in the time domain. Set to ``0`` if `magnitude` is ``'perfect_linear'``. + in the time domain. Set to ``0`` if `magnitude` is + ``'perfect_linear'``. stop_margin : int, float Time in samples, at which the sweep stops. This is relative to `n_samples`, e.g., a stop margin of 100 samples means that the sweep @@ -373,7 +374,7 @@ def _sweep_synthesis_freq( frequency_range : array like Frequency range of the sweep given by the lower and upper cut-off frequency in Hz. The restriction of the frequency range is realized - by appling a Butterworth band-pass with the specified frequencies. + by applying a Butterworth band-pass with the specified frequencies. Ignored if `magnitude` is ``'perfect_linear'`` or `signal`. butterworth_order : int, None The order of the Butterworth filters that are applied to limit the @@ -449,8 +450,8 @@ def _sweep_synthesis_freq( "The exponential sweep has a 1/frequency magnitude spectrum. " "The magnitude is set to 0 at 0 Hz to avoid division by zero.")) if magnitude == 'perfect_linear' and \ - (start_margin != 0 or stop_margin != 0 or double or \ - frequency_range[0] != 0 or \ + (start_margin != 0 or stop_margin != 0 or double or + frequency_range[0] != 0 or frequency_range[1] != sampling_rate / 2): # internal warning. Users will not call this function directly # and can not cause this error. From c6cd5586b7301e50c8e37a25f1028ebcab7c2f79 Mon Sep 17 00:00:00 2001 From: Fabian Brinkmann Date: Thu, 19 Oct 2023 10:55:26 +0200 Subject: [PATCH 17/42] add fade in and out --- pyfar/signals/deterministic.py | 27 +++++++++++++++++++++++---- 1 file changed, 23 insertions(+), 4 deletions(-) diff --git a/pyfar/signals/deterministic.py b/pyfar/signals/deterministic.py index abd35acf8..a1ff74dec 100644 --- a/pyfar/signals/deterministic.py +++ b/pyfar/signals/deterministic.py @@ -315,8 +315,8 @@ def perfect_linear_sweep( def _sweep_synthesis_freq( - n_samples, magnitude, start_margin, stop_margin, frequency_range, - butterworth_order, double, sampling_rate): + n_samples, magnitude, frequency_range, butterworth_order, double, + start_margin, stop_margin, fade_in, fade_out, sampling_rate): """ Frequency domain sweep synthesis with arbitrary magnitude response. @@ -558,13 +558,32 @@ def _sweep_synthesis_freq( stop_margin -= n_samples sweep.time = sweep.time[..., :n_samples] - # TODO fade-in/out + # apply window + if fade_in and fade_out: + fade = [start_margin, + start_margin + fade_in, + n_samples - stop_margin - fade_out, + n_samples - stop_margin] + shape = 'symmetric' + elif fade_in: + fade = [start_margin, + start_margin + fade_in] + shape = 'left' + elif fade_out: + fade = [n_samples - stop_margin - fade_out, + n_samples - stop_margin] + shape = 'right' + else: + fade = None + + if fade is not None: + sweep = pyfar.dsp.time_window(sweep, fade, shape=shape) # 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) - return sweep, sweep_gd + return sweep, sweep_gd, sweep_abs def _time_domain_sweep(n_samples, frequency_range, n_fade_out, amplitude, From 24bcd7452c6c1df8cc4a0b4592edbd2736b86b5c Mon Sep 17 00:00:00 2001 From: Fabian Brinkmann Date: Fri, 20 Oct 2023 11:08:06 +0200 Subject: [PATCH 18/42] renaming for consistency --- pyfar/signals/deterministic.py | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/pyfar/signals/deterministic.py b/pyfar/signals/deterministic.py index a1ff74dec..e72d7f788 100644 --- a/pyfar/signals/deterministic.py +++ b/pyfar/signals/deterministic.py @@ -206,7 +206,7 @@ def linear_sweep_freq( n_samples, start_margin=None, stop_margin=None, frequency_range=None, butterworth_order=8, double=True, sampling_rate=44100): - signal, group_delay = _sweep_synthesis_freq( + signal, group_delay = _frequency_domain_sweep( n_samples, "linear", start_margin, stop_margin, frequency_range, butterworth_order, double, sampling_rate) @@ -286,7 +286,7 @@ def exponential_sweep_freq( n_samples, start_margin=None, stop_margin=None, frequency_range=None, butterworth_order=8, double=True, sampling_rate=44100): - signal, group_delay = _sweep_synthesis_freq( + signal, group_delay = _frequency_domain_sweep( n_samples, "exponential", start_margin, stop_margin, frequency_range, butterworth_order, double, sampling_rate) @@ -297,7 +297,7 @@ def magnitude_spectrum_weighted_sweep( n_samples, magnitude, start_margin=None, stop_margin=None, double=True, sampling_rate=44100): - signal, group_delay = _sweep_synthesis_freq( + signal, group_delay = _frequency_domain_sweep( n_samples, magnitude, start_margin, stop_margin, None, None, double, sampling_rate) @@ -307,14 +307,14 @@ def magnitude_spectrum_weighted_sweep( def perfect_linear_sweep( n_samples, sampling_rate=44100): - signal, group_delay = _sweep_synthesis_freq( + signal, group_delay = _frequency_domain_sweep( n_samples, "perfect_linear", None, None, None, None, False, sampling_rate) return signal, group_delay -def _sweep_synthesis_freq( +def _frequency_domain_sweep( n_samples, magnitude, frequency_range, butterworth_order, double, start_margin, stop_margin, fade_in, fade_out, sampling_rate): """ From f23913737f4b630b90d5bb298a6e59e6bd978046 Mon Sep 17 00:00:00 2001 From: Fabian Brinkmann Date: Fri, 20 Oct 2023 11:55:41 +0200 Subject: [PATCH 19/42] first working version --- pyfar/signals/__init__.py | 6 +- pyfar/signals/deterministic.py | 124 +++++++++++++++++++++------------ 2 files changed, 84 insertions(+), 46 deletions(-) diff --git a/pyfar/signals/__init__.py b/pyfar/signals/__init__.py index 34ec00c61..460973d47 100644 --- a/pyfar/signals/__init__.py +++ b/pyfar/signals/__init__.py @@ -11,7 +11,9 @@ """ from .deterministic import ( - sine, impulse, linear_sweep_time, exponential_sweep_time) + sine, impulse, linear_sweep_time, exponential_sweep_time, + linear_sweep_freq, exponential_sweep_freq, + perfect_linear_sweep, magnitude_spectrum_weighted_sweep) from .stochastic import ( noise, pulsed_noise) @@ -21,4 +23,6 @@ __all__ = [ 'sine', 'impulse', 'noise', 'pulsed_noise', 'linear_sweep_time', 'exponential_sweep_time', + 'linear_sweep_freq', 'exponential_sweep_freq', + 'perfect_linear_sweep', 'magnitude_spectrum_weighted_sweep,' 'files'] diff --git a/pyfar/signals/deterministic.py b/pyfar/signals/deterministic.py index e72d7f788..ae5a06a1a 100644 --- a/pyfar/signals/deterministic.py +++ b/pyfar/signals/deterministic.py @@ -203,12 +203,20 @@ def linear_sweep_time(n_samples, frequency_range, n_fade_out=90, amplitude=1, def linear_sweep_freq( - n_samples, start_margin=None, stop_margin=None, frequency_range=None, - butterworth_order=8, double=True, sampling_rate=44100): + n_samples, frequency_range, start_margin, stop_margin, fade_in=0, + fade_out=0, butterworth_order=8, double=True, sampling_rate=44100): signal, group_delay = _frequency_domain_sweep( - n_samples, "linear", start_margin, stop_margin, - frequency_range, butterworth_order, double, sampling_rate) + n_samples=n_samples, + sweep_type='linear', + frequency_range=frequency_range, + butterworth_order=butterworth_order, + double=double, + start_margin=start_margin, + stop_margin=stop_margin, + fade_in=fade_in, + fade_out=fade_out, + sampling_rate=sampling_rate) return signal, group_delay @@ -283,39 +291,62 @@ def exponential_sweep_time(n_samples, frequency_range, n_fade_out=90, def exponential_sweep_freq( - n_samples, start_margin=None, stop_margin=None, frequency_range=None, - butterworth_order=8, double=True, sampling_rate=44100): + n_samples, frequency_range, start_margin, stop_margin, fade_in=0, + fade_out=0, butterworth_order=8, double=True, sampling_rate=44100): signal, group_delay = _frequency_domain_sweep( - n_samples, "exponential", start_margin, stop_margin, - frequency_range, butterworth_order, double, sampling_rate) + n_samples=n_samples, + sweep_type='exponential', + frequency_range=frequency_range, + butterworth_order=butterworth_order, + double=double, + start_margin=start_margin, + stop_margin=stop_margin, + fade_in=fade_in, + fade_out=fade_out, + sampling_rate=sampling_rate) return signal, group_delay def magnitude_spectrum_weighted_sweep( - n_samples, magnitude, start_margin=None, stop_margin=None, + n_samples, magnitude_spectrum, start_margin, stop_margin, double=True, sampling_rate=44100): signal, group_delay = _frequency_domain_sweep( - n_samples, magnitude, start_margin, stop_margin, - None, None, double, sampling_rate) + n_samples=n_samples, + sweep_type=magnitude_spectrum, + frequency_range=[0, sampling_rate / 2], + butterworth_order=0, + double=double, + start_margin=start_margin, + stop_margin=stop_margin, + fade_in=0, + fade_out=0, + sampling_rate=sampling_rate) return signal, group_delay -def perfect_linear_sweep( - n_samples, sampling_rate=44100): +def perfect_linear_sweep(n_samples, sampling_rate=44100): signal, group_delay = _frequency_domain_sweep( - n_samples, "perfect_linear", None, None, None, None, False, - sampling_rate) + n_samples=n_samples, + sweep_type='perfect_linear', + frequency_range=[0, sampling_rate / 2], + butterworth_order=0, + double=False, + start_margin=0, + stop_margin=0, + fade_in=0, + fade_out=0, + sampling_rate=sampling_rate) return signal, group_delay def _frequency_domain_sweep( - n_samples, magnitude, frequency_range, butterworth_order, double, + n_samples, sweep_type, frequency_range, butterworth_order, double, start_margin, stop_margin, fade_in, fade_out, sampling_rate): """ Frequency domain sweep synthesis with arbitrary magnitude response. @@ -340,7 +371,7 @@ def _frequency_domain_sweep( ---------- n_samples : int The length of the sweep in samples. - magnitude : Signal, string + sweep_type : Signal, string Specify the magnitude response of the sweep. signal @@ -363,26 +394,26 @@ def _frequency_domain_sweep( start_margin : int, float The time in samples, at which the sweep starts. The start margin is required because the frequency domain sweep synthesis has pre-ringing - in the time domain. Set to ``0`` if `magnitude` is + in the time domain. Set to ``0`` if `sweep_type` is ``'perfect_linear'``. stop_margin : int, float Time in samples, at which the sweep stops. This is relative to `n_samples`, e.g., a stop margin of 100 samples means that the sweep ends at sample ``n_samples-10``. This is required, because the frequency domain sweep synthesis has post-ringing in the time domain. - Set to ``0`` if `magnitude` is ``'perfect_linear'``. + Set to ``0`` if `sweep_type` is ``'perfect_linear'``. frequency_range : array like Frequency range of the sweep given by the lower and upper cut-off frequency in Hz. The restriction of the frequency range is realized by applying a Butterworth band-pass with the specified frequencies. - Ignored if `magnitude` is ``'perfect_linear'`` or `signal`. + Ignored if `sweep_type` is ``'perfect_linear'`` or `signal`. 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``. double : bool Double `n_samples` during the sweep calculation (recommended). Set to - ``False`` if `magnitude` is ``'perfect_linear'``. + ``False`` if `sweep_type` is ``'perfect_linear'``. sampling_rate : int The sampling rate in Hz. @@ -406,7 +437,7 @@ def _frequency_domain_sweep( ----- The envelope of the sweep time signal should be constant, appart from slight overshoots at the beginning and end. If this is not the case, try to - provide a smoother spectrum (if `magnitude` is `signal`) or increase + provide a smoother spectrum (if `sweep_type` is `signal`) or increase `n_samples`. References @@ -421,36 +452,38 @@ def _frequency_domain_sweep( Examples -------- - TODO Example with magnitude=singal + TODO Example with sweep_type=singal (e.g., Bass emphasis by means of low shelve filter) - TODO Examples with magnitude="linear" + TODO Examples with sweep_type="linear" - TODO Examples with magnitude="exponential" + TODO Examples with sweep_type="exponential" - TODO Examples with magnitude="perfect_linear" + TODO Examples with sweep_type="perfect_linear" """ # check input ------------------------------------------------------------- - if not isinstance(magnitude, (pyfar.Signal, str)): - raise TypeError("Magnitude must be type Signal or str.") - if isinstance(magnitude, str) and magnitude not in ['linear', - 'exponential', - 'perfect_linear']: - raise ValueError("Magnitude must be 'linear', 'exponential' or", - "'perfect_linear', when its type is str.") + if not isinstance(sweep_type, (pyfar.Signal, str)): + raise TypeError("sweep_type must be type Signal or str.") + if isinstance(sweep_type, pyfar.Signal): + magnitude = sweep_type + sweep_type = 'signal' + if sweep_type not in ['linear', 'exponential', 'perfect_linear', 'signal']: + raise ValueError("sweep_type must be 'linear', 'exponential' or", + "'perfect_linear', when it is a str.") if np.atleast_1d(frequency_range).size != 2: raise ValueError( "Frequency_range must be an array like with two elements.") if frequency_range[1] > sampling_rate/2: raise ValueError( "Upper frequency limit is larger than half the sampling rate.") - if frequency_range[0] == 0 and magnitude == "exponential": + if frequency_range[0] == 0 and sweep_type == "exponential": Warning(( "The exponential sweep has a 1/frequency magnitude spectrum. " "The magnitude is set to 0 at 0 Hz to avoid division by zero.")) - if magnitude == 'perfect_linear' and \ + if sweep_type == 'perfect_linear' and \ (start_margin != 0 or stop_margin != 0 or double or + fade_in != 0 or fade_out != 0 or frequency_range[0] != 0 or frequency_range[1] != sampling_rate / 2): # internal warning. Users will not call this function directly @@ -459,7 +492,7 @@ def _frequency_domain_sweep( # initialize basic parameters --------------------------------------------- # double n_samples - if double and magnitude != 'perfect_linear': + if double and sweep_type != 'perfect_linear': stop_margin += n_samples n_samples *= 2 @@ -470,25 +503,26 @@ def _frequency_domain_sweep( n_bins = n_samples // 2 + 1 # compute magnitude spectrum ---------------------------------------------- - if isinstance(magnitude, pyfar.Signal): + if sweep_type == 'signal': # zero pad magnitude Signal or raise error if needed if n_samples > magnitude.n_samples: magnitude = pyfar.dsp.pad_zeros( magnitude, n_samples - magnitude.n_samples) elif magnitude.n_samples > n_samples: - raise ValueError((f'magnitue can has {magnitude.n_samples} samples' - f' but must not be longer than {n_samples}')) - sweep_abs = np.abs(magnitude.freq_raw) - elif magnitude in ['linear', 'perfect_linear']: + raise ValueError( + (f'magnitude_spectrum has {magnitude.n_samples} samples ' + f'but must not be longer than {n_samples}')) + sweep_abs = np.abs(magnitude.freq_raw.flatten()) + elif sweep_type in ['linear', 'perfect_linear']: # constant spectrum sweep_abs = np.ones(n_bins) - elif magnitude == 'exponential': + elif sweep_type == 'exponential': # 1/f spectrum sweep_abs = np.zeros(n_bins) sweep_abs[1:] = 1 / np.sqrt(2 * np.pi * np.arange(1, n_bins) * df) # band limit to magnitude spectrum - if magnitude in ['linear', 'exponential']: + if sweep_type in ['linear', 'exponential']: if frequency_range[0] > 0 and frequency_range[1] < sampling_rate / 2: band_limit = pyfar.dsp.filter.butterworth( pyfar.signals.impulse(n_samples, sampling_rate=sampling_rate), @@ -533,7 +567,7 @@ def _frequency_domain_sweep( sweep_ang = pyfar.dsp.wrap_to_2pi(sweep_ang) sweep_ang[sweep_ang > np.pi] -= 2*np.pi - if magnitude == 'perfect_linear': + if sweep_type == 'perfect_linear': sweep_ang[-1] = 0 elif sweep_ang[-1] != 0 and not n_samples % 2: factor = np.cumsum(np.ones_like(sweep_ang)) - 1 @@ -583,7 +617,7 @@ def _frequency_domain_sweep( # (to avoid clipping if written to fixed point wav file) sweep = pyfar.dsp.normalize(sweep) * (1 - 2**-15) - return sweep, sweep_gd, sweep_abs + return sweep, sweep_gd def _time_domain_sweep(n_samples, frequency_range, n_fade_out, amplitude, From 717cafdf24f444cd2feace6e128e5850002d95a1 Mon Sep 17 00:00:00 2001 From: Fabian Brinkmann Date: Fri, 20 Oct 2023 12:17:05 +0200 Subject: [PATCH 20/42] fix flake8 error --- pyfar/signals/__init__.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pyfar/signals/__init__.py b/pyfar/signals/__init__.py index 460973d47..2a3a002a0 100644 --- a/pyfar/signals/__init__.py +++ b/pyfar/signals/__init__.py @@ -24,5 +24,5 @@ 'sine', 'impulse', 'noise', 'pulsed_noise', 'linear_sweep_time', 'exponential_sweep_time', 'linear_sweep_freq', 'exponential_sweep_freq', - 'perfect_linear_sweep', 'magnitude_spectrum_weighted_sweep,' + 'perfect_linear_sweep', 'magnitude_spectrum_weighted_sweep', 'files'] From e5d7c2e3f1e8378fd8840515e0ff508ccb24cc86 Mon Sep 17 00:00:00 2001 From: Fabian Brinkmann Date: Thu, 26 Oct 2023 12:28:45 +0200 Subject: [PATCH 21/42] bug fix --- pyfar/signals/deterministic.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/pyfar/signals/deterministic.py b/pyfar/signals/deterministic.py index ae5a06a1a..e7c9123e4 100644 --- a/pyfar/signals/deterministic.py +++ b/pyfar/signals/deterministic.py @@ -580,7 +580,8 @@ def _frequency_domain_sweep( sweep = sweep_abs * np.exp(1j * sweep_ang) # put sweep in pyfar.Signal an transform to time domain - sweep = pyfar.Signal(sweep, sampling_rate, n_samples, 'freq', 'rms') + sweep = pyfar.Signal(sweep, sampling_rate, n_samples, 'freq', 'none') + sweep.fft_norm = 'rms' # put group delay on pyfar FrequencyData sweep_gd = pyfar.FrequencyData( From d7a5fb65073dc2e05b29b1e2bac6a75cf76a9ab9 Mon Sep 17 00:00:00 2001 From: Fabian Brinkmann Date: Thu, 26 Oct 2023 13:13:20 +0200 Subject: [PATCH 22/42] improve windowing --- pyfar/signals/deterministic.py | 32 +++++++++++++++++++------------- 1 file changed, 19 insertions(+), 13 deletions(-) diff --git a/pyfar/signals/deterministic.py b/pyfar/signals/deterministic.py index e7c9123e4..80f2f6249 100644 --- a/pyfar/signals/deterministic.py +++ b/pyfar/signals/deterministic.py @@ -576,6 +576,10 @@ def _frequency_domain_sweep( sweep_ang[-1] = np.abs(sweep_ang[-1]) # compute and finalize return data ---------------------------------------- + # put group delay on pyfar FrequencyData + sweep_gd = pyfar.FrequencyData( + sweep_gd, pyfar.dsp.fft.rfftfreq(n_samples, sampling_rate)) + # combine magnitude and phase of sweep sweep = sweep_abs * np.exp(1j * sweep_ang) @@ -583,34 +587,36 @@ def _frequency_domain_sweep( sweep = pyfar.Signal(sweep, sampling_rate, n_samples, 'freq', 'none') sweep.fft_norm = 'rms' - # put group delay on pyfar FrequencyData - sweep_gd = pyfar.FrequencyData( - sweep_gd, pyfar.dsp.fft.rfftfreq(n_samples, sampling_rate)) - # cut to originally desired length if double: n_samples = n_samples // 2 stop_margin -= n_samples sweep.time = sweep.time[..., :n_samples] - # apply window + # find windowing end points heuristically. The start and stop margin are + # not reliable, because freq. domain synthesis causes pre/post-ringing + if fade_in or fade_out: + threshold = 10**(-60/20) * np.max(np.abs(sweep.time)) + fade_start = np.argmax(np.abs(sweep.time) > threshold) + fade_end = n_samples - \ + np.argmax(np.abs(sweep.time.flatten()[::-1]) > threshold) + print([fade_start, fade_end]) + + # generate windows for fade in and/or out if fade_in and fade_out: - fade = [start_margin, - start_margin + fade_in, - n_samples - stop_margin - fade_out, - n_samples - stop_margin] + fade = [fade_start, fade_start + fade_in, + fade_end - fade_out, fade_end] shape = 'symmetric' elif fade_in: - fade = [start_margin, - start_margin + fade_in] + fade = [fade_start, fade_start + fade_in] shape = 'left' elif fade_out: - fade = [n_samples - stop_margin - fade_out, - n_samples - stop_margin] + fade = [fade_end - fade_out, fade_end] shape = 'right' else: fade = None + # apply window if fade is not None: sweep = pyfar.dsp.time_window(sweep, fade, shape=shape) From f0fb9e92abe41326490e562932095eb6f4dfbb96 Mon Sep 17 00:00:00 2001 From: Fabian Brinkmann Date: Thu, 26 Oct 2023 14:13:03 +0200 Subject: [PATCH 23/42] [skip ci] initial version of docstrings --- pyfar/signals/deterministic.py | 315 +++++++++++++++++++++++++++------ 1 file changed, 259 insertions(+), 56 deletions(-) diff --git a/pyfar/signals/deterministic.py b/pyfar/signals/deterministic.py index 80f2f6249..592324bef 100644 --- a/pyfar/signals/deterministic.py +++ b/pyfar/signals/deterministic.py @@ -158,7 +158,7 @@ def linear_sweep_time(n_samples, frequency_range, n_fade_out=90, amplitude=1, .. note:: The linear sweep can also be generated in the frequency domain - (see :py:func:`~general_sweep_synthesis`). Time domain synthesis + (see :py:func:`~exponential_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 of a slightly irregular temporal @@ -205,6 +205,82 @@ def linear_sweep_time(n_samples, frequency_range, n_fade_out=90, amplitude=1, def linear_sweep_freq( n_samples, frequency_range, start_margin, stop_margin, fade_in=0, fade_out=0, butterworth_order=8, double=True, sampling_rate=44100): + """ + Generate single channel sine sweep with linearlly increasing frequency. + + Sweep sweep synthesis according to [#]_ + + .. note:: + The linear sweep can also be generated in the time domain + (:py:func:`~linear_sweep_time`). Frequency domain synthesis exhibits + smooth magnitude spectra in trade of a slightly irregular temporal + envelope. Time domain synthesis exhibits a constant temporal envelope + in trade of slight ripples in the magnitude response. + + Parameters + ---------- + n_samples : int + The length of the sweep in samples. + frequency_range : array like + Frequency range of the sweep given by the lower and upper cut-off + frequency in Hz. The restriction of the frequency range is realized + by applying a Butterworth band-pass with the specified frequencies. + start_margin : int, float + The time in samples, at which the sweep starts. The start margin is + required because the frequency domain sweep synthesis has pre-ringing + in the time domain. + stop_margin : int, float + Time in samples, at which the sweep stops. This is relative to + `n_samples`, e.g., a stop margin of 100 samples means that the sweep + ends at sample ``n_samples-10``. This is required, because the + frequency domain sweep synthesis has post-ringing in the time domain. + 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 + a fade-in. + 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 + a fade-out. + 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``. + double : bool + Double `n_samples` during the sweep calculation. The default is + ``True``. + sampling_rate : int + The sampling rate in Hz. The default is ``44100``. + + Returns + ------- + sweep : Signal + The sweep signal. The Signal is in the time domain, has a maximum + absolute amplitude of 1 and the ``rms`` FFT normalization + (see :py:func:`~pyfar.dsp.fft.normalization`). + group_delay_sweep : FrequencyData + The group delay of the sweep in seconds as a single sided spectrum. + + Notes + ----- + The envelope of the sweep time signal should be constant, apart from + slight overshoots at the beginning and end. If this is not the case, try + to increase `n_samples`, `start_margin`, `stop_margin`, `fade_in` or + `fade_out`. + + References + ---------- + .. [#] 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. + + Examples + -------- + TODO + """ signal, group_delay = _frequency_domain_sweep( n_samples=n_samples, @@ -242,7 +318,7 @@ def exponential_sweep_time(n_samples, frequency_range, n_fade_out=90, .. note:: The exponential sweep can also be generated in the frequency domain - (see :py:func:`~general_sweep_synthesis`). Time domain synthesis + (see :py:func:`~exponential_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 of a slightly irregular temporal @@ -293,6 +369,82 @@ def exponential_sweep_time(n_samples, frequency_range, n_fade_out=90, def exponential_sweep_freq( n_samples, frequency_range, start_margin, stop_margin, fade_in=0, fade_out=0, butterworth_order=8, double=True, sampling_rate=44100): + """ + Generate single channel sine sweep with exponentially increasing frequency. + + Sweep sweep synthesis according to [#]_ + + .. note:: + 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 + temporal envelope. Time domain synthesis exhibits a constant temporal + envelope in trade of slight ripples in the magnitude response. + + Parameters + ---------- + n_samples : int + The length of the sweep in samples. + frequency_range : array like + Frequency range of the sweep given by the lower and upper cut-off + frequency in Hz. The restriction of the frequency range is realized + by applying a Butterworth band-pass with the specified frequencies. + start_margin : int, float + The time in samples, at which the sweep starts. The start margin is + required because the frequency domain sweep synthesis has pre-ringing + in the time domain. + stop_margin : int, float + Time in samples, at which the sweep stops. This is relative to + `n_samples`, e.g., a stop margin of 100 samples means that the sweep + ends at sample ``n_samples-10``. This is required, because the + frequency domain sweep synthesis has post-ringing in the time domain. + 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 + a fade-in. + 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 + a fade-out. + 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``. + double : bool + Double `n_samples` during the sweep calculation. The default is + ``True``. + sampling_rate : int + The sampling rate in Hz. The default is ``44100``. + + Returns + ------- + sweep : Signal + The sweep signal. The Signal is in the time domain, has a maximum + absolute amplitude of 1 and the ``rms`` FFT normalization + (see :py:func:`~pyfar.dsp.fft.normalization`). + group_delay_sweep : FrequencyData + The group delay of the sweep in seconds as a single sided spectrum. + + Notes + ----- + The envelope of the sweep time signal should be constant, apart from + slight overshoots at the beginning and end. If this is not the case, try + to increase `n_samples`, `start_margin`, `stop_margin`, `fade_in` or + `fade_out`. + + References + ---------- + .. [#] 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. + + Examples + -------- + TODO + """ signal, group_delay = _frequency_domain_sweep( n_samples=n_samples, @@ -311,7 +463,77 @@ def exponential_sweep_freq( def magnitude_spectrum_weighted_sweep( n_samples, magnitude_spectrum, start_margin, stop_margin, - double=True, sampling_rate=44100): + fade_in, fade_out, double=True, sampling_rate=44100): + """ + Generate single channel sine sweep with arbitrary magnitude spectrum. + + Sweep sweep synthesis according to [#]_ + + .. note:: + The linear sweep can also be generated in the time domain + (:py:func:`~linear_sweep_time`). Frequency domain synthesis exhibits + smooth magnitude spectra in trade of a slightly irregular temporal + envelope. Time domain synthesis exhibits a constant temporal envelope + in trade of slight ripples in the magnitude response. + + Parameters + ---------- + n_samples : int + The length of the sweep in samples. + magnitude_spectrum : Signal + The magnitude spectrum as a Signal object. Must be at least `n_samples` + long and is zero-padded to `n_samples` otherwise. + start_margin : int, float + The time in samples, at which the sweep starts. The start margin is + required because the frequency domain sweep synthesis has pre-ringing + in the time domain. + stop_margin : int, float + Time in samples, at which the sweep stops. This is relative to + `n_samples`, e.g., a stop margin of 100 samples means that the sweep + ends at sample ``n_samples-10``. This is required, because the + frequency domain sweep synthesis has post-ringing in the time domain. + 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 + a fade-in. + 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 + a fade-out. + double : bool + Double `n_samples` during the sweep calculation. The default is + ``True``. + sampling_rate : int + The sampling rate in Hz. The default is ``44100``. + + Returns + ------- + sweep : Signal + The sweep signal. The Signal is in the time domain, has a maximum + absolute amplitude of 1 and the ``rms`` FFT normalization + (see :py:func:`~pyfar.dsp.fft.normalization`). + group_delay_sweep : FrequencyData + The group delay of the sweep in seconds as a single sided spectrum. + + Notes + ----- + The envelope of the sweep time signal should be constant, apart from + slight overshoots at the beginning and end. If this is not the case, + provide a smother spectrum in `magnitude_spectrum` or try to increase + `n_samples`, `start_margin`, `stop_margin`, `fade_in` or `fade_out`. + + References + ---------- + .. [#] 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. + + Examples + -------- + TODO + """ signal, group_delay = _frequency_domain_sweep( n_samples=n_samples, @@ -321,14 +543,45 @@ def magnitude_spectrum_weighted_sweep( double=double, start_margin=start_margin, stop_margin=stop_margin, - fade_in=0, - fade_out=0, + fade_in=fade_in, + fade_out=fade_out, sampling_rate=sampling_rate) return signal, group_delay -def perfect_linear_sweep(n_samples, sampling_rate=44100): +def linear_perfect_sweep(n_samples, sampling_rate=44100): + """ + Generate a perfect linear sweep. + + The perfect sweep is generated according to [#]_ and is used for adaptive + system identification. It is orthogonal to delayed versions of itself and + can be played back in a loop. + + Parameters + ---------- + n_samples : int + The length of the sweep in samples. + sampling_rate : int + The sampling rate in Hz. The default is ``44100``. + + Returns + ------- + sweep : Signal + The sweep signal. The Signal is in the time domain, has a maximum + absolute amplitude of 1 and the ``rms`` FFT normalization + (see :py:func:`~pyfar.dsp.fft.normalization`). + group_delay_sweep : FrequencyData + The group delay of the sweep in seconds as a single sided spectrum. + + References + ---------- + [#] C. Antweiler, A. Telle, P. Vary, and G. Enzner, “Perfect-sweep NLMS for + time-variant acoustic system identification,” in IEEE International + Conference on Acoustics, Speech and Signal Processing (ICASSP), Prague, + Czech Republic, May 2011. doi: 10.1109/ICASSP.2012.6287930. + + """ signal, group_delay = _frequency_domain_sweep( n_samples=n_samples, @@ -351,22 +604,6 @@ def _frequency_domain_sweep( """ Frequency domain sweep synthesis with arbitrary magnitude response. - TODO Link to fractional octave smoothing in Notes - - TODO Implement calculation of group delay of impulse responses - - TODO Examples - - Sweep sweep synthesis according to [#]_ - - .. note:: - The linear and exponential sweep can also be generated in the time - domain (see :py:func:`~linear_sweep`, :py:func:`~exponential_sweep`). - Frequency domain synthesis exhibits smooth magnitude spectra and in - trade of a slightly irregular temporal envelope. Time domain synthesis - exhibits a constant temporal envelope in trade of slight ripples in the - magnitude response. - Parameters ---------- n_samples : int @@ -426,40 +663,6 @@ def _frequency_domain_sweep( group_delay_sweep : FrequencyData The group delay of the sweep as a single sided spectrum between 0 Hz and ``sampling_rate/2``. - - TODO add this after implementation is complete: - - This can be used to calculate the group delay of the impulse responses - of linear and harmonic distortion products after deconvoloution (see - :py:func:`~pyfar.dsp...`). - - Notes - ----- - The envelope of the sweep time signal should be constant, appart from - slight overshoots at the beginning and end. If this is not the case, try to - provide a smoother spectrum (if `sweep_type` is `signal`) or increase - `n_samples`. - - References - ---------- - .. [#] 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. - .. [#] C. Antweiler, A. Telle, P. Vary, G. Enzner. 'Perfect-Sweep NLMS for - Time-Variant Acoustic System Identification,' IEEE Int. Conf. - Acoustics, Speech and Signal Processing (ICASSP), - Prague, Czech Republic, 2011. doi: 10.1109/ICASSP.2012.6287930. - - Examples - -------- - TODO Example with sweep_type=singal - (e.g., Bass emphasis by means of low shelve filter) - - TODO Examples with sweep_type="linear" - - TODO Examples with sweep_type="exponential" - - TODO Examples with sweep_type="perfect_linear" """ # check input ------------------------------------------------------------- From 67e64037c1a89270029ab231f29bce84c29fa859 Mon Sep 17 00:00:00 2001 From: Fabian Brinkmann Date: Thu, 26 Oct 2023 14:14:59 +0200 Subject: [PATCH 24/42] [skip ci] improve docstring --- pyfar/signals/deterministic.py | 2 ++ 1 file changed, 2 insertions(+) diff --git a/pyfar/signals/deterministic.py b/pyfar/signals/deterministic.py index 592324bef..98ac5baf2 100644 --- a/pyfar/signals/deterministic.py +++ b/pyfar/signals/deterministic.py @@ -389,6 +389,8 @@ def exponential_sweep_freq( Frequency range of the sweep given by the lower and upper cut-off frequency in Hz. The restriction of the frequency range is realized by applying a Butterworth band-pass with the specified frequencies. + Note that the exponential sweep can not start at 0 Hz, because its + magnitude is defined by 1/frequency. start_margin : int, float The time in samples, at which the sweep starts. The start margin is required because the frequency domain sweep synthesis has pre-ringing From 788050077d684eef94c663049ae9f5307e909f57 Mon Sep 17 00:00:00 2001 From: Fabian Brinkmann Date: Fri, 27 Oct 2023 15:58:27 +0200 Subject: [PATCH 25/42] add plot examples --- pyfar/signals/__init__.py | 4 +- pyfar/signals/deterministic.py | 126 +++++++++++++++++++++++++++------ 2 files changed, 107 insertions(+), 23 deletions(-) diff --git a/pyfar/signals/__init__.py b/pyfar/signals/__init__.py index 2a3a002a0..73fefd866 100644 --- a/pyfar/signals/__init__.py +++ b/pyfar/signals/__init__.py @@ -13,7 +13,7 @@ from .deterministic import ( sine, impulse, linear_sweep_time, exponential_sweep_time, linear_sweep_freq, exponential_sweep_freq, - perfect_linear_sweep, magnitude_spectrum_weighted_sweep) + linear_perfect_sweep, magnitude_spectrum_weighted_sweep) from .stochastic import ( noise, pulsed_noise) @@ -24,5 +24,5 @@ 'sine', 'impulse', 'noise', 'pulsed_noise', 'linear_sweep_time', 'exponential_sweep_time', 'linear_sweep_freq', 'exponential_sweep_freq', - 'perfect_linear_sweep', 'magnitude_spectrum_weighted_sweep', + 'linear_perfect_sweep', 'magnitude_spectrum_weighted_sweep', 'files'] diff --git a/pyfar/signals/deterministic.py b/pyfar/signals/deterministic.py index 98ac5baf2..4b27a6fa8 100644 --- a/pyfar/signals/deterministic.py +++ b/pyfar/signals/deterministic.py @@ -193,6 +193,16 @@ def linear_sweep_time(n_samples, frequency_range, n_fade_out=90, amplitude=1, .. [#] Farina, Angelo (2000): "Simultaneous measurement of impulse response and distortion with a swept-sine technique." 108th AES Convention, Paris: France. + + Examples + -------- + Linear sweep between 0 and 22050 Hz + + .. plot:: + + >>> import pyfar as pf + >>> sweep = pf.signals.linear_sweep_time(2**16, [0, 22050]) + >>> pf.plot.time_freq(sweep) """ signal = _time_domain_sweep( @@ -203,8 +213,8 @@ def linear_sweep_time(n_samples, frequency_range, n_fade_out=90, amplitude=1, def linear_sweep_freq( - n_samples, frequency_range, start_margin, stop_margin, fade_in=0, - fade_out=0, butterworth_order=8, double=True, sampling_rate=44100): + n_samples, frequency_range, start_margin, stop_margin, fade_in=None, + fade_out=None, butterworth_order=8, double=True, sampling_rate=44100): """ Generate single channel sine sweep with linearlly increasing frequency. @@ -275,11 +285,18 @@ def linear_sweep_freq( ---------- .. [#] 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. + Corrections. J. Audio Eng. Soc. 2001, 49 (6), 443-471. Examples -------- - TODO + Linear sweep between 0 and 22050 Hz + + .. plot:: + + >>> import pyfar as pf + >>> sweep, _ = pf.signals.linear_sweep_freq( + ... 2**16, [0, 22050], 1000, 1000) + >>> pf.plot.time_freq(sweep) """ signal, group_delay = _frequency_domain_sweep( @@ -357,6 +374,17 @@ def exponential_sweep_time(n_samples, frequency_range, n_fade_out=90, .. [#] Farina, Angelo (2000): "Simultaneous measurement of impulse response and distortion with a swept-sine technique." 108th AES Convention, Paris: France. + + Examples + -------- + + Exponential sweep between 50 and 22050 Hz + + .. plot:: + + >>> import pyfar as pf + >>> sweep = pf.signals.exponential_sweep_time(2**16, [50, 22050]) + >>> pf.plot.time_freq(sweep) """ signal = _time_domain_sweep( @@ -367,8 +395,8 @@ def exponential_sweep_time(n_samples, frequency_range, n_fade_out=90, def exponential_sweep_freq( - n_samples, frequency_range, start_margin, stop_margin, fade_in=0, - fade_out=0, butterworth_order=8, double=True, sampling_rate=44100): + n_samples, frequency_range, start_margin, stop_margin, fade_in=None, + fade_out=None, butterworth_order=8, double=True, sampling_rate=44100): """ Generate single channel sine sweep with exponentially increasing frequency. @@ -441,11 +469,19 @@ def exponential_sweep_freq( ---------- .. [#] 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. + Corrections. J. Audio Eng. Soc. 2001, 49 (6), 443-471. Examples -------- - TODO + + Exponential sweep between 50 and 22050 Hz + + .. plot:: + + >>> import pyfar as pf + >>> sweep, _ = pf.signals.exponential_sweep_freq( + ... 2**16, [50, 22050], 5000, 1000) + >>> pf.plot.time_freq(sweep) """ signal, group_delay = _frequency_domain_sweep( @@ -465,7 +501,7 @@ def exponential_sweep_freq( def magnitude_spectrum_weighted_sweep( n_samples, magnitude_spectrum, start_margin, stop_margin, - fade_in, fade_out, double=True, sampling_rate=44100): + fade_in=None, fade_out=None, double=True, sampling_rate=44100): """ Generate single channel sine sweep with arbitrary magnitude spectrum. @@ -530,11 +566,20 @@ def magnitude_spectrum_weighted_sweep( ---------- .. [#] 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. + Corrections. J. Audio Eng. Soc. 2001, 49 (6), 443-471. Examples -------- - TODO + + .. plot:: + + >>> 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) """ signal, group_delay = _frequency_domain_sweep( @@ -578,11 +623,50 @@ def linear_perfect_sweep(n_samples, sampling_rate=44100): References ---------- - [#] C. Antweiler, A. Telle, P. Vary, and G. Enzner, “Perfect-sweep NLMS for - time-variant acoustic system identification,” in IEEE International - Conference on Acoustics, Speech and Signal Processing (ICASSP), Prague, - Czech Republic, May 2011. doi: 10.1109/ICASSP.2012.6287930. + .. [#] C. Antweiler, A. Telle, P. Vary, and G. Enzner, “Perfect-sweep NLMS + for time-variant acoustic system identification,” in IEEE + Int. Conf. Acoustics, Speech and Signal Processing (ICASSP), Prague, + Czech Republic, May 2011. doi: 10.1109/ICASSP.2012.6287930. + Examples + -------- + Plot a shifted perfect sweep to show that it can be looped + + .. plot:: + + >>> import pyfar as pf + >>> sweep, _ = pf.signals.linear_perfect_sweep(2**8) + >>> sweep_shifted = pf.dsp.time_shift(sweep, 2**7) + >>> ax = pf.plot.time_freq(sweep_shifted, unit='samples') + >>> ax[0].axvline(2**7, ls='--', c='r', label='start/end') + >>> ax[0].legend(loc=4) + + Show that the sweep is orthogonal to all shifted version of itself. This + property is important for adaptive system identification. + + .. plot:: + + >>> import pyfar as pf + >>> import numpy as np + >>> import matplotlib.pyplot as plt + >>> + >>> sweep, _ = pf.signals.linear_perfect_sweep(2**8) + >>> + >>> # compute auto-correlation + >>> 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) + >>> + >>> # 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') """ signal, group_delay = _frequency_domain_sweep( @@ -593,8 +677,8 @@ def linear_perfect_sweep(n_samples, sampling_rate=44100): double=False, start_margin=0, stop_margin=0, - fade_in=0, - fade_out=0, + fade_in=None, + fade_out=None, sampling_rate=sampling_rate) return signal, group_delay @@ -626,9 +710,9 @@ def _frequency_domain_sweep( magnitude decreases by 3 dB per frequency doubling and has constant energy in fiters of relative constant bandwidth (e.g. octaves). ``'perfect_linear'`` - Perfect linear sweep according to [#]_. Note that the parameters - `start_margin`, `stop_margin`, `frequency_range` and `double` are - not required in this case. + Perfect linear sweep. Note that the parameters `start_margin`, + `stop_margin`, `frequency_range` and `double` are not required in + this case. start_margin : int, float The time in samples, at which the sweep starts. The start margin is @@ -688,7 +772,7 @@ def _frequency_domain_sweep( "The magnitude is set to 0 at 0 Hz to avoid division by zero.")) if sweep_type == 'perfect_linear' and \ (start_margin != 0 or stop_margin != 0 or double or - fade_in != 0 or fade_out != 0 or + fade_in is not None or fade_out is not None or frequency_range[0] != 0 or frequency_range[1] != sampling_rate / 2): # internal warning. Users will not call this function directly From 1f4f97493d61faefd7273dcd72745bf3fb2c7bcc Mon Sep 17 00:00:00 2001 From: Fabian Brinkmann Date: Sat, 28 Oct 2023 09:24:30 +0200 Subject: [PATCH 26/42] add tests for linear perfect sweep --- tests/test_signals.py | 40 +++++++++++++++++++++++++++++++++++----- 1 file changed, 35 insertions(+), 5 deletions(-) diff --git a/tests/test_signals.py b/tests/test_signals.py index a141f8817..fa4349739 100644 --- a/tests/test_signals.py +++ b/tests/test_signals.py @@ -2,7 +2,7 @@ import numpy.testing as npt import pytest import os -from pyfar import Signal +import pyfar as pf import pyfar.dsp.filter as pff import pyfar.signals as pfs from pyfar.signals.deterministic import _match_shape @@ -14,7 +14,7 @@ def test_sine_with_defaults(): signal = pfs.sine(99, 441) sin = np.sin(np.arange(441) / 44100 * 2 * np.pi * 99) - assert isinstance(signal, Signal) + assert isinstance(signal, pf.Signal) assert signal.comment == ("Sine signal (f = [99] Hz, amplitude = [1], " "phase = [0] rad)") assert signal.sampling_rate == 44100 @@ -70,7 +70,7 @@ def test_sine_assertions(): def test_impulse_with_defaults(): """Test impulse with default parameters""" signal = pfs.impulse(3) - assert isinstance(signal, Signal) + assert isinstance(signal, pf.Signal) npt.assert_allclose(signal.time, np.atleast_2d([1, 0, 0])) assert signal.sampling_rate == 44100 assert signal.fft_norm == 'none' @@ -113,7 +113,7 @@ def test_noise_with_defaults(): """Test noise with default parameters.""" signal = pfs.noise(100) - assert isinstance(signal, Signal) + assert isinstance(signal, pf.Signal) assert signal.sampling_rate == 44100 assert signal.fft_norm == "rms" assert signal.comment == "white noise signal (rms = [1])" @@ -178,7 +178,7 @@ def test_pulsed_noise_with_defaults(): """Test pulsed noise signal generation with default values.""" signal = pfs.pulsed_noise(n_pulse=200, n_pause=100) - assert isinstance(signal, Signal) + assert isinstance(signal, pf.Signal) assert signal.sampling_rate == 44100 assert signal.fft_norm == "rms" assert signal.n_samples == 5 * 200 + 4 * 100 @@ -352,3 +352,33 @@ def test_match_shape(): assert cshape == (2, 3) npt.assert_allclose(np.ones(cshape), a_match) npt.assert_allclose(b, b_match) + + +@pytest.mark.parametrize('n_samples', [2**8, 2**8 + 1]) +@pytest.mark.parametrize('sampling_rate', [44100, 48000]) +def test_linear_perfect_sweep(n_samples, sampling_rate): + '''Test the perfect sweep generation''' + + sweep, group_delay = pfs.linear_perfect_sweep(n_samples, sampling_rate) + + # basic checks + assert type(sweep) is pf.Signal + assert sweep.n_samples == n_samples + assert sweep.sampling_rate == sampling_rate + assert type(group_delay) is pf.FrequencyData + + # assert magnitude response + npt.assert_almost_equal(np.abs(sweep.freq_raw), + np.mean(np.abs(sweep.freq_raw))) + + # compute normalized auto-correlation + auto_correlation = np.empty(n_samples-1) + for shift in range(1, n_samples): + auto_correlation[shift - 1] = np.dot( + sweep.time.flatten(), + np.roll(sweep.time.flatten(), shift)) + auto_correlation /= pf.dsp.energy(sweep) + + # must be 0 - sweep is orthogonal to shifted versions of itself + assert np.all( + np.abs(auto_correlation) < 5 * np.finfo(sweep.time.dtype).eps) From b46bd7dd8119f6ba6fabac6847176e547b0bc5c9 Mon Sep 17 00:00:00 2001 From: Fabian Brinkmann Date: Sat, 28 Oct 2023 14:09:30 +0200 Subject: [PATCH 27/42] add testing for public functions --- pyfar/signals/deterministic.py | 61 +++++++---------------- tests/test_signals.py | 90 +++++++++++++++++++++++++++++++++- 2 files changed, 106 insertions(+), 45 deletions(-) diff --git a/pyfar/signals/deterministic.py b/pyfar/signals/deterministic.py index 4b27a6fa8..74a57a1cc 100644 --- a/pyfar/signals/deterministic.py +++ b/pyfar/signals/deterministic.py @@ -214,7 +214,7 @@ def linear_sweep_time(n_samples, frequency_range, n_fade_out=90, amplitude=1, def linear_sweep_freq( n_samples, frequency_range, start_margin, stop_margin, fade_in=None, - fade_out=None, butterworth_order=8, double=True, sampling_rate=44100): + fade_out=None, butterworth_order=8, sampling_rate=44100): """ Generate single channel sine sweep with linearlly increasing frequency. @@ -259,9 +259,6 @@ def linear_sweep_freq( 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``. - double : bool - Double `n_samples` during the sweep calculation. The default is - ``True``. sampling_rate : int The sampling rate in Hz. The default is ``44100``. @@ -304,7 +301,7 @@ def linear_sweep_freq( sweep_type='linear', frequency_range=frequency_range, butterworth_order=butterworth_order, - double=double, + double=False, start_margin=start_margin, stop_margin=stop_margin, fade_in=fade_in, @@ -396,7 +393,7 @@ def exponential_sweep_time(n_samples, frequency_range, n_fade_out=90, def exponential_sweep_freq( n_samples, frequency_range, start_margin, stop_margin, fade_in=None, - fade_out=None, butterworth_order=8, double=True, sampling_rate=44100): + fade_out=None, butterworth_order=8, sampling_rate=44100): """ Generate single channel sine sweep with exponentially increasing frequency. @@ -443,9 +440,6 @@ def exponential_sweep_freq( 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``. - double : bool - Double `n_samples` during the sweep calculation. The default is - ``True``. sampling_rate : int The sampling rate in Hz. The default is ``44100``. @@ -489,7 +483,7 @@ def exponential_sweep_freq( sweep_type='exponential', frequency_range=frequency_range, butterworth_order=butterworth_order, - double=double, + double=False, start_margin=start_margin, stop_margin=stop_margin, fade_in=fade_in, @@ -501,7 +495,7 @@ def exponential_sweep_freq( def magnitude_spectrum_weighted_sweep( n_samples, magnitude_spectrum, start_margin, stop_margin, - fade_in=None, fade_out=None, double=True, sampling_rate=44100): + fade_in=None, fade_out=None, sampling_rate=44100): """ Generate single channel sine sweep with arbitrary magnitude spectrum. @@ -540,9 +534,6 @@ def magnitude_spectrum_weighted_sweep( 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 a fade-out. - double : bool - Double `n_samples` during the sweep calculation. The default is - ``True``. sampling_rate : int The sampling rate in Hz. The default is ``44100``. @@ -587,7 +578,7 @@ def magnitude_spectrum_weighted_sweep( sweep_type=magnitude_spectrum, frequency_range=[0, sampling_rate / 2], butterworth_order=0, - double=double, + double=False, start_margin=start_margin, stop_margin=stop_margin, fade_in=fade_in, @@ -690,10 +681,12 @@ def _frequency_domain_sweep( """ Frequency domain sweep synthesis with arbitrary magnitude response. + Non-unique parameters are documented in ``linear_sweep_freq``, + ``exponential_sweep_freq``, ``linear_perfect_sweep``, and + ``magnitude_spectrum_weighted_sweep``. + Parameters ---------- - n_samples : int - The length of the sweep in samples. sweep_type : Signal, string Specify the magnitude response of the sweep. @@ -708,37 +701,17 @@ def _frequency_domain_sweep( ``'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). + energy in filters of relative constant bandwidth (e.g. octaves). ``'perfect_linear'`` Perfect linear sweep. Note that the parameters `start_margin`, `stop_margin`, `frequency_range` and `double` are not required in this case. - - start_margin : int, float - The time in samples, at which the sweep starts. The start margin is - required because the frequency domain sweep synthesis has pre-ringing - in the time domain. Set to ``0`` if `sweep_type` is - ``'perfect_linear'``. - stop_margin : int, float - Time in samples, at which the sweep stops. This is relative to - `n_samples`, e.g., a stop margin of 100 samples means that the sweep - ends at sample ``n_samples-10``. This is required, because the - frequency domain sweep synthesis has post-ringing in the time domain. - Set to ``0`` if `sweep_type` is ``'perfect_linear'``. - frequency_range : array like - Frequency range of the sweep given by the lower and upper cut-off - frequency in Hz. The restriction of the frequency range is realized - by applying a Butterworth band-pass with the specified frequencies. - Ignored if `sweep_type` is ``'perfect_linear'`` or `signal`. - 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``. double : bool - Double `n_samples` during the sweep calculation (recommended). Set to - ``False`` if `sweep_type` is ``'perfect_linear'``. - sampling_rate : int - The sampling rate in Hz. + Double `n_samples` during the sweep calculation. This parameter seems + not to be required in the current implementation. It was kept if cases + come up, where it could help. It is contained in the AKtools version + from which this function was ported. Must be ``False`` if `sweep_type` + is ``'perfect_linear'``. Returns ------- @@ -806,7 +779,7 @@ def _frequency_domain_sweep( # constant spectrum sweep_abs = np.ones(n_bins) elif sweep_type == 'exponential': - # 1/f spectrum + # 1/sqrt(f) spectrum sweep_abs = np.zeros(n_bins) sweep_abs[1:] = 1 / np.sqrt(2 * np.pi * np.arange(1, n_bins) * df) diff --git a/tests/test_signals.py b/tests/test_signals.py index fa4349739..2a0d9d8f7 100644 --- a/tests/test_signals.py +++ b/tests/test_signals.py @@ -357,7 +357,10 @@ def test_match_shape(): @pytest.mark.parametrize('n_samples', [2**8, 2**8 + 1]) @pytest.mark.parametrize('sampling_rate', [44100, 48000]) def test_linear_perfect_sweep(n_samples, sampling_rate): - '''Test the perfect sweep generation''' + ''' + Test the perfect sweep generation. Test only specifics general tests + for frequency domain sweep synthesis in test_frequency_domain_sweep. + ''' sweep, group_delay = pfs.linear_perfect_sweep(n_samples, sampling_rate) @@ -365,6 +368,7 @@ def test_linear_perfect_sweep(n_samples, sampling_rate): assert type(sweep) is pf.Signal assert sweep.n_samples == n_samples assert sweep.sampling_rate == sampling_rate + assert sweep.fft_norm == 'rms' assert type(group_delay) is pf.FrequencyData # assert magnitude response @@ -382,3 +386,87 @@ def test_linear_perfect_sweep(n_samples, sampling_rate): # must be 0 - sweep is orthogonal to shifted versions of itself assert np.all( np.abs(auto_correlation) < 5 * np.finfo(sweep.time.dtype).eps) + + +@pytest.mark.parametrize('n_samples', [2**16, 2**16 + 1]) +@pytest.mark.parametrize('sampling_rate', [44100, 48000]) +def test_linear_sweep_freq(n_samples, sampling_rate): + ''' + Test the linear sweep generation. Test only specifics general tests + for frequency domain sweep synthesis in test_frequency_domain_sweep. + ''' + + sweep, group_delay = pfs.linear_sweep_freq( + n_samples, [200, 16e3], 5000, 5000, sampling_rate=sampling_rate) + + # basic checks + assert type(sweep) is pf.Signal + assert sweep.n_samples == n_samples + assert sweep.sampling_rate == sampling_rate + assert sweep.fft_norm == 'rms' + assert type(group_delay) is pf.FrequencyData + + # assert constant magnitude response within frequency range + idx = sweep.find_nearest_frequency([300, 15000]) + freq_db = pf.dsp.decibel(sweep).flatten()[idx[0]:idx[1]] + npt.assert_allclose(freq_db[idx[0]:idx[1]], + np.mean(freq_db[idx[0]:idx[1]]), + atol=.4) + + +@pytest.mark.parametrize('n_samples', [2**16, 2**16 + 1]) +@pytest.mark.parametrize('sampling_rate', [44100, 48000]) +def test_exponential_sweep_freq(n_samples, sampling_rate): + ''' + Test the exponential sweep generation. Test only specifics general tests + for frequency domain sweep synthesis in test_frequency_domain_sweep. + ''' + + sweep, group_delay = pfs.exponential_sweep_freq( + n_samples, [200, 16e3], 5000, 5000, sampling_rate=sampling_rate) + + # basic checks + assert type(sweep) is pf.Signal + assert sweep.n_samples == n_samples + assert sweep.sampling_rate == sampling_rate + assert sweep.fft_norm == 'rms' + assert type(group_delay) is pf.FrequencyData + + # assert 1/sqrt(f) magnitude response within frequency range + idx = sweep.find_nearest_frequency([300, 15000]) + freq_db = np.abs(sweep.freq_raw).flatten() + freq_db *= np.sqrt(sweep.frequencies) + freq_db = 20 * np.log10(freq_db[idx[0]:idx[1]]) + npt.assert_allclose(freq_db, np.mean(freq_db), atol=.4) + + +@pytest.mark.parametrize('n_samples', [2**16, 2**16 + 1]) +@pytest.mark.parametrize('sampling_rate', [44100, 48000]) +def test_magnitude_spectrum_weighted_sweep(n_samples, sampling_rate): + ''' + Test the magnitude weighted sweep. Test only specifics general tests + for frequency domain sweep synthesis in test_frequency_domain_sweep. + ''' + + # magnitude spectrum + magnitude = pf.dsp.filter.low_shelve( + pf.signals.impulse(n_samples, sampling_rate=sampling_rate), 500, 20, 2) + magnitude = pf.dsp.filter.butterworth(magnitude, 8, 50, 'highpass') + # sweep + sweep, group_delay = pf.signals.magnitude_spectrum_weighted_sweep( + n_samples, magnitude, 5000, 1000, sampling_rate=sampling_rate) + + # basic checks + assert type(sweep) is pf.Signal + assert sweep.n_samples == n_samples + assert sweep.sampling_rate == sampling_rate + assert sweep.fft_norm == 'rms' + assert type(group_delay) is pf.FrequencyData + + # assert magnitude response against reference within sweep frequency range + idx = sweep.find_nearest_frequency([50, sampling_rate / 2]) + freq_db = pf.dsp.decibel(sweep).flatten()[idx[0]:idx[1]] + freq_db -= np.mean(freq_db) + ref_db = pf.dsp.decibel(magnitude).flatten()[idx[0]:idx[1]] + ref_db -= np.mean(ref_db) + npt.assert_allclose(freq_db, ref_db, atol=.1) From edf1da8b222aaab7d13b2c54f8aef0f59b8ec671 Mon Sep 17 00:00:00 2001 From: Fabian Brinkmann Date: Sat, 28 Oct 2023 15:05:15 +0200 Subject: [PATCH 28/42] add more tests --- pyfar/signals/deterministic.py | 28 +++------ tests/test_signals.py | 105 ++++++++++++++++++++++++++++++--- 2 files changed, 105 insertions(+), 28 deletions(-) diff --git a/pyfar/signals/deterministic.py b/pyfar/signals/deterministic.py index 74a57a1cc..5cbb45a5f 100644 --- a/pyfar/signals/deterministic.py +++ b/pyfar/signals/deterministic.py @@ -724,33 +724,23 @@ def _frequency_domain_sweep( and ``sampling_rate/2``. """ - # check input ------------------------------------------------------------- - if not isinstance(sweep_type, (pyfar.Signal, str)): - raise TypeError("sweep_type must be type Signal or str.") + # check input (only checks for user input, no checks for calling the + # private function from the public functions) if isinstance(sweep_type, pyfar.Signal): magnitude = sweep_type sweep_type = 'signal' - if sweep_type not in ['linear', 'exponential', 'perfect_linear', 'signal']: - raise ValueError("sweep_type must be 'linear', 'exponential' or", - "'perfect_linear', when it is a str.") + if np.atleast_1d(frequency_range).size != 2: raise ValueError( "Frequency_range must be an array like with two elements.") - if frequency_range[1] > sampling_rate/2: + if frequency_range[0] < 0 or frequency_range[1] > sampling_rate/2: raise ValueError( - "Upper frequency limit is larger than half the sampling rate.") + ("Lower frequency limit must be at least 0 Hz and upper frequency " + "limit must be below half the sampling rate.")) if frequency_range[0] == 0 and sweep_type == "exponential": - Warning(( - "The exponential sweep has a 1/frequency magnitude spectrum. " - "The magnitude is set to 0 at 0 Hz to avoid division by zero.")) - if sweep_type == 'perfect_linear' and \ - (start_margin != 0 or stop_margin != 0 or double or - fade_in is not None or fade_out is not None or - frequency_range[0] != 0 or - frequency_range[1] != sampling_rate / 2): - # internal warning. Users will not call this function directly - # and can not cause this error. - raise ValueError(('Found conflicting parameters')) + Warning(UserWarning, ( + "The exponential sweep has a 1/sqrt(frequency) magnitude spectrum." + " The magnitude is set to 0 at 0 Hz to avoid division by zero.")) # initialize basic parameters --------------------------------------------- # double n_samples diff --git a/tests/test_signals.py b/tests/test_signals.py index 2a0d9d8f7..8302a0a51 100644 --- a/tests/test_signals.py +++ b/tests/test_signals.py @@ -358,8 +358,8 @@ def test_match_shape(): @pytest.mark.parametrize('sampling_rate', [44100, 48000]) def test_linear_perfect_sweep(n_samples, sampling_rate): ''' - Test the perfect sweep generation. Test only specifics general tests - for frequency domain sweep synthesis in test_frequency_domain_sweep. + Test the perfect sweep generation. General tests for frequency domain sweep + synthesis in test_frequency_domain_sweep. ''' sweep, group_delay = pfs.linear_perfect_sweep(n_samples, sampling_rate) @@ -392,8 +392,8 @@ def test_linear_perfect_sweep(n_samples, sampling_rate): @pytest.mark.parametrize('sampling_rate', [44100, 48000]) def test_linear_sweep_freq(n_samples, sampling_rate): ''' - Test the linear sweep generation. Test only specifics general tests - for frequency domain sweep synthesis in test_frequency_domain_sweep. + Test the linear sweep generation. General tests for frequency domain sweep + synthesis in test_frequency_domain_sweep. ''' sweep, group_delay = pfs.linear_sweep_freq( @@ -418,8 +418,8 @@ def test_linear_sweep_freq(n_samples, sampling_rate): @pytest.mark.parametrize('sampling_rate', [44100, 48000]) def test_exponential_sweep_freq(n_samples, sampling_rate): ''' - Test the exponential sweep generation. Test only specifics general tests - for frequency domain sweep synthesis in test_frequency_domain_sweep. + Test the exponential sweep generation. General tests for frequency domain + sweep synthesis in test_frequency_domain_sweep. ''' sweep, group_delay = pfs.exponential_sweep_freq( @@ -444,8 +444,8 @@ def test_exponential_sweep_freq(n_samples, sampling_rate): @pytest.mark.parametrize('sampling_rate', [44100, 48000]) def test_magnitude_spectrum_weighted_sweep(n_samples, sampling_rate): ''' - Test the magnitude weighted sweep. Test only specifics general tests - for frequency domain sweep synthesis in test_frequency_domain_sweep. + Test the magnitude weighted sweep. General tests for frequency domain sweep + synthesis in test_frequency_domain_sweep. ''' # magnitude spectrum @@ -453,7 +453,7 @@ def test_magnitude_spectrum_weighted_sweep(n_samples, sampling_rate): pf.signals.impulse(n_samples, sampling_rate=sampling_rate), 500, 20, 2) magnitude = pf.dsp.filter.butterworth(magnitude, 8, 50, 'highpass') # sweep - sweep, group_delay = pf.signals.magnitude_spectrum_weighted_sweep( + sweep, group_delay = pfs.magnitude_spectrum_weighted_sweep( n_samples, magnitude, 5000, 1000, sampling_rate=sampling_rate) # basic checks @@ -470,3 +470,90 @@ def test_magnitude_spectrum_weighted_sweep(n_samples, sampling_rate): ref_db = pf.dsp.decibel(magnitude).flatten()[idx[0]:idx[1]] ref_db -= np.mean(ref_db) npt.assert_allclose(freq_db, ref_db, atol=.1) + + +def test_magnitude_spectrum_weigthed_sweep_input_length(): + ''' + Test length of input magnitude spectrum. General tests for frequency domain + sweep synthesis in test_frequency_domain_sweep.''' + + # sorter signal will be zero padded + sweep_a, _ = pfs.magnitude_spectrum_weighted_sweep( + 2**10, pfs.impulse(2**9), 100, 100) + sweep_b, _ = pfs.magnitude_spectrum_weighted_sweep( + 2**10, pfs.impulse(2**10), 100, 100) + assert sweep_a == sweep_b + + with pytest.raises(ValueError, match='magnitude_spectrum has'): + pfs.magnitude_spectrum_weighted_sweep( + 2**9, pfs.impulse(2**10), 100, 100) + + +def test_frequency_domain_sweep(): + ''' + Test general frequency domain sweep synthesis. Specific tests for public + functions are contained in ``test_linear_perfect_sweep``, + ``test_linear_sweep_freq``, ``test_exponential_sweep_freq``, and + ``test_magnitude_spectrum_weighted_sweep``, + ''' + + # test frequency range - sweep with smaller range contains less energy + sweep_a, _ = pfs.linear_sweep_freq(2**14, [200, 20e3], 5000, 5000) + sweep_b, _ = pfs.linear_sweep_freq(2**14, [500, 15e3], 5000, 5000) + idx = sweep_a.find_nearest_frequency([200, 20e3]) + assert np.abs(sweep_a.freq[0, idx[0]]) > np.abs(sweep_b.freq[0, idx[0]]) + assert np.abs(sweep_a.freq[0, idx[1]]) > np.abs(sweep_b.freq[0, idx[1]]) + + # test single sided restriction of frequency range + sweep_a, _ = pfs.linear_sweep_freq(2**14, [0, 22050], 5000, 5000) + sweep_b, _ = pfs.linear_sweep_freq(2**14, [200, 20050], 5000, 5000) + assert np.abs(sweep_a.freq[0, idx[0]]) > np.abs(sweep_b.freq[0, idx[0]]) + + # test single sided restriction of frequency range + sweep_a, _ = pfs.linear_sweep_freq(2**14, [0, 22050], 5000, 5000) + sweep_b, _ = pfs.linear_sweep_freq(2**14, [0, 20000], 5000, 5000) + assert np.abs(sweep_a.freq[0, idx[1]]) > np.abs(sweep_b.freq[0, idx[1]]) + + # test single sided restriction of frequency range + sweep_a, _ = pfs.linear_sweep_freq(2**14, [0, 22050], 5000, 5000) + sweep_b, _ = pfs.linear_sweep_freq(2**14, [200, 22050], 5000, 5000) + assert np.abs(sweep_a.freq[0, idx[0]]) > np.abs(sweep_b.freq[0, idx[0]]) + + # test butterworth order - sweep with steeper filters contains less energy + sweep_a, _ = pfs.linear_sweep_freq( + 2**14, [400, 10e3], 5000, 5000, butterworth_order=4) + sweep_b, _ = pfs.linear_sweep_freq( + 2**14, [400, 10e3], 5000, 5000, butterworth_order=8) + assert np.abs(sweep_a.freq[0, idx[0]]) > np.abs(sweep_b.freq[0, idx[0]]) + assert np.abs(sweep_a.freq[0, idx[1]]) > np.abs(sweep_b.freq[0, idx[1]]) + + # test start and stop margin - shorter sweep contains less energy + sweep_a, _ = pfs.linear_sweep_freq(2**14, [200, 20e3], 5000, 5000) + sweep_b, _ = pfs.linear_sweep_freq(2**14, [200, 20e3], 7000, 7000) + assert pf.dsp.energy(sweep_a)[0] > pf.dsp.energy(sweep_b)[0] + + # test fade in and out - sweep with fades starts/ends with exact zeros + sweep_a, _ = pfs.linear_sweep_freq(2**14, [200, 20e3], 5000, 5000) + sweep_b, _ = pfs.linear_sweep_freq(2**14, [200, 20e3], 5000, 5000, 20, 20) + assert np.abs(sweep_a.time[0, 0]) > np.abs(sweep_b.time[0, 0]) + assert np.abs(sweep_a.time[0, -1]) > np.abs(sweep_b.time[0, -1]) + + # test fade in only - sweep with fades starts/ends with exact zeros + sweep_a, _ = pfs.linear_sweep_freq(2**14, [200, 20e3], 5000, 5000) + sweep_b, _ = pfs.linear_sweep_freq(2**14, [200, 20e3], 5000, 5000, 20, 0) + assert np.abs(sweep_a.time[0, 0]) > np.abs(sweep_b.time[0, 0]) + + # test fade out only - sweep with fades starts/ends with exact zeros + sweep_a, _ = pfs.linear_sweep_freq(2**14, [200, 20e3], 5000, 5000) + sweep_b, _ = pfs.linear_sweep_freq(2**14, [200, 20e3], 5000, 5000, 0, 20) + assert np.abs(sweep_a.time[0, -1]) > np.abs(sweep_b.time[0, -1]) + + # assertions for invalid frequency ranges + with pytest.raises(ValueError, match='with two elements'): + pfs.linear_sweep_freq(2**10, [0, 1, 2], 100, 100) + + with pytest.raises(ValueError, match='at least 0 Hz'): + pfs.linear_sweep_freq(2**10, [-1, 1], 100, 100) + + with pytest.raises(ValueError, match='below half the sampling rate'): + pfs.linear_sweep_freq(2**10, [0, 40000], 100, 100) From 4f2e289b43ebe7e744e2f2fc534d613945e6744f Mon Sep 17 00:00:00 2001 From: Fabian Brinkmann Date: Sat, 28 Oct 2023 20:03:50 +0200 Subject: [PATCH 29/42] finalize tests and functions --- pyfar/signals/deterministic.py | 30 +++++------------------------- tests/test_signals.py | 3 +++ 2 files changed, 8 insertions(+), 25 deletions(-) diff --git a/pyfar/signals/deterministic.py b/pyfar/signals/deterministic.py index 5cbb45a5f..b2cbca485 100644 --- a/pyfar/signals/deterministic.py +++ b/pyfar/signals/deterministic.py @@ -1,4 +1,5 @@ import numpy as np +import warnings import pyfar @@ -301,7 +302,6 @@ def linear_sweep_freq( sweep_type='linear', frequency_range=frequency_range, butterworth_order=butterworth_order, - double=False, start_margin=start_margin, stop_margin=stop_margin, fade_in=fade_in, @@ -483,7 +483,6 @@ def exponential_sweep_freq( sweep_type='exponential', frequency_range=frequency_range, butterworth_order=butterworth_order, - double=False, start_margin=start_margin, stop_margin=stop_margin, fade_in=fade_in, @@ -578,7 +577,6 @@ def magnitude_spectrum_weighted_sweep( sweep_type=magnitude_spectrum, frequency_range=[0, sampling_rate / 2], butterworth_order=0, - double=False, start_margin=start_margin, stop_margin=stop_margin, fade_in=fade_in, @@ -665,7 +663,6 @@ def linear_perfect_sweep(n_samples, sampling_rate=44100): sweep_type='perfect_linear', frequency_range=[0, sampling_rate / 2], butterworth_order=0, - double=False, start_margin=0, stop_margin=0, fade_in=None, @@ -676,7 +673,7 @@ def linear_perfect_sweep(n_samples, sampling_rate=44100): def _frequency_domain_sweep( - n_samples, sweep_type, frequency_range, butterworth_order, double, + n_samples, sweep_type, frequency_range, butterworth_order, start_margin, stop_margin, fade_in, fade_out, sampling_rate): """ Frequency domain sweep synthesis with arbitrary magnitude response. @@ -704,14 +701,7 @@ def _frequency_domain_sweep( energy in filters of relative constant bandwidth (e.g. octaves). ``'perfect_linear'`` Perfect linear sweep. Note that the parameters `start_margin`, - `stop_margin`, `frequency_range` and `double` are not required in - this case. - double : bool - Double `n_samples` during the sweep calculation. This parameter seems - not to be required in the current implementation. It was kept if cases - come up, where it could help. It is contained in the AKtools version - from which this function was ported. Must be ``False`` if `sweep_type` - is ``'perfect_linear'``. + `stop_margin`, and `frequency_range` are not required in this case. Returns ------- @@ -738,15 +728,11 @@ def _frequency_domain_sweep( ("Lower frequency limit must be at least 0 Hz and upper frequency " "limit must be below half the sampling rate.")) if frequency_range[0] == 0 and sweep_type == "exponential": - Warning(UserWarning, ( + warnings.warn(( "The exponential sweep has a 1/sqrt(frequency) magnitude spectrum." - " The magnitude is set to 0 at 0 Hz to avoid division by zero.")) + " The magnitude is set to 0 at 0 Hz to avoid a division by zero.")) # initialize basic parameters --------------------------------------------- - # double n_samples - if double and sweep_type != 'perfect_linear': - stop_margin += n_samples - n_samples *= 2 # spacing between frequency bins of FFT df = sampling_rate / n_samples @@ -839,12 +825,6 @@ def _frequency_domain_sweep( sweep = pyfar.Signal(sweep, sampling_rate, n_samples, 'freq', 'none') sweep.fft_norm = 'rms' - # cut to originally desired length - if double: - n_samples = n_samples // 2 - stop_margin -= n_samples - sweep.time = sweep.time[..., :n_samples] - # find windowing end points heuristically. The start and stop margin are # not reliable, because freq. domain synthesis causes pre/post-ringing if fade_in or fade_out: diff --git a/tests/test_signals.py b/tests/test_signals.py index 8302a0a51..ad0da5b10 100644 --- a/tests/test_signals.py +++ b/tests/test_signals.py @@ -557,3 +557,6 @@ def test_frequency_domain_sweep(): with pytest.raises(ValueError, match='below half the sampling rate'): pfs.linear_sweep_freq(2**10, [0, 40000], 100, 100) + + with pytest.warns(UserWarning, match='avoid a division by zero'): + pfs.exponential_sweep_freq(2*10, [0, 22050], 32, 32) From 89e4bbbe2de047542108f49c8ef224ecf6e6e239 Mon Sep 17 00:00:00 2001 From: Fabian Brinkmann Date: Sun, 29 Oct 2023 11:01:19 +0100 Subject: [PATCH 30/42] Fix docstring --- pyfar/signals/deterministic.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pyfar/signals/deterministic.py b/pyfar/signals/deterministic.py index b2cbca485..9a88546a2 100644 --- a/pyfar/signals/deterministic.py +++ b/pyfar/signals/deterministic.py @@ -159,7 +159,7 @@ def linear_sweep_time(n_samples, frequency_range, n_fade_out=90, amplitude=1, .. note:: The linear sweep can also be generated in the frequency domain - (see :py:func:`~exponential_sweep_freq`). Time domain synthesis + (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 of a slightly irregular temporal From 0f3693c27a9fd5c7ba99276821c1a9060c1d032d Mon Sep 17 00:00:00 2001 From: Fabian Brinkmann Date: Sat, 2 Dec 2023 09:23:35 +0100 Subject: [PATCH 31/42] Changes suggested by pingelit --- pyfar/signals/deterministic.py | 37 ++++++++++++++++++++-------------- 1 file changed, 22 insertions(+), 15 deletions(-) diff --git a/pyfar/signals/deterministic.py b/pyfar/signals/deterministic.py index 9a88546a2..b4f5cef9a 100644 --- a/pyfar/signals/deterministic.py +++ b/pyfar/signals/deterministic.py @@ -217,9 +217,9 @@ def linear_sweep_freq( n_samples, frequency_range, start_margin, stop_margin, fade_in=None, fade_out=None, butterworth_order=8, sampling_rate=44100): """ - Generate single channel sine sweep with linearlly increasing frequency. + Generate single channel sine sweep with linearly increasing frequency. - Sweep sweep synthesis according to [#]_ + Sine sweep sweep synthesis according to [#]_. .. note:: The linear sweep can also be generated in the time domain @@ -243,7 +243,7 @@ def linear_sweep_freq( stop_margin : int, float Time in samples, at which the sweep stops. This is relative to `n_samples`, e.g., a stop margin of 100 samples means that the sweep - ends at sample ``n_samples-10``. This is required, because the + ends at sample ``n_samples-100``. This is required, because the frequency domain sweep synthesis has post-ringing in the time domain. fade_in : int Duration of a squared sine fade-in in samples. The fade starts at the @@ -397,7 +397,7 @@ def exponential_sweep_freq( """ Generate single channel sine sweep with exponentially increasing frequency. - Sweep sweep synthesis according to [#]_ + Sweep sweep synthesis according to [#]_. .. note:: The exponential sweep can also be generated in the time domain @@ -423,7 +423,7 @@ def exponential_sweep_freq( stop_margin : int, float Time in samples, at which the sweep stops. This is relative to `n_samples`, e.g., a stop margin of 100 samples means that the sweep - ends at sample ``n_samples-10``. This is required, because the + ends at sample ``n_samples-100``. This is required, because the frequency domain sweep synthesis has post-ringing in the time domain. fade_in : int Duration of a squared sine fade-in in samples. The fade starts at the @@ -498,14 +498,14 @@ def magnitude_spectrum_weighted_sweep( """ Generate single channel sine sweep with arbitrary magnitude spectrum. - Sweep sweep synthesis according to [#]_ + Sine sweep sweep synthesis according to [#]_. .. note:: - The linear sweep can also be generated in the time domain - (:py:func:`~linear_sweep_time`). Frequency domain synthesis exhibits - smooth magnitude spectra in trade of a slightly irregular temporal - envelope. Time domain synthesis exhibits a constant temporal envelope - in trade of slight ripples in the magnitude response. + Frequency domain synthesis exhibits smooth magnitude spectra in trade + of a slightly irregular temporal envelope. Time domain synthesis + exhibits a constant temporal envelope in trade of slight ripples in the + magnitude response. But there is currently no method to design sine + sweeps with arbitrary magnitude response in the time domain. Parameters ---------- @@ -521,7 +521,7 @@ def magnitude_spectrum_weighted_sweep( stop_margin : int, float Time in samples, at which the sweep stops. This is relative to `n_samples`, e.g., a stop margin of 100 samples means that the sweep - ends at sample ``n_samples-10``. This is required, because the + ends at sample ``n_samples-100``. This is required, because the frequency domain sweep synthesis has post-ringing in the time domain. fade_in : int Duration of a squared sine fade-in in samples. The fade starts at the @@ -706,12 +706,19 @@ def _frequency_domain_sweep( Returns ------- sweep : Signal - The sweep signal. The Signal is in the time domain and has the ``none`` + The sweep signal. The Signal is in the time domain and has the ``rms`` FFT normalization (see :py:func:`~pyfar.dsp.fft.normalization`). The sweep parameters are written to `comment`. group_delay_sweep : FrequencyData The group delay of the sweep as a single sided spectrum between 0 Hz and ``sampling_rate/2``. + + References + ---------- + 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). """ # check input (only checks for user input, no checks for calling the @@ -790,11 +797,11 @@ def _frequency_domain_sweep( # group delay at Nyquist equals stopping time sweep_gd[-1] = (n_samples - stop_margin) / sampling_rate - # FORMULA (11, p.40 ) + # FORMULA (11, p.40, in Müller and Massarani 2001) sweep_power = np.sum(np.abs(sweep_abs**2)) C = (sweep_gd[-1] - sweep_gd[1]) / sweep_power - # FORMULA (10, p.40 ) + # FORMULA (10, p.40, in Müller and Massarani 2001) for k in range(tg_start, n_bins): # index 2 to nyq sweep_gd[k] = sweep_gd[k-1] + C * np.abs(sweep_abs[k])**2 From 8276679158dba3e1c7c66aec4fb6a2c49e4efe79 Mon Sep 17 00:00:00 2001 From: Fabian Brinkmann Date: Mon, 4 Dec 2023 08:53:18 +0100 Subject: [PATCH 32/42] Apply suggestion from pingelit --- pyfar/signals/deterministic.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/pyfar/signals/deterministic.py b/pyfar/signals/deterministic.py index b4f5cef9a..25e300e29 100644 --- a/pyfar/signals/deterministic.py +++ b/pyfar/signals/deterministic.py @@ -287,13 +287,13 @@ def linear_sweep_freq( Examples -------- - Linear sweep between 0 and 22050 Hz + Linear sweep between 50 and 22050 Hz .. plot:: >>> import pyfar as pf >>> sweep, _ = pf.signals.linear_sweep_freq( - ... 2**16, [0, 22050], 1000, 1000) + ... 2**16, [50, 22050], 1000, 1000) >>> pf.plot.time_freq(sweep) """ From 876ceb10c0da8a4defb76dacc623cf212fe2c401 Mon Sep 17 00:00:00 2001 From: Fabian Brinkmann Date: Tue, 12 Dec 2023 10:05:55 +0100 Subject: [PATCH 33/42] apply comment from ahms5 --- pyfar/signals/deterministic.py | 77 ++++++++++++++++++---------------- 1 file changed, 40 insertions(+), 37 deletions(-) diff --git a/pyfar/signals/deterministic.py b/pyfar/signals/deterministic.py index 25e300e29..282517fe7 100644 --- a/pyfar/signals/deterministic.py +++ b/pyfar/signals/deterministic.py @@ -157,13 +157,11 @@ def linear_sweep_time(n_samples, frequency_range, n_fade_out=90, amplitude=1, seconds, and the frequency limits :math:`f_\\mathrm{low}` and :math:`f_\\mathrm{high}`. - .. note:: - 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 of a slightly irregular temporal - envelope. + he 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 + of a slightly irregular temporal envelope. Parameters ---------- @@ -221,12 +219,17 @@ def linear_sweep_freq( Sine sweep sweep synthesis according to [#]_. + The linear sweep can also be generated in the time domain + (:py:func:`~linear_sweep_time`). Frequency domain synthesis exhibits + smooth magnitude spectra in trade of a slightly irregular temporal + envelope. Time domain synthesis exhibits a constant temporal envelope in + trade of slight ripples in the magnitude response. + .. note:: - The linear sweep can also be generated in the time domain - (:py:func:`~linear_sweep_time`). Frequency domain synthesis exhibits - smooth magnitude spectra in trade of a slightly irregular temporal - envelope. Time domain synthesis exhibits a constant temporal envelope - in trade of slight ripples in the magnitude response. + The envelope of the sweep time signal should be constant, apart from + slight overshoots at the beginning and end. If this is not the case, + try to increase `n_samples`, `start_margin`, `stop_margin`, `fade_in` + or `fade_out`. Parameters ---------- @@ -330,13 +333,11 @@ def exponential_sweep_time(n_samples, frequency_range, n_fade_out=90, seconds, and the frequency limits :math:`f_\\mathrm{low}` and :math:`f_\\mathrm{high}`. - .. note:: - The exponential sweep can also be generated in the frequency domain - (see :py:func:`~exponential_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 of a slightly irregular temporal - envelope. + The exponential sweep can also be generated in the frequency domain (see + see :py:func:`~exponential_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 of a slightly irregular temporal envelope. Parameters ---------- @@ -399,12 +400,17 @@ def exponential_sweep_freq( Sweep sweep synthesis according to [#]_. + 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 + temporal envelope. Time domain synthesis exhibits a constant temporal + envelope in trade of slight ripples in the magnitude response. + .. note:: - 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 - temporal envelope. Time domain synthesis exhibits a constant temporal - envelope in trade of slight ripples in the magnitude response. + The envelope of the sweep time signal should be constant, apart from + slight overshoots at the beginning and end. If this is not the case, + try to increase `n_samples`, `start_margin`, `stop_margin`, `fade_in` + or `fade_out`. Parameters ---------- @@ -452,13 +458,6 @@ def exponential_sweep_freq( group_delay_sweep : FrequencyData The group delay of the sweep in seconds as a single sided spectrum. - Notes - ----- - The envelope of the sweep time signal should be constant, apart from - slight overshoots at the beginning and end. If this is not the case, try - to increase `n_samples`, `start_margin`, `stop_margin`, `fade_in` or - `fade_out`. - References ---------- .. [#] S. Müller, P. Massarani. 'Transfer Function Measurement with Sweeps. @@ -498,14 +497,18 @@ def magnitude_spectrum_weighted_sweep( """ Generate single channel sine sweep with arbitrary magnitude spectrum. - Sine sweep sweep synthesis according to [#]_. + Sine sweep sweep synthesis according to [#]_. Frequency domain synthesis + exhibits smooth magnitude spectra in trade of a slightly irregular temporal + envelope. Time domain synthesis exhibits a constant temporal envelope in + trade of slight ripples in the magnitude response. But there is currently + no method to design sine sweeps with arbitrary magnitude response in the + time domain. .. note:: - Frequency domain synthesis exhibits smooth magnitude spectra in trade - of a slightly irregular temporal envelope. Time domain synthesis - exhibits a constant temporal envelope in trade of slight ripples in the - magnitude response. But there is currently no method to design sine - sweeps with arbitrary magnitude response in the time domain. + The envelope of the sweep time signal should be constant, apart from + slight overshoots at the beginning and end. If this is not the case, + try to increase `n_samples`, `start_margin`, `stop_margin`, `fade_in` + or `fade_out`, or provide a more smooth magnitude spectrum. Parameters ---------- From 6cc1af2139557bf69c71c0b4e7f208598174cd15 Mon Sep 17 00:00:00 2001 From: Fabian Brinkmann Date: Wed, 13 Dec 2023 17:16:23 +0100 Subject: [PATCH 34/42] Account for comments by ahms5 --- pyfar/signals/deterministic.py | 16 +--------------- 1 file changed, 1 insertion(+), 15 deletions(-) diff --git a/pyfar/signals/deterministic.py b/pyfar/signals/deterministic.py index 282517fe7..0b96bb28f 100644 --- a/pyfar/signals/deterministic.py +++ b/pyfar/signals/deterministic.py @@ -157,7 +157,7 @@ def linear_sweep_time(n_samples, frequency_range, n_fade_out=90, amplitude=1, seconds, and the frequency limits :math:`f_\\mathrm{low}` and :math:`f_\\mathrm{high}`. - he linear sweep can also be generated in the frequency domain (see + 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 @@ -275,13 +275,6 @@ def linear_sweep_freq( group_delay_sweep : FrequencyData The group delay of the sweep in seconds as a single sided spectrum. - Notes - ----- - The envelope of the sweep time signal should be constant, apart from - slight overshoots at the beginning and end. If this is not the case, try - to increase `n_samples`, `start_margin`, `stop_margin`, `fade_in` or - `fade_out`. - References ---------- .. [#] S. Müller, P. Massarani. 'Transfer Function Measurement with Sweeps. @@ -548,13 +541,6 @@ def magnitude_spectrum_weighted_sweep( group_delay_sweep : FrequencyData The group delay of the sweep in seconds as a single sided spectrum. - Notes - ----- - The envelope of the sweep time signal should be constant, apart from - slight overshoots at the beginning and end. If this is not the case, - provide a smother spectrum in `magnitude_spectrum` or try to increase - `n_samples`, `start_margin`, `stop_margin`, `fade_in` or `fade_out`. - References ---------- .. [#] S. Müller, P. Massarani. 'Transfer Function Measurement with Sweeps. From 1b90fc168639dc1d0074f376d00f15fafb67abda Mon Sep 17 00:00:00 2001 From: Fabian Brinkmann Date: Mon, 18 Dec 2023 10:42:58 +0100 Subject: [PATCH 35/42] update accroding to pingelit Co-authored-by: Pascal Palenda --- pyfar/signals/deterministic.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pyfar/signals/deterministic.py b/pyfar/signals/deterministic.py index 0b96bb28f..5d65f8b11 100644 --- a/pyfar/signals/deterministic.py +++ b/pyfar/signals/deterministic.py @@ -217,7 +217,7 @@ def linear_sweep_freq( """ Generate single channel sine sweep with linearly increasing frequency. - Sine sweep sweep synthesis according to [#]_. + Sine sweep synthesis according to [#]_. The linear sweep can also be generated in the time domain (:py:func:`~linear_sweep_time`). Frequency domain synthesis exhibits From c954e59e79ff65cfe604bfd269589543b723f2c7 Mon Sep 17 00:00:00 2001 From: Fabian Brinkmann Date: Thu, 21 Dec 2023 15:01:18 +0100 Subject: [PATCH 36/42] correct typos --- pyfar/signals/deterministic.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/pyfar/signals/deterministic.py b/pyfar/signals/deterministic.py index 5d65f8b11..4efef8b5a 100644 --- a/pyfar/signals/deterministic.py +++ b/pyfar/signals/deterministic.py @@ -391,7 +391,7 @@ def exponential_sweep_freq( """ Generate single channel sine sweep with exponentially increasing frequency. - Sweep sweep synthesis according to [#]_. + Sweep synthesis according to [#]_. The exponential sweep can also be generated in the time domain (:py:func:`~exponential_sweep_time`). Frequency domain synthesis @@ -490,7 +490,7 @@ def magnitude_spectrum_weighted_sweep( """ Generate single channel sine sweep with arbitrary magnitude spectrum. - Sine sweep sweep synthesis according to [#]_. Frequency domain synthesis + Sine sweep synthesis according to [#]_. Frequency domain synthesis exhibits smooth magnitude spectra in trade of a slightly irregular temporal envelope. Time domain synthesis exhibits a constant temporal envelope in trade of slight ripples in the magnitude response. But there is currently From e1d3bff1e3a695efa1531122981d88844e2065f2 Mon Sep 17 00:00:00 2001 From: Fabian Brinkmann Date: Thu, 4 Apr 2024 18:39:21 +0200 Subject: [PATCH 37/42] account for review by --- pyfar/signals/deterministic.py | 19 ++++++++++++------- 1 file changed, 12 insertions(+), 7 deletions(-) diff --git a/pyfar/signals/deterministic.py b/pyfar/signals/deterministic.py index 4efef8b5a..9af141efb 100644 --- a/pyfar/signals/deterministic.py +++ b/pyfar/signals/deterministic.py @@ -145,7 +145,8 @@ def impulse(n_samples, delay=0, amplitude=1, sampling_rate=44100): def linear_sweep_time(n_samples, frequency_range, n_fade_out=90, amplitude=1, sampling_rate=44100): - """Generate single channel sine sweep with linearly increasing frequency. + """ + Generate sine sweep with linearly increasing frequency in the time domain. Time domain sweep generation according to [#]_: @@ -215,7 +216,8 @@ def linear_sweep_freq( n_samples, frequency_range, start_margin, stop_margin, fade_in=None, fade_out=None, butterworth_order=8, sampling_rate=44100): """ - Generate single channel sine sweep with linearly increasing frequency. + Generate sine sweep with linearly increasing frequency in the frequency + domain. Sine sweep synthesis according to [#]_. @@ -310,7 +312,8 @@ def linear_sweep_freq( def exponential_sweep_time(n_samples, frequency_range, n_fade_out=90, amplitude=1, sweep_rate=None, sampling_rate=44100): """ - Generate single channel sine sweep with exponentially increasing frequency. + Generate sine sweep with exponentially increasing frequency in the time + domain. Time domain sweep generation according to [#]_: @@ -327,7 +330,7 @@ def exponential_sweep_time(n_samples, frequency_range, n_fade_out=90, :math:`f_\\mathrm{high}`. The exponential sweep can also be generated in the frequency domain (see - see :py:func:`~exponential_sweep.freq`). Time domain synthesis exhibits a + see :py:func:`~exponential_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 of a slightly irregular temporal envelope. @@ -389,7 +392,8 @@ def exponential_sweep_freq( n_samples, frequency_range, start_margin, stop_margin, fade_in=None, fade_out=None, butterworth_order=8, sampling_rate=44100): """ - Generate single channel sine sweep with exponentially increasing frequency. + Generate sine sweep with exponentially increasing frequency in the + frequency domain. Sweep synthesis according to [#]_. @@ -488,7 +492,8 @@ def magnitude_spectrum_weighted_sweep( n_samples, magnitude_spectrum, start_margin, stop_margin, fade_in=None, fade_out=None, sampling_rate=44100): """ - Generate single channel sine sweep with arbitrary magnitude spectrum. + Generate sine sweep with arbitrary magnitude spectrum in the frequency + domain. Sine sweep synthesis according to [#]_. Frequency domain synthesis exhibits smooth magnitude spectra in trade of a slightly irregular temporal @@ -577,7 +582,7 @@ def magnitude_spectrum_weighted_sweep( def linear_perfect_sweep(n_samples, sampling_rate=44100): """ - Generate a perfect linear sweep. + Generate a perfect linear sweep in the frequency domain. The perfect sweep is generated according to [#]_ and is used for adaptive system identification. It is orthogonal to delayed versions of itself and From 9916ec8ad71f57af4bc21f97bcb345b2c1e8ef5c Mon Sep 17 00:00:00 2001 From: Fabian Brinkmann Date: Fri, 5 Apr 2024 10:27:15 +0200 Subject: [PATCH 38/42] update according to review --- pyfar/signals/deterministic.py | 62 ++++++++++++++++++---------------- tests/test_signals.py | 8 ++--- 2 files changed, 36 insertions(+), 34 deletions(-) diff --git a/pyfar/signals/deterministic.py b/pyfar/signals/deterministic.py index c5138cda3..ee5120c58 100644 --- a/pyfar/signals/deterministic.py +++ b/pyfar/signals/deterministic.py @@ -161,7 +161,7 @@ def linear_sweep_time(n_samples, frequency_range, n_fade_out=90, amplitude=1, 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 + Frequency domain synthesis exhibits smooth magnitude spectra in trade of a slightly irregular temporal envelope. Parameters @@ -240,7 +240,8 @@ def linear_sweep_freq( frequency_range : array like Frequency range of the sweep given by the lower and upper cut-off frequency in Hz. The restriction of the frequency range is realized - by applying a Butterworth band-pass with the specified frequencies. + by applying a Butterworth high-pass if ``frequency_range[0] > 0`` + and/or by a low-pass if ``frequency_range[1] < sampling_rate / 2``. start_margin : int, float The time in samples, at which the sweep starts. The start margin is required because the frequency domain sweep synthesis has pre-ringing @@ -253,18 +254,16 @@ def linear_sweep_freq( 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 + maximum of the sweep time signal. The default ``None`` does not apply a fade-in. 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 + maximum of the sweep time signal. The default ``None`` does not apply a fade-out. 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``. + frequency range (see above). The default is ``8``. sampling_rate : int The sampling rate in Hz. The default is ``44100``. @@ -272,10 +271,12 @@ def linear_sweep_freq( ------- sweep : Signal The sweep signal. The Signal is in the time domain, has a maximum - absolute amplitude of 1 and the ``rms`` FFT normalization + absolute amplitude of 1 and the ``none`` FFT normalization (see :py:func:`~pyfar.dsp.fft.normalization`). group_delay_sweep : FrequencyData - The group delay of the sweep in seconds as a single sided spectrum. + The analytical group delay of the sweep in seconds as a single sided + spectrum. This can be used to compute the times at which distortion + products appear. References ---------- @@ -332,7 +333,7 @@ def exponential_sweep_time(n_samples, frequency_range, n_fade_out=90, The exponential sweep can also be generated in the frequency domain (see see :py:func:`~exponential_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 + response. Frequency domain synthesis exhibits smooth magnitude spectra in trade of a slightly irregular temporal envelope. Parameters @@ -417,7 +418,8 @@ def exponential_sweep_freq( frequency_range : array like Frequency range of the sweep given by the lower and upper cut-off frequency in Hz. The restriction of the frequency range is realized - by applying a Butterworth band-pass with the specified frequencies. + by applying a Butterworth high-pass if ``frequency_range[0] > 0`` + and/or by a low-pass if ``frequency_range[1] < sampling_rate / 2``. Note that the exponential sweep can not start at 0 Hz, because its magnitude is defined by 1/frequency. start_margin : int, float @@ -432,18 +434,16 @@ def exponential_sweep_freq( 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 + maximum of the sweep time signal. The default ``None`` does not apply a fade-in. 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 + maximum of the sweep time signal. The default ``None`` does not apply a fade-out. 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``. + frequency range (see above). The default is ``8``. sampling_rate : int The sampling rate in Hz. The default is ``44100``. @@ -451,10 +451,12 @@ def exponential_sweep_freq( ------- sweep : Signal The sweep signal. The Signal is in the time domain, has a maximum - absolute amplitude of 1 and the ``rms`` FFT normalization + absolute amplitude of 1 and the ``none`` FFT normalization (see :py:func:`~pyfar.dsp.fft.normalization`). group_delay_sweep : FrequencyData - The group delay of the sweep in seconds as a single sided spectrum. + The analytical group delay of the sweep in seconds as a single sided + spectrum. This can be used to compute the times at which distortion + products appear. References ---------- @@ -496,10 +498,7 @@ def magnitude_spectrum_weighted_sweep( Generate sine sweep with arbitrary magnitude spectrum in the frequency domain. - Sine sweep synthesis according to [#]_. Frequency domain synthesis - exhibits smooth magnitude spectra in trade of a slightly irregular temporal - envelope. Time domain synthesis exhibits a constant temporal envelope in - trade of slight ripples in the magnitude response. But there is currently + Frequency domain sine sweep synthesis according to [#]_. There is currently no method to design sine sweeps with arbitrary magnitude response in the time domain. @@ -528,12 +527,12 @@ def magnitude_spectrum_weighted_sweep( 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 + maximum of the sweep time signal. The default ``None`` does not apply a fade-in. 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 + maximum of the sweep time signal. The default ``None`` does not apply a fade-out. sampling_rate : int The sampling rate in Hz. The default is ``44100``. @@ -542,10 +541,12 @@ def magnitude_spectrum_weighted_sweep( ------- sweep : Signal The sweep signal. The Signal is in the time domain, has a maximum - absolute amplitude of 1 and the ``rms`` FFT normalization + absolute amplitude of 1 and the ``none`` FFT normalization (see :py:func:`~pyfar.dsp.fft.normalization`). group_delay_sweep : FrequencyData - The group delay of the sweep in seconds as a single sided spectrum. + The analytical group delay of the sweep in seconds as a single sided + spectrum. This can be used to compute the times at which distortion + products appear. References ---------- @@ -600,10 +601,12 @@ def linear_perfect_sweep(n_samples, sampling_rate=44100): ------- sweep : Signal The sweep signal. The Signal is in the time domain, has a maximum - absolute amplitude of 1 and the ``rms`` FFT normalization + absolute amplitude of 1 and the ``none`` FFT normalization (see :py:func:`~pyfar.dsp.fft.normalization`). group_delay_sweep : FrequencyData - The group delay of the sweep in seconds as a single sided spectrum. + The analytical group delay of the sweep in seconds as a single sided + spectrum. This can be used to compute the times at which distortion + products appear. References ---------- @@ -701,7 +704,7 @@ def _frequency_domain_sweep( Returns ------- sweep : Signal - The sweep signal. The Signal is in the time domain and has the ``rms`` + The sweep signal. The Signal is in the time domain and has the ``none`` FFT normalization (see :py:func:`~pyfar.dsp.fft.normalization`). The sweep parameters are written to `comment`. group_delay_sweep : FrequencyData @@ -825,7 +828,6 @@ def _frequency_domain_sweep( # put sweep in pyfar.Signal an transform to time domain sweep = pyfar.Signal(sweep, sampling_rate, n_samples, 'freq', 'none') - sweep.fft_norm = 'rms' # find windowing end points heuristically. The start and stop margin are # not reliable, because freq. domain synthesis causes pre/post-ringing diff --git a/tests/test_signals.py b/tests/test_signals.py index ad0da5b10..0a12a4235 100644 --- a/tests/test_signals.py +++ b/tests/test_signals.py @@ -368,7 +368,7 @@ def test_linear_perfect_sweep(n_samples, sampling_rate): assert type(sweep) is pf.Signal assert sweep.n_samples == n_samples assert sweep.sampling_rate == sampling_rate - assert sweep.fft_norm == 'rms' + assert sweep.fft_norm == 'none' assert type(group_delay) is pf.FrequencyData # assert magnitude response @@ -403,7 +403,7 @@ def test_linear_sweep_freq(n_samples, sampling_rate): assert type(sweep) is pf.Signal assert sweep.n_samples == n_samples assert sweep.sampling_rate == sampling_rate - assert sweep.fft_norm == 'rms' + assert sweep.fft_norm == 'none' assert type(group_delay) is pf.FrequencyData # assert constant magnitude response within frequency range @@ -429,7 +429,7 @@ def test_exponential_sweep_freq(n_samples, sampling_rate): assert type(sweep) is pf.Signal assert sweep.n_samples == n_samples assert sweep.sampling_rate == sampling_rate - assert sweep.fft_norm == 'rms' + assert sweep.fft_norm == 'none' assert type(group_delay) is pf.FrequencyData # assert 1/sqrt(f) magnitude response within frequency range @@ -460,7 +460,7 @@ def test_magnitude_spectrum_weighted_sweep(n_samples, sampling_rate): assert type(sweep) is pf.Signal assert sweep.n_samples == n_samples assert sweep.sampling_rate == sampling_rate - assert sweep.fft_norm == 'rms' + assert sweep.fft_norm == 'none' assert type(group_delay) is pf.FrequencyData # assert magnitude response against reference within sweep frequency range From 68669cad6c91d59ada7ba1a7b40823d83970476b Mon Sep 17 00:00:00 2001 From: Fabian Brinkmann Date: Fri, 5 Apr 2024 10:33:45 +0200 Subject: [PATCH 39/42] cleaner return structure --- pyfar/signals/deterministic.py | 16 ++++------------ 1 file changed, 4 insertions(+), 12 deletions(-) diff --git a/pyfar/signals/deterministic.py b/pyfar/signals/deterministic.py index ee5120c58..91abf1006 100644 --- a/pyfar/signals/deterministic.py +++ b/pyfar/signals/deterministic.py @@ -296,7 +296,7 @@ def linear_sweep_freq( >>> pf.plot.time_freq(sweep) """ - signal, group_delay = _frequency_domain_sweep( + return _frequency_domain_sweep( n_samples=n_samples, sweep_type='linear', frequency_range=frequency_range, @@ -307,8 +307,6 @@ def linear_sweep_freq( fade_out=fade_out, sampling_rate=sampling_rate) - return signal, group_delay - def exponential_sweep_time(n_samples, frequency_range, n_fade_out=90, amplitude=1, sweep_rate=None, sampling_rate=44100): @@ -477,7 +475,7 @@ def exponential_sweep_freq( >>> pf.plot.time_freq(sweep) """ - signal, group_delay = _frequency_domain_sweep( + return _frequency_domain_sweep( n_samples=n_samples, sweep_type='exponential', frequency_range=frequency_range, @@ -488,8 +486,6 @@ def exponential_sweep_freq( fade_out=fade_out, sampling_rate=sampling_rate) - return signal, group_delay - def magnitude_spectrum_weighted_sweep( n_samples, magnitude_spectrum, start_margin, stop_margin, @@ -568,7 +564,7 @@ def magnitude_spectrum_weighted_sweep( >>> pf.plot.time_freq(sweep) """ - signal, group_delay = _frequency_domain_sweep( + return _frequency_domain_sweep( n_samples=n_samples, sweep_type=magnitude_spectrum, frequency_range=[0, sampling_rate / 2], @@ -579,8 +575,6 @@ def magnitude_spectrum_weighted_sweep( fade_out=fade_out, sampling_rate=sampling_rate) - return signal, group_delay - def linear_perfect_sweep(n_samples, sampling_rate=44100): """ @@ -656,7 +650,7 @@ def linear_perfect_sweep(n_samples, sampling_rate=44100): >>> plt.gca().set_ylabel('auto correlation') """ - signal, group_delay = _frequency_domain_sweep( + return _frequency_domain_sweep( n_samples=n_samples, sweep_type='perfect_linear', frequency_range=[0, sampling_rate / 2], @@ -667,8 +661,6 @@ def linear_perfect_sweep(n_samples, sampling_rate=44100): fade_out=None, sampling_rate=sampling_rate) - return signal, group_delay - def _frequency_domain_sweep( n_samples, sweep_type, frequency_range, butterworth_order, From be1b624ed38f9e31b18fca5dba5572b107207864 Mon Sep 17 00:00:00 2001 From: Fabian Brinkmann Date: Fri, 5 Apr 2024 14:10:03 +0200 Subject: [PATCH 40/42] make group delay optional return parameter and improve docstrings --- pyfar/signals/deterministic.py | 97 +++++++++++++++++++++------------- tests/test_signals.py | 83 +++++++++++++++++++---------- 2 files changed, 115 insertions(+), 65 deletions(-) diff --git a/pyfar/signals/deterministic.py b/pyfar/signals/deterministic.py index 91abf1006..d114731e4 100644 --- a/pyfar/signals/deterministic.py +++ b/pyfar/signals/deterministic.py @@ -214,7 +214,8 @@ def linear_sweep_time(n_samples, frequency_range, n_fade_out=90, amplitude=1, def linear_sweep_freq( n_samples, frequency_range, start_margin, stop_margin, fade_in=None, - fade_out=None, butterworth_order=8, sampling_rate=44100): + fade_out=None, butterworth_order=8, sampling_rate=44100, + group_delay=False): """ Generate sine sweep with linearly increasing frequency in the frequency domain. @@ -251,21 +252,25 @@ def linear_sweep_freq( `n_samples`, e.g., a stop margin of 100 samples means that the sweep ends at sample ``n_samples-100``. This is required, because the frequency domain sweep synthesis has post-ringing in the time domain. - fade_in : int + fade_in : int, optional 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 ``None`` does not apply a fade-in. - fade_out : int + fade_out : int, optional 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 ``None`` does not apply a fade-out. - butterworth_order : int, None + butterworth_order : int, optional The order of the Butterworth filters that are applied to limit the frequency range (see above). The default is ``8``. - sampling_rate : int + sampling_rate : int, optional The sampling rate in Hz. The default is ``44100``. + 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``. Returns ------- @@ -275,8 +280,7 @@ def linear_sweep_freq( (see :py:func:`~pyfar.dsp.fft.normalization`). group_delay_sweep : FrequencyData The analytical group delay of the sweep in seconds as a single sided - spectrum. This can be used to compute the times at which distortion - products appear. + spectrum. Only returned if `group_delay` is ``True``. References ---------- @@ -291,7 +295,7 @@ def linear_sweep_freq( .. plot:: >>> import pyfar as pf - >>> sweep, _ = pf.signals.linear_sweep_freq( + >>> sweep = pf.signals.linear_sweep_freq( ... 2**16, [50, 22050], 1000, 1000) >>> pf.plot.time_freq(sweep) """ @@ -305,7 +309,8 @@ def linear_sweep_freq( stop_margin=stop_margin, fade_in=fade_in, fade_out=fade_out, - sampling_rate=sampling_rate) + sampling_rate=sampling_rate, + group_delay=group_delay) def exponential_sweep_time(n_samples, frequency_range, n_fade_out=90, @@ -390,7 +395,8 @@ def exponential_sweep_time(n_samples, frequency_range, n_fade_out=90, def exponential_sweep_freq( n_samples, frequency_range, start_margin, stop_margin, fade_in=None, - fade_out=None, butterworth_order=8, sampling_rate=44100): + fade_out=None, butterworth_order=8, sampling_rate=44100, + group_delay=False): """ Generate sine sweep with exponentially increasing frequency in the frequency domain. @@ -429,21 +435,25 @@ def exponential_sweep_freq( `n_samples`, e.g., a stop margin of 100 samples means that the sweep ends at sample ``n_samples-100``. This is required, because the frequency domain sweep synthesis has post-ringing in the time domain. - fade_in : int + fade_in : int, optional 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 ``None`` does not apply a fade-in. - fade_out : int + fade_out : int, optional 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 ``None`` does not apply a fade-out. - butterworth_order : int, None + butterworth_order : int, optional The order of the Butterworth filters that are applied to limit the frequency range (see above). The default is ``8``. - sampling_rate : int + sampling_rate : int, optional The sampling rate in Hz. The default is ``44100``. + 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``. Returns ------- @@ -453,8 +463,7 @@ def exponential_sweep_freq( (see :py:func:`~pyfar.dsp.fft.normalization`). group_delay_sweep : FrequencyData The analytical group delay of the sweep in seconds as a single sided - spectrum. This can be used to compute the times at which distortion - products appear. + spectrum. Only returned if `group_delay` is ``True``. References ---------- @@ -470,7 +479,7 @@ def exponential_sweep_freq( .. plot:: >>> import pyfar as pf - >>> sweep, _ = pf.signals.exponential_sweep_freq( + >>> sweep = pf.signals.exponential_sweep_freq( ... 2**16, [50, 22050], 5000, 1000) >>> pf.plot.time_freq(sweep) """ @@ -484,12 +493,14 @@ def exponential_sweep_freq( stop_margin=stop_margin, fade_in=fade_in, fade_out=fade_out, - sampling_rate=sampling_rate) + sampling_rate=sampling_rate, + group_delay=group_delay) def magnitude_spectrum_weighted_sweep( n_samples, magnitude_spectrum, start_margin, stop_margin, - fade_in=None, fade_out=None, sampling_rate=44100): + fade_in=None, fade_out=None, sampling_rate=44100, + group_delay=False): """ Generate sine sweep with arbitrary magnitude spectrum in the frequency domain. @@ -520,18 +531,22 @@ def magnitude_spectrum_weighted_sweep( `n_samples`, e.g., a stop margin of 100 samples means that the sweep ends at sample ``n_samples-100``. This is required, because the frequency domain sweep synthesis has post-ringing in the time domain. - fade_in : int + fade_in : int, optional 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 ``None`` does not apply a fade-in. - fade_out : int + fade_out : int, optional 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 ``None`` does not apply a fade-out. - sampling_rate : int + sampling_rate : int, optional The sampling rate in Hz. The default is ``44100``. + 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``. Returns ------- @@ -541,8 +556,7 @@ def magnitude_spectrum_weighted_sweep( (see :py:func:`~pyfar.dsp.fft.normalization`). group_delay_sweep : FrequencyData The analytical group delay of the sweep in seconds as a single sided - spectrum. This can be used to compute the times at which distortion - products appear. + spectrum. Only returned if `group_delay` is ``True``. References ---------- @@ -559,7 +573,7 @@ def magnitude_spectrum_weighted_sweep( >>> 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( + >>> sweep = pf.signals.magnitude_spectrum_weighted_sweep( ... 2**16, magnitude, 5000, 1000) >>> pf.plot.time_freq(sweep) """ @@ -573,10 +587,11 @@ def magnitude_spectrum_weighted_sweep( stop_margin=stop_margin, fade_in=fade_in, fade_out=fade_out, - sampling_rate=sampling_rate) + sampling_rate=sampling_rate, + group_delay=group_delay) -def linear_perfect_sweep(n_samples, sampling_rate=44100): +def linear_perfect_sweep(n_samples, sampling_rate=44100, group_delay=False): """ Generate a perfect linear sweep in the frequency domain. @@ -588,8 +603,12 @@ def linear_perfect_sweep(n_samples, sampling_rate=44100): ---------- n_samples : int The length of the sweep in samples. - sampling_rate : int + sampling_rate : int, optional The sampling rate in Hz. The default is ``44100``. + 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``. Returns ------- @@ -599,8 +618,7 @@ def linear_perfect_sweep(n_samples, sampling_rate=44100): (see :py:func:`~pyfar.dsp.fft.normalization`). group_delay_sweep : FrequencyData The analytical group delay of the sweep in seconds as a single sided - spectrum. This can be used to compute the times at which distortion - products appear. + spectrum. Only returned if `group_delay` is ``True``. References ---------- @@ -616,7 +634,7 @@ def linear_perfect_sweep(n_samples, sampling_rate=44100): .. plot:: >>> import pyfar as pf - >>> sweep, _ = pf.signals.linear_perfect_sweep(2**8) + >>> sweep = pf.signals.linear_perfect_sweep(2**8) >>> sweep_shifted = pf.dsp.time_shift(sweep, 2**7) >>> ax = pf.plot.time_freq(sweep_shifted, unit='samples') >>> ax[0].axvline(2**7, ls='--', c='r', label='start/end') @@ -631,7 +649,7 @@ def linear_perfect_sweep(n_samples, sampling_rate=44100): >>> import numpy as np >>> import matplotlib.pyplot as plt >>> - >>> sweep, _ = pf.signals.linear_perfect_sweep(2**8) + >>> sweep = pf.signals.linear_perfect_sweep(2**8) >>> >>> # compute auto-correlation >>> auto_correlation = np.empty(2**8) @@ -659,12 +677,14 @@ def linear_perfect_sweep(n_samples, sampling_rate=44100): stop_margin=0, fade_in=None, fade_out=None, - sampling_rate=sampling_rate) + sampling_rate=sampling_rate, + group_delay=group_delay) def _frequency_domain_sweep( n_samples, sweep_type, frequency_range, butterworth_order, - start_margin, stop_margin, fade_in, fade_out, sampling_rate): + start_margin, stop_margin, fade_in, fade_out, sampling_rate, + group_delay): """ Frequency domain sweep synthesis with arbitrary magnitude response. @@ -700,8 +720,8 @@ def _frequency_domain_sweep( FFT normalization (see :py:func:`~pyfar.dsp.fft.normalization`). The sweep parameters are written to `comment`. group_delay_sweep : FrequencyData - The group delay of the sweep as a single sided spectrum between 0 Hz - and ``sampling_rate/2``. + The analytical group delay of the sweep in seconds as a single sided + spectrum. Only returned if `group_delay` is ``True``. References ---------- @@ -852,7 +872,10 @@ def _frequency_domain_sweep( # (to avoid clipping if written to fixed point wav file) sweep = pyfar.dsp.normalize(sweep) * (1 - 2**-15) - return sweep, sweep_gd + if group_delay: + return sweep, sweep_gd + else: + return sweep def _time_domain_sweep(n_samples, frequency_range, n_fade_out, amplitude, diff --git a/tests/test_signals.py b/tests/test_signals.py index 0a12a4235..41ec47744 100644 --- a/tests/test_signals.py +++ b/tests/test_signals.py @@ -362,19 +362,25 @@ def test_linear_perfect_sweep(n_samples, sampling_rate): synthesis in test_frequency_domain_sweep. ''' - sweep, group_delay = pfs.linear_perfect_sweep(n_samples, sampling_rate) + sweep = pfs.linear_perfect_sweep(n_samples, sampling_rate) # basic checks assert type(sweep) is pf.Signal assert sweep.n_samples == n_samples assert sweep.sampling_rate == sampling_rate assert sweep.fft_norm == 'none' - assert type(group_delay) is pf.FrequencyData # assert magnitude response npt.assert_almost_equal(np.abs(sweep.freq_raw), np.mean(np.abs(sweep.freq_raw))) + # test returning the group delay + sweep_2, group_delay = pfs.linear_perfect_sweep( + n_samples, sampling_rate, group_delay=True) + + assert type(group_delay) is pf.FrequencyData + assert sweep_2 == sweep + # compute normalized auto-correlation auto_correlation = np.empty(n_samples-1) for shift in range(1, n_samples): @@ -396,7 +402,7 @@ def test_linear_sweep_freq(n_samples, sampling_rate): synthesis in test_frequency_domain_sweep. ''' - sweep, group_delay = pfs.linear_sweep_freq( + sweep = pfs.linear_sweep_freq( n_samples, [200, 16e3], 5000, 5000, sampling_rate=sampling_rate) # basic checks @@ -404,7 +410,6 @@ def test_linear_sweep_freq(n_samples, sampling_rate): assert sweep.n_samples == n_samples assert sweep.sampling_rate == sampling_rate assert sweep.fft_norm == 'none' - assert type(group_delay) is pf.FrequencyData # assert constant magnitude response within frequency range idx = sweep.find_nearest_frequency([300, 15000]) @@ -413,6 +418,14 @@ def test_linear_sweep_freq(n_samples, sampling_rate): np.mean(freq_db[idx[0]:idx[1]]), atol=.4) + # test returning the group delay + sweep_2, group_delay = pfs.linear_sweep_freq( + n_samples, [200, 16e3], 5000, 5000, sampling_rate=sampling_rate, + group_delay=True) + + assert type(group_delay) is pf.FrequencyData + assert sweep == sweep_2 + @pytest.mark.parametrize('n_samples', [2**16, 2**16 + 1]) @pytest.mark.parametrize('sampling_rate', [44100, 48000]) @@ -422,7 +435,7 @@ def test_exponential_sweep_freq(n_samples, sampling_rate): sweep synthesis in test_frequency_domain_sweep. ''' - sweep, group_delay = pfs.exponential_sweep_freq( + sweep = pfs.exponential_sweep_freq( n_samples, [200, 16e3], 5000, 5000, sampling_rate=sampling_rate) # basic checks @@ -430,7 +443,6 @@ def test_exponential_sweep_freq(n_samples, sampling_rate): assert sweep.n_samples == n_samples assert sweep.sampling_rate == sampling_rate assert sweep.fft_norm == 'none' - assert type(group_delay) is pf.FrequencyData # assert 1/sqrt(f) magnitude response within frequency range idx = sweep.find_nearest_frequency([300, 15000]) @@ -439,6 +451,14 @@ def test_exponential_sweep_freq(n_samples, sampling_rate): freq_db = 20 * np.log10(freq_db[idx[0]:idx[1]]) npt.assert_allclose(freq_db, np.mean(freq_db), atol=.4) + # test returning the group delay + sweep_2, group_delay = pfs.exponential_sweep_freq( + n_samples, [200, 16e3], 5000, 5000, sampling_rate=sampling_rate, + group_delay=True) + + assert type(group_delay) is pf.FrequencyData + assert sweep == sweep_2 + @pytest.mark.parametrize('n_samples', [2**16, 2**16 + 1]) @pytest.mark.parametrize('sampling_rate', [44100, 48000]) @@ -453,7 +473,7 @@ def test_magnitude_spectrum_weighted_sweep(n_samples, sampling_rate): pf.signals.impulse(n_samples, sampling_rate=sampling_rate), 500, 20, 2) magnitude = pf.dsp.filter.butterworth(magnitude, 8, 50, 'highpass') # sweep - sweep, group_delay = pfs.magnitude_spectrum_weighted_sweep( + sweep = pfs.magnitude_spectrum_weighted_sweep( n_samples, magnitude, 5000, 1000, sampling_rate=sampling_rate) # basic checks @@ -461,7 +481,6 @@ def test_magnitude_spectrum_weighted_sweep(n_samples, sampling_rate): assert sweep.n_samples == n_samples assert sweep.sampling_rate == sampling_rate assert sweep.fft_norm == 'none' - assert type(group_delay) is pf.FrequencyData # assert magnitude response against reference within sweep frequency range idx = sweep.find_nearest_frequency([50, sampling_rate / 2]) @@ -471,6 +490,14 @@ def test_magnitude_spectrum_weighted_sweep(n_samples, sampling_rate): ref_db -= np.mean(ref_db) npt.assert_allclose(freq_db, ref_db, atol=.1) + # test returning the group delay + sweep_2, group_delay = pfs.magnitude_spectrum_weighted_sweep( + n_samples, magnitude, 5000, 1000, sampling_rate=sampling_rate, + group_delay=True) + + assert type(group_delay) is pf.FrequencyData + assert sweep == sweep_2 + def test_magnitude_spectrum_weigthed_sweep_input_length(): ''' @@ -478,9 +505,9 @@ def test_magnitude_spectrum_weigthed_sweep_input_length(): sweep synthesis in test_frequency_domain_sweep.''' # sorter signal will be zero padded - sweep_a, _ = pfs.magnitude_spectrum_weighted_sweep( + sweep_a = pfs.magnitude_spectrum_weighted_sweep( 2**10, pfs.impulse(2**9), 100, 100) - sweep_b, _ = pfs.magnitude_spectrum_weighted_sweep( + sweep_b = pfs.magnitude_spectrum_weighted_sweep( 2**10, pfs.impulse(2**10), 100, 100) assert sweep_a == sweep_b @@ -498,54 +525,54 @@ def test_frequency_domain_sweep(): ''' # test frequency range - sweep with smaller range contains less energy - sweep_a, _ = pfs.linear_sweep_freq(2**14, [200, 20e3], 5000, 5000) - sweep_b, _ = pfs.linear_sweep_freq(2**14, [500, 15e3], 5000, 5000) + sweep_a = pfs.linear_sweep_freq(2**14, [200, 20e3], 5000, 5000) + sweep_b = pfs.linear_sweep_freq(2**14, [500, 15e3], 5000, 5000) idx = sweep_a.find_nearest_frequency([200, 20e3]) assert np.abs(sweep_a.freq[0, idx[0]]) > np.abs(sweep_b.freq[0, idx[0]]) assert np.abs(sweep_a.freq[0, idx[1]]) > np.abs(sweep_b.freq[0, idx[1]]) # test single sided restriction of frequency range - sweep_a, _ = pfs.linear_sweep_freq(2**14, [0, 22050], 5000, 5000) - sweep_b, _ = pfs.linear_sweep_freq(2**14, [200, 20050], 5000, 5000) + sweep_a = pfs.linear_sweep_freq(2**14, [0, 22050], 5000, 5000) + sweep_b = pfs.linear_sweep_freq(2**14, [200, 20050], 5000, 5000) assert np.abs(sweep_a.freq[0, idx[0]]) > np.abs(sweep_b.freq[0, idx[0]]) # test single sided restriction of frequency range - sweep_a, _ = pfs.linear_sweep_freq(2**14, [0, 22050], 5000, 5000) - sweep_b, _ = pfs.linear_sweep_freq(2**14, [0, 20000], 5000, 5000) + sweep_a = pfs.linear_sweep_freq(2**14, [0, 22050], 5000, 5000) + sweep_b = pfs.linear_sweep_freq(2**14, [0, 20000], 5000, 5000) assert np.abs(sweep_a.freq[0, idx[1]]) > np.abs(sweep_b.freq[0, idx[1]]) # test single sided restriction of frequency range - sweep_a, _ = pfs.linear_sweep_freq(2**14, [0, 22050], 5000, 5000) - sweep_b, _ = pfs.linear_sweep_freq(2**14, [200, 22050], 5000, 5000) + sweep_a = pfs.linear_sweep_freq(2**14, [0, 22050], 5000, 5000) + sweep_b = pfs.linear_sweep_freq(2**14, [200, 22050], 5000, 5000) assert np.abs(sweep_a.freq[0, idx[0]]) > np.abs(sweep_b.freq[0, idx[0]]) # test butterworth order - sweep with steeper filters contains less energy - sweep_a, _ = pfs.linear_sweep_freq( + sweep_a = pfs.linear_sweep_freq( 2**14, [400, 10e3], 5000, 5000, butterworth_order=4) - sweep_b, _ = pfs.linear_sweep_freq( + sweep_b = pfs.linear_sweep_freq( 2**14, [400, 10e3], 5000, 5000, butterworth_order=8) assert np.abs(sweep_a.freq[0, idx[0]]) > np.abs(sweep_b.freq[0, idx[0]]) assert np.abs(sweep_a.freq[0, idx[1]]) > np.abs(sweep_b.freq[0, idx[1]]) # test start and stop margin - shorter sweep contains less energy - sweep_a, _ = pfs.linear_sweep_freq(2**14, [200, 20e3], 5000, 5000) - sweep_b, _ = pfs.linear_sweep_freq(2**14, [200, 20e3], 7000, 7000) + sweep_a = pfs.linear_sweep_freq(2**14, [200, 20e3], 5000, 5000) + sweep_b = pfs.linear_sweep_freq(2**14, [200, 20e3], 7000, 7000) assert pf.dsp.energy(sweep_a)[0] > pf.dsp.energy(sweep_b)[0] # test fade in and out - sweep with fades starts/ends with exact zeros - sweep_a, _ = pfs.linear_sweep_freq(2**14, [200, 20e3], 5000, 5000) - sweep_b, _ = pfs.linear_sweep_freq(2**14, [200, 20e3], 5000, 5000, 20, 20) + sweep_a = pfs.linear_sweep_freq(2**14, [200, 20e3], 5000, 5000) + sweep_b = pfs.linear_sweep_freq(2**14, [200, 20e3], 5000, 5000, 20, 20) assert np.abs(sweep_a.time[0, 0]) > np.abs(sweep_b.time[0, 0]) assert np.abs(sweep_a.time[0, -1]) > np.abs(sweep_b.time[0, -1]) # test fade in only - sweep with fades starts/ends with exact zeros - sweep_a, _ = pfs.linear_sweep_freq(2**14, [200, 20e3], 5000, 5000) - sweep_b, _ = pfs.linear_sweep_freq(2**14, [200, 20e3], 5000, 5000, 20, 0) + sweep_a = pfs.linear_sweep_freq(2**14, [200, 20e3], 5000, 5000) + sweep_b = pfs.linear_sweep_freq(2**14, [200, 20e3], 5000, 5000, 20, 0) assert np.abs(sweep_a.time[0, 0]) > np.abs(sweep_b.time[0, 0]) # test fade out only - sweep with fades starts/ends with exact zeros - sweep_a, _ = pfs.linear_sweep_freq(2**14, [200, 20e3], 5000, 5000) - sweep_b, _ = pfs.linear_sweep_freq(2**14, [200, 20e3], 5000, 5000, 0, 20) + sweep_a = pfs.linear_sweep_freq(2**14, [200, 20e3], 5000, 5000) + sweep_b = pfs.linear_sweep_freq(2**14, [200, 20e3], 5000, 5000, 0, 20) assert np.abs(sweep_a.time[0, -1]) > np.abs(sweep_b.time[0, -1]) # assertions for invalid frequency ranges From 13f12b76b5730ab7bbf5b0ca6643b5dfae665849 Mon Sep 17 00:00:00 2001 From: Fabian Brinkmann Date: Wed, 24 Apr 2024 16:42:34 +0200 Subject: [PATCH 41/42] update according to review from mberz --- pyfar/signals/deterministic.py | 101 ++++++++++++++++----------------- tests/test_signals.py | 12 ++-- 2 files changed, 56 insertions(+), 57 deletions(-) diff --git a/pyfar/signals/deterministic.py b/pyfar/signals/deterministic.py index d114731e4..930c36f74 100644 --- a/pyfar/signals/deterministic.py +++ b/pyfar/signals/deterministic.py @@ -160,8 +160,8 @@ def linear_sweep_time(n_samples, frequency_range, n_fade_out=90, amplitude=1, 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 in trade + temporal envelope at the cost of slight ripples in the magnitude response. + Frequency domain synthesis exhibits smooth magnitude spectra at the cost of a slightly irregular temporal envelope. Parameters @@ -214,8 +214,8 @@ def linear_sweep_time(n_samples, frequency_range, n_fade_out=90, amplitude=1, def linear_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): + fade_out=None, bandpass_order=8, sampling_rate=44100, + return_group_delay=False): """ Generate sine sweep with linearly increasing frequency in the frequency domain. @@ -224,7 +224,7 @@ def linear_sweep_freq( The linear sweep can also be generated in the time domain (:py:func:`~linear_sweep_time`). Frequency domain synthesis exhibits - smooth magnitude spectra in trade of a slightly irregular temporal + smooth magnitude spectra at the cost of a slightly irregular temporal envelope. Time domain synthesis exhibits a constant temporal envelope in trade of slight ripples in the magnitude response. @@ -262,12 +262,12 @@ def linear_sweep_freq( last sample of the sweep that is closer than 60 dB to the absolute maximum of the sweep time signal. The default ``None`` does not apply a fade-out. - butterworth_order : int, optional + bandpass_order : int, optional The order of the Butterworth filters that are applied to limit the frequency range (see above). The default is ``8``. sampling_rate : int, optional The sampling rate in Hz. The default is ``44100``. - group_delay : boolean, optional + return_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``. @@ -280,13 +280,13 @@ def linear_sweep_freq( (see :py:func:`~pyfar.dsp.fft.normalization`). group_delay_sweep : FrequencyData The analytical group delay of the sweep in seconds as a single sided - spectrum. Only returned if `group_delay` is ``True``. + spectrum. Only returned if `return_group_delay` is ``True``. References ---------- .. [#] 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. + Corrections.' J. Audio Eng. Soc. 2001, 49 (6), 443-471. Examples -------- @@ -304,13 +304,13 @@ def linear_sweep_freq( n_samples=n_samples, sweep_type='linear', frequency_range=frequency_range, - butterworth_order=butterworth_order, + bandpass_order=bandpass_order, start_margin=start_margin, stop_margin=stop_margin, fade_in=fade_in, fade_out=fade_out, sampling_rate=sampling_rate, - group_delay=group_delay) + return_group_delay=return_group_delay) def exponential_sweep_time(n_samples, frequency_range, n_fade_out=90, @@ -335,9 +335,9 @@ def exponential_sweep_time(n_samples, frequency_range, n_fade_out=90, The exponential sweep can also be generated in the frequency domain (see see :py:func:`~exponential_sweep_freq`). Time domain synthesis exhibits a - constant temporal envelope in trade of slight ripples in the magnitude + constant temporal envelope at the cost of slight ripples in the magnitude response. Frequency domain synthesis exhibits smooth magnitude spectra - in trade of a slightly irregular temporal envelope. + at the cost of a slightly irregular temporal envelope. Parameters ---------- @@ -395,8 +395,8 @@ def exponential_sweep_time(n_samples, frequency_range, n_fade_out=90, 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): + fade_out=None, bandpass_order=8, sampling_rate=44100, + return_group_delay=False): """ Generate sine sweep with exponentially increasing frequency in the frequency domain. @@ -405,9 +405,9 @@ def exponential_sweep_freq( 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 + exhibits smooth magnitude spectra at the cost of a slightly irregular temporal envelope. Time domain synthesis exhibits a constant temporal - envelope in trade of slight ripples in the magnitude response. + envelope at the cost of slight ripples in the magnitude response. .. note:: The envelope of the sweep time signal should be constant, apart from @@ -445,12 +445,12 @@ def exponential_sweep_freq( last sample of the sweep that is closer than 60 dB to the absolute maximum of the sweep time signal. The default ``None`` does not apply a fade-out. - butterworth_order : int, optional + bandpass_order : int, optional The order of the Butterworth filters that are applied to limit the frequency range (see above). The default is ``8``. sampling_rate : int, optional The sampling rate in Hz. The default is ``44100``. - group_delay : boolean, optional + return_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``. @@ -463,13 +463,13 @@ def exponential_sweep_freq( (see :py:func:`~pyfar.dsp.fft.normalization`). group_delay_sweep : FrequencyData The analytical group delay of the sweep in seconds as a single sided - spectrum. Only returned if `group_delay` is ``True``. + spectrum. Only returned if `return_group_delay` is ``True``. References ---------- .. [#] 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. + Corrections.' J. Audio Eng. Soc. 2001, 49 (6), 443-471. Examples -------- @@ -488,19 +488,19 @@ def exponential_sweep_freq( n_samples=n_samples, sweep_type='exponential', frequency_range=frequency_range, - butterworth_order=butterworth_order, + bandpass_order=bandpass_order, start_margin=start_margin, stop_margin=stop_margin, fade_in=fade_in, fade_out=fade_out, sampling_rate=sampling_rate, - group_delay=group_delay) + return_group_delay=return_group_delay) def magnitude_spectrum_weighted_sweep( n_samples, magnitude_spectrum, start_margin, stop_margin, fade_in=None, fade_out=None, sampling_rate=44100, - group_delay=False): + return_group_delay=False): """ Generate sine sweep with arbitrary magnitude spectrum in the frequency domain. @@ -543,7 +543,7 @@ def magnitude_spectrum_weighted_sweep( a fade-out. sampling_rate : int, optional The sampling rate in Hz. The default is ``44100``. - group_delay : boolean, optional + return_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``. @@ -556,13 +556,13 @@ def magnitude_spectrum_weighted_sweep( (see :py:func:`~pyfar.dsp.fft.normalization`). group_delay_sweep : FrequencyData The analytical group delay of the sweep in seconds as a single sided - spectrum. Only returned if `group_delay` is ``True``. + spectrum. Only returned if `return_group_delay` is ``True``. References ---------- .. [#] 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. + Corrections.' J. Audio Eng. Soc. 2001, 49 (6), 443-471. Examples -------- @@ -582,22 +582,22 @@ def magnitude_spectrum_weighted_sweep( n_samples=n_samples, sweep_type=magnitude_spectrum, frequency_range=[0, sampling_rate / 2], - butterworth_order=0, + bandpass_order=0, start_margin=start_margin, stop_margin=stop_margin, fade_in=fade_in, fade_out=fade_out, sampling_rate=sampling_rate, - group_delay=group_delay) + return_group_delay=return_group_delay) -def linear_perfect_sweep(n_samples, sampling_rate=44100, group_delay=False): +def linear_perfect_sweep( + n_samples, sampling_rate=44100, return_group_delay=False): """ Generate a perfect linear sweep in the frequency domain. The perfect sweep is generated according to [#]_ and is used for adaptive - system identification. It is orthogonal to delayed versions of itself and - can be played back in a loop. + system identification. It is orthogonal to delayed versions of itself. Parameters ---------- @@ -605,7 +605,7 @@ def linear_perfect_sweep(n_samples, sampling_rate=44100, group_delay=False): The length of the sweep in samples. sampling_rate : int, optional The sampling rate in Hz. The default is ``44100``. - group_delay : boolean, optional + return_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``. @@ -618,7 +618,7 @@ def linear_perfect_sweep(n_samples, sampling_rate=44100, group_delay=False): (see :py:func:`~pyfar.dsp.fft.normalization`). group_delay_sweep : FrequencyData The analytical group delay of the sweep in seconds as a single sided - spectrum. Only returned if `group_delay` is ``True``. + spectrum. Only returned if `return_group_delay` is ``True``. References ---------- @@ -672,21 +672,21 @@ def linear_perfect_sweep(n_samples, sampling_rate=44100, group_delay=False): n_samples=n_samples, sweep_type='perfect_linear', frequency_range=[0, sampling_rate / 2], - butterworth_order=0, + bandpass_order=0, start_margin=0, stop_margin=0, fade_in=None, fade_out=None, sampling_rate=sampling_rate, - group_delay=group_delay) + return_group_delay=return_group_delay) def _frequency_domain_sweep( - n_samples, sweep_type, frequency_range, butterworth_order, + n_samples, sweep_type, frequency_range, bandpass_order, start_margin, stop_margin, fade_in, fade_out, sampling_rate, - group_delay): + return_group_delay): """ - Frequency domain sweep synthesis with arbitrary magnitude response. + Frequency domain sweep synthesis with arbitrary magnitude response [#]_. Non-unique parameters are documented in ``linear_sweep_freq``, ``exponential_sweep_freq``, ``linear_perfect_sweep``, and @@ -721,14 +721,14 @@ def _frequency_domain_sweep( sweep parameters are written to `comment`. group_delay_sweep : FrequencyData The analytical group delay of the sweep in seconds as a single sided - spectrum. Only returned if `group_delay` is ``True``. + spectrum. Only returned if `return_group_delay` is ``True``. References ---------- - 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). + .. [#] 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). """ # check input (only checks for user input, no checks for calling the @@ -781,15 +781,15 @@ def _frequency_domain_sweep( if frequency_range[0] > 0 and frequency_range[1] < sampling_rate / 2: band_limit = pyfar.dsp.filter.butterworth( pyfar.signals.impulse(n_samples, sampling_rate=sampling_rate), - butterworth_order, frequency_range, 'bandpass') + bandpass_order, frequency_range, 'bandpass') elif frequency_range[0] > 0: band_limit = pyfar.dsp.filter.butterworth( pyfar.signals.impulse(n_samples, sampling_rate=sampling_rate), - butterworth_order, frequency_range[0], 'highpass') + bandpass_order, frequency_range[0], 'highpass') elif frequency_range[1] < sampling_rate / 2: band_limit = pyfar.dsp.filter.butterworth( pyfar.signals.impulse(n_samples, sampling_rate=sampling_rate), - butterworth_order, frequency_range[1], 'lowpass') + bandpass_order, frequency_range[1], 'lowpass') else: band_limit = pyfar.Signal(np.ones_like(sweep_abs), sampling_rate, n_samples, 'freq') @@ -868,11 +868,10 @@ def _frequency_domain_sweep( if fade is not None: sweep = pyfar.dsp.time_window(sweep, fade, shape=shape) - # 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) + # normalize to time domain amplitude + sweep = pyfar.dsp.normalize(sweep) - if group_delay: + if return_group_delay: return sweep, sweep_gd else: return sweep diff --git a/tests/test_signals.py b/tests/test_signals.py index 41ec47744..a9c6827eb 100644 --- a/tests/test_signals.py +++ b/tests/test_signals.py @@ -376,7 +376,7 @@ def test_linear_perfect_sweep(n_samples, sampling_rate): # test returning the group delay sweep_2, group_delay = pfs.linear_perfect_sweep( - n_samples, sampling_rate, group_delay=True) + n_samples, sampling_rate, return_group_delay=True) assert type(group_delay) is pf.FrequencyData assert sweep_2 == sweep @@ -421,7 +421,7 @@ def test_linear_sweep_freq(n_samples, sampling_rate): # test returning the group delay sweep_2, group_delay = pfs.linear_sweep_freq( n_samples, [200, 16e3], 5000, 5000, sampling_rate=sampling_rate, - group_delay=True) + return_group_delay=True) assert type(group_delay) is pf.FrequencyData assert sweep == sweep_2 @@ -454,7 +454,7 @@ def test_exponential_sweep_freq(n_samples, sampling_rate): # test returning the group delay sweep_2, group_delay = pfs.exponential_sweep_freq( n_samples, [200, 16e3], 5000, 5000, sampling_rate=sampling_rate, - group_delay=True) + return_group_delay=True) assert type(group_delay) is pf.FrequencyData assert sweep == sweep_2 @@ -493,7 +493,7 @@ def test_magnitude_spectrum_weighted_sweep(n_samples, sampling_rate): # test returning the group delay sweep_2, group_delay = pfs.magnitude_spectrum_weighted_sweep( n_samples, magnitude, 5000, 1000, sampling_rate=sampling_rate, - group_delay=True) + return_group_delay=True) assert type(group_delay) is pf.FrequencyData assert sweep == sweep_2 @@ -548,9 +548,9 @@ def test_frequency_domain_sweep(): # test butterworth order - sweep with steeper filters contains less energy sweep_a = pfs.linear_sweep_freq( - 2**14, [400, 10e3], 5000, 5000, butterworth_order=4) + 2**14, [400, 10e3], 5000, 5000, bandpass_order=4) sweep_b = pfs.linear_sweep_freq( - 2**14, [400, 10e3], 5000, 5000, butterworth_order=8) + 2**14, [400, 10e3], 5000, 5000, bandpass_order=8) assert np.abs(sweep_a.freq[0, idx[0]]) > np.abs(sweep_b.freq[0, idx[0]]) assert np.abs(sweep_a.freq[0, idx[1]]) > np.abs(sweep_b.freq[0, idx[1]]) From fb4d817443dc1e59604176c3cb14a9a4f9368f23 Mon Sep 17 00:00:00 2001 From: Fabian Brinkmann Date: Thu, 25 Apr 2024 10:24:14 +0200 Subject: [PATCH 42/42] rename fade parameters --- pyfar/signals/deterministic.py | 80 +++++++++++++++++----------------- 1 file changed, 40 insertions(+), 40 deletions(-) diff --git a/pyfar/signals/deterministic.py b/pyfar/signals/deterministic.py index 930c36f74..25c2ac675 100644 --- a/pyfar/signals/deterministic.py +++ b/pyfar/signals/deterministic.py @@ -213,8 +213,8 @@ def linear_sweep_time(n_samples, frequency_range, n_fade_out=90, amplitude=1, def linear_sweep_freq( - 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, n_fade_in=0, + n_fade_out=0, bandpass_order=8, sampling_rate=44100, return_group_delay=False): """ Generate sine sweep with linearly increasing frequency in the frequency @@ -231,8 +231,8 @@ def linear_sweep_freq( .. note:: The envelope of the sweep time signal should be constant, apart from slight overshoots at the beginning and end. If this is not the case, - try to increase `n_samples`, `start_margin`, `stop_margin`, `fade_in` - or `fade_out`. + try to increase `n_samples`, `start_margin`, `stop_margin`, `n_fade_in` + or `n_fade_out`. Parameters ---------- @@ -252,15 +252,15 @@ def linear_sweep_freq( `n_samples`, e.g., a stop margin of 100 samples means that the sweep ends at sample ``n_samples-100``. This is required, because the frequency domain sweep synthesis has post-ringing in the time domain. - fade_in : int, optional + n_fade_in : int, optional 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 ``None`` does not apply + maximum of the sweep time signal. The default ``0`` does not apply a fade-in. - fade_out : int, optional + n_fade_out : int, optional 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 ``None`` does not apply + maximum of the sweep time signal. The default ``0`` does not apply a fade-out. bandpass_order : int, optional The order of the Butterworth filters that are applied to limit the @@ -307,8 +307,8 @@ def linear_sweep_freq( bandpass_order=bandpass_order, start_margin=start_margin, stop_margin=stop_margin, - fade_in=fade_in, - fade_out=fade_out, + n_fade_in=n_fade_in, + n_fade_out=n_fade_out, sampling_rate=sampling_rate, return_group_delay=return_group_delay) @@ -394,8 +394,8 @@ def exponential_sweep_time(n_samples, frequency_range, n_fade_out=90, def exponential_sweep_freq( - 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, n_fade_in=0, + n_fade_out=0, bandpass_order=8, sampling_rate=44100, return_group_delay=False): """ Generate sine sweep with exponentially increasing frequency in the @@ -412,8 +412,8 @@ def exponential_sweep_freq( .. note:: The envelope of the sweep time signal should be constant, apart from slight overshoots at the beginning and end. If this is not the case, - try to increase `n_samples`, `start_margin`, `stop_margin`, `fade_in` - or `fade_out`. + try to increase `n_samples`, `start_margin`, `stop_margin`, `n_fade_in` + or `n_fade_out`. Parameters ---------- @@ -435,15 +435,15 @@ def exponential_sweep_freq( `n_samples`, e.g., a stop margin of 100 samples means that the sweep ends at sample ``n_samples-100``. This is required, because the frequency domain sweep synthesis has post-ringing in the time domain. - fade_in : int, optional + n_fade_in : int, optional 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 ``None`` does not apply + maximum of the sweep time signal. The default ``0`` does not apply a fade-in. - fade_out : int, optional + n_fade_out : int, optional 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 ``None`` does not apply + maximum of the sweep time signal. The default ``0`` does not apply a fade-out. bandpass_order : int, optional The order of the Butterworth filters that are applied to limit the @@ -491,15 +491,15 @@ def exponential_sweep_freq( bandpass_order=bandpass_order, start_margin=start_margin, stop_margin=stop_margin, - fade_in=fade_in, - fade_out=fade_out, + n_fade_in=n_fade_in, + n_fade_out=n_fade_out, sampling_rate=sampling_rate, return_group_delay=return_group_delay) def magnitude_spectrum_weighted_sweep( n_samples, magnitude_spectrum, start_margin, stop_margin, - fade_in=None, fade_out=None, sampling_rate=44100, + n_fade_in=0, n_fade_out=0, sampling_rate=44100, return_group_delay=False): """ Generate sine sweep with arbitrary magnitude spectrum in the frequency @@ -512,8 +512,8 @@ def magnitude_spectrum_weighted_sweep( .. note:: The envelope of the sweep time signal should be constant, apart from slight overshoots at the beginning and end. If this is not the case, - try to increase `n_samples`, `start_margin`, `stop_margin`, `fade_in` - or `fade_out`, or provide a more smooth magnitude spectrum. + try to increase `n_samples`, `start_margin`, `stop_margin`, `n_fade_in` + or `n_fade_out`, or provide a more smooth magnitude spectrum. Parameters ---------- @@ -531,15 +531,15 @@ def magnitude_spectrum_weighted_sweep( `n_samples`, e.g., a stop margin of 100 samples means that the sweep ends at sample ``n_samples-100``. This is required, because the frequency domain sweep synthesis has post-ringing in the time domain. - fade_in : int, optional + n_fade_in : int, optional 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 ``None`` does not apply + maximum of the sweep time signal. The default ``0`` does not apply a fade-in. - fade_out : int, optional + n_fade_out : int, optional 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 ``None`` does not apply + maximum of the sweep time signal. The default ``0`` does not apply a fade-out. sampling_rate : int, optional The sampling rate in Hz. The default is ``44100``. @@ -585,8 +585,8 @@ def magnitude_spectrum_weighted_sweep( bandpass_order=0, start_margin=start_margin, stop_margin=stop_margin, - fade_in=fade_in, - fade_out=fade_out, + n_fade_in=n_fade_in, + n_fade_out=n_fade_out, sampling_rate=sampling_rate, return_group_delay=return_group_delay) @@ -675,15 +675,15 @@ def linear_perfect_sweep( bandpass_order=0, start_margin=0, stop_margin=0, - fade_in=None, - fade_out=None, + n_fade_in=0, + n_fade_out=0, sampling_rate=sampling_rate, return_group_delay=return_group_delay) def _frequency_domain_sweep( n_samples, sweep_type, frequency_range, bandpass_order, - start_margin, stop_margin, fade_in, fade_out, sampling_rate, + start_margin, stop_margin, n_fade_in, n_fade_out, sampling_rate, return_group_delay): """ Frequency domain sweep synthesis with arbitrary magnitude response [#]_. @@ -843,7 +843,7 @@ def _frequency_domain_sweep( # find windowing end points heuristically. The start and stop margin are # not reliable, because freq. domain synthesis causes pre/post-ringing - if fade_in or fade_out: + if n_fade_in or n_fade_out: threshold = 10**(-60/20) * np.max(np.abs(sweep.time)) fade_start = np.argmax(np.abs(sweep.time) > threshold) fade_end = n_samples - \ @@ -851,15 +851,15 @@ def _frequency_domain_sweep( print([fade_start, fade_end]) # generate windows for fade in and/or out - if fade_in and fade_out: - fade = [fade_start, fade_start + fade_in, - fade_end - fade_out, fade_end] + if n_fade_in and n_fade_out: + fade = [fade_start, fade_start + n_fade_in, + fade_end - n_fade_out, fade_end] shape = 'symmetric' - elif fade_in: - fade = [fade_start, fade_start + fade_in] + elif n_fade_in: + fade = [fade_start, fade_start + n_fade_in] shape = 'left' - elif fade_out: - fade = [fade_end - fade_out, fade_end] + elif n_fade_out: + fade = [fade_end - n_fade_out, fade_end] shape = 'right' else: fade = None