From 96d56d95bb204000c627523b841be314ba0e5f8b Mon Sep 17 00:00:00 2001 From: Dietrich Brunn Date: Fri, 15 Nov 2024 00:18:50 +0100 Subject: [PATCH] ENH: signal.chirp(): Add functionality of complex-valued waveforms. * Added boolean parameter `complex` * Simplified `signal.chirp()` code * Added unit tests (including switching over to `xp_assert_allclose`) * Improved docstr Co-authored-by: Anil Gurses --- scipy/signal/_waveforms.py | 192 ++++++++++++++------------- scipy/signal/tests/test_waveforms.py | 27 ++++ 2 files changed, 129 insertions(+), 90 deletions(-) diff --git a/scipy/signal/_waveforms.py b/scipy/signal/_waveforms.py index a1a6888763ca..a6be46cfd386 100644 --- a/scipy/signal/_waveforms.py +++ b/scipy/signal/_waveforms.py @@ -257,8 +257,9 @@ def gausspulse(t, fc=1000, bw=0.5, bwr=-6, tpr=-60, retquad=False, return yI, yQ, yenv -def chirp(t, f0, t1, f1, method='linear', phi=0, vertex_zero=True): - """Frequency-swept cosine generator. +def chirp(t, f0, t1, f1, method='linear', phi=0, vertex_zero=True, *, + complex=False): + r"""Frequency-swept cosine generator. In the following, 'Hz' should be interpreted as 'cycles per unit'; there is no requirement here that the unit is one second. The @@ -284,14 +285,21 @@ def chirp(t, f0, t1, f1, method='linear', phi=0, vertex_zero=True): This parameter is only used when `method` is 'quadratic'. It determines whether the vertex of the parabola that is the graph of the frequency is at t=0 or t=t1. + complex : bool, optional + This parameter creates a complex-valued analytic signal instead of a + real-valued signal. It allows the use of complex baseband (in communications + domain). Default is False. + + .. versionadded:: 1.15.0 Returns ------- y : ndarray - A numpy array containing the signal evaluated at `t` with the - requested time-varying frequency. More precisely, the function - returns ``cos(phase + (pi/180)*phi)`` where `phase` is the integral - (from 0 to `t`) of ``2*pi*f(t)``. ``f(t)`` is defined below. + A numpy array containing the signal evaluated at `t` with the requested + time-varying frequency. More precisely, the function returns + ``exp(1j*phase + 1j*(pi/180)*phi) if complex else cos(phase + (pi/180)*phi)`` + where `phase` is the integral (from 0 to `t`) of ``2*pi*f(t)``. + The instantaneous frequency ``f(t)`` is defined below. See Also -------- @@ -299,120 +307,124 @@ def chirp(t, f0, t1, f1, method='linear', phi=0, vertex_zero=True): Notes ----- - There are four options for the `method`. The following formulas give - the instantaneous frequency (in Hz) of the signal generated by - `chirp()`. For convenience, the shorter names shown below may also be - used. - - linear, lin, li: - - ``f(t) = f0 + (f1 - f0) * t / t1`` - - quadratic, quad, q: + There are four possible options for the parameter `method`, which have a (long) + standard form and some allowed abbreviations. The formulas for the instantaneous + frequency :math:`f(t)` of the generated signal are as follows: - The graph of the frequency f(t) is a parabola through (0, f0) and - (t1, f1). By default, the vertex of the parabola is at (0, f0). - If `vertex_zero` is False, then the vertex is at (t1, f1). The - formula is: + 1. Parameter `method` in ``('linear', 'lin', 'li')``: - if vertex_zero is True: + .. math:: + f(t) = f_0 + \beta\, t \quad\text{with}\quad + \beta = \frac{f_1 - f_0}{t_1} - ``f(t) = f0 + (f1 - f0) * t**2 / t1**2`` + Frequency :math:`f(t)` varies linearly over time with a constant rate + :math:`\beta`. - else: + 2. Parameter `method` in ``('quadratic', 'quad', 'q')``: - ``f(t) = f1 - (f1 - f0) * (t1 - t)**2 / t1**2`` + .. math:: + f(t) = + \begin{cases} + f_0 + \beta\, t^2 & \text{if vertex_zero is True,}\\ + f_1 + \beta\, (t_1 - t)^2 & \text{otherwise,} + \end{cases} + \quad\text{with}\quad + \beta = \frac{f_1 - f_0}{t_1^2} - To use a more general quadratic function, or an arbitrary - polynomial, use the function `scipy.signal.sweep_poly`. + The graph of the frequency f(t) is a parabola through :math:`(0, f_0)` and + :math:`(t_1, f_1)`. By default, the vertex of the parabola is at + :math:`(0, f_0)`. If `vertex_zero` is ``False``, then the vertex is at + :math:`(t_1, f_1)`. + To use a more general quadratic function, or an arbitrary + polynomial, use the function `scipy.signal.sweep_poly`. - logarithmic, log, lo: + 3. Parameter `method` in ``('logarithmic', 'log', 'lo')``: - ``f(t) = f0 * (f1/f0)**(t/t1)`` + .. math:: + f(t) = f_0 \left(\frac{f_1}{f_0}\right)^{t/t_1} - f0 and f1 must be nonzero and have the same sign. + :math:`f_0` and :math:`f_1` must be nonzero and have the same sign. + This signal is also known as a geometric or exponential chirp. - This signal is also known as a geometric or exponential chirp. + 4. Parameter `method` in ``('hyperbolic', 'hyp')``: - hyperbolic, hyp: + .. math:: + f(t) = \frac{\alpha}{\beta\, t + \gamma} \quad\text{with}\quad + \alpha = f_0 f_1 t_1, \ \beta = f_0 - f_1, \ \gamma = f_1 t_1 - ``f(t) = f0*f1*t1 / ((f0 - f1)*t + f1*t1)`` + :math:`f_0` and :math:`f_1` must be nonzero. - f0 and f1 must be nonzero. Examples -------- - The following will be used in the examples: + For the first example, a linear chirp ranging from 6 Hz to 1 Hz over 10 seconds is + plotted: >>> import numpy as np - >>> from scipy.signal import chirp, spectrogram + >>> from matplotlib.pyplot import tight_layout + >>> from scipy.signal import chirp, square, ShortTimeFFT + >>> from scipy.signal.windows import gaussian >>> import matplotlib.pyplot as plt - - For the first example, we'll plot the waveform for a linear chirp - from 6 Hz to 1 Hz over 10 seconds: - - >>> t = np.linspace(0, 10, 1500) - >>> w = chirp(t, f0=6, f1=1, t1=10, method='linear') - >>> plt.plot(t, w) - >>> plt.title("Linear Chirp, f(0)=6, f(10)=1") - >>> plt.xlabel('t (sec)') - >>> plt.show() - - For the remaining examples, we'll use higher frequency ranges, - and demonstrate the result using `scipy.signal.spectrogram`. - We'll use a 4 second interval sampled at 7200 Hz. - - >>> fs = 7200 - >>> T = 4 - >>> t = np.arange(0, int(T*fs)) / fs - - We'll use this function to plot the spectrogram in each example. - - >>> def plot_spectrogram(title, w, fs): - ... ff, tt, Sxx = spectrogram(w, fs=fs, nperseg=256, nfft=576) - ... fig, ax = plt.subplots() - ... ax.pcolormesh(tt, ff[:145], Sxx[:145], cmap='gray_r', - ... shading='gouraud') - ... ax.set_title(title) - ... ax.set_xlabel('t (sec)') - ... ax.set_ylabel('Frequency (Hz)') - ... ax.grid(True) ... - - Quadratic chirp from 1500 Hz to 250 Hz - (vertex of the parabolic curve of the frequency is at t=0): - - >>> w = chirp(t, f0=1500, f1=250, t1=T, method='quadratic') - >>> plot_spectrogram(f'Quadratic Chirp, f(0)=1500, f({T})=250', w, fs) - >>> plt.show() - - Quadratic chirp from 1500 Hz to 250 Hz - (vertex of the parabolic curve of the frequency is at t=T): - - >>> w = chirp(t, f0=1500, f1=250, t1=T, method='quadratic', - ... vertex_zero=False) - >>> plot_spectrogram(f'Quadratic Chirp, f(0)=1500, f({T})=250\\n' + - ... '(vertex_zero=False)', w, fs) + >>> N, T = 1000, 0.01 # number of samples and sampling interval for 10 s signal + >>> t = np.arange(N) * T # timestamps + ... + >>> x_lin = chirp(t, f0=6, f1=1, t1=10, method='linear') + ... + >>> fg0, ax0 = plt.subplots() + >>> ax0.set_title(r"Linear Chirp from $f(0)=6\,$Hz to $f(10)=1\,$Hz") + >>> ax0.set(xlabel="Time $t$ in Seconds", ylabel=r"Amplitude $x_\text{lin}(t)$") + >>> ax0.plot(t, x_lin) >>> plt.show() - Logarithmic chirp from 1500 Hz to 250 Hz: + The following four plots each show the short-time Fourier transform of a chirp + ranging from 45 Hz to 5 Hz with different values for the parameter `method` + (and `vertex_zero`): - >>> w = chirp(t, f0=1500, f1=250, t1=T, method='logarithmic') - >>> plot_spectrogram(f'Logarithmic Chirp, f(0)=1500, f({T})=250', w, fs) + >>> x_qu0 = chirp(t, f0=45, f1=5, t1=N*T, method='quadratic', vertex_zero=True) + >>> x_qu1 = chirp(t, f0=45, f1=5, t1=N*T, method='quadratic', vertex_zero=False) + >>> x_log = chirp(t, f0=45, f1=5, t1=N*T, method='logarithmic') + >>> x_hyp = chirp(t, f0=45, f1=5, t1=N*T, method='hyperbolic') + ... + >>> win = gaussian(50, std=12, sym=True) + >>> SFT = ShortTimeFFT(win, hop=2, fs=1/T, mfft=800, scale_to='magnitude') + >>> ts = ("'quadratic', vertex_zero=True", "'quadratic', vertex_zero=False", + ... "'logarithmic'", "'hyperbolic'") + >>> fg1, ax1s = plt.subplots(2, 2, sharex='all', sharey='all', + ... figsize=(6, 5), layout="constrained") + >>> for x_, ax_, t_ in zip([x_qu0, x_qu1, x_log, x_hyp], ax1s.ravel(), ts): + ... aSx = abs(SFT.stft(x_)) + ... im_ = ax_.imshow(aSx, origin='lower', aspect='auto', extent=SFT.extent(N), + ... cmap='plasma') + ... ax_.set_title(t_) + ... if t_ == "'hyperbolic'": + ... fg1.colorbar(im_, ax=ax1s, label='Magnitude $|S_z(t,f)|$') + >>> _ = fg1.supxlabel("Time $t$ in Seconds") # `_ =` is needed to pass doctests + >>> _ = fg1.supylabel("Frequency $f$ in Hertz") >>> plt.show() - Hyperbolic chirp from 1500 Hz to 250 Hz: + Finally, the short-time Fourier transform of a complex-valued linear chirp + ranging from -30 Hz to 30 Hz is depicted: - >>> w = chirp(t, f0=1500, f1=250, t1=T, method='hyperbolic') - >>> plot_spectrogram(f'Hyperbolic Chirp, f(0)=1500, f({T})=250', w, fs) + >>> z_lin = chirp(t, f0=-30, f1=30, t1=N*T, method="linear", complex=True) + >>> SFT.fft_mode = 'centered' # needed to work with complex signals + >>> aSz = abs(SFT.stft(z_lin)) + ... + >>> fg2, ax2 = plt.subplots() + >>> ax2.set_title(r"Linear Chirp from $-30\,$Hz to $30\,$Hz") + >>> ax2.set(xlabel="Time $t$ in Seconds", ylabel="Frequency $f$ in Hertz") + >>> im2 = ax2.imshow(aSz, origin='lower', aspect='auto', + ... extent=SFT.extent(N), cmap='viridis') + >>> fg2.colorbar(im2, label='Magnitude $|S_z(t,f)|$') >>> plt.show() + Note that using negative frequencies makes only sense with complex-valued signals. + Furthermore, the magnitude of the complex exponential function is one whereas the + magnitude of the real-valued cosine function is only 1/2. """ # 'phase' is computed in _chirp_phase, to make testing easier. - phase = _chirp_phase(t, f0, t1, f1, method, vertex_zero) - # Convert phi to radians. - phi *= pi / 180 - return cos(phase + phi) + phase = _chirp_phase(t, f0, t1, f1, method, vertex_zero) + np.deg2rad(phi) + return np.exp(1j*phase) if complex else np.cos(phase) def _chirp_phase(t, f0, t1, f1, method='linear', vertex_zero=True): diff --git a/scipy/signal/tests/test_waveforms.py b/scipy/signal/tests/test_waveforms.py index ae5b9cd15b0b..b30f7b9ceba8 100644 --- a/scipy/signal/tests/test_waveforms.py +++ b/scipy/signal/tests/test_waveforms.py @@ -73,6 +73,28 @@ def test_linear_freq_02(self): abserr = np.max(np.abs(f - chirp_linear(tf, f0, f1, t1))) assert abserr < 1e-6 + def test_linear_complex_power(self): + method = 'linear' + f0 = 1.0 + f1 = 2.0 + t1 = 1.0 + t = np.linspace(0, t1, 100) + w_real = waveforms.chirp(t, f0, t1, f1, method, complex=False) + w_complex = waveforms.chirp(t, f0, t1, f1, method, complex=True) + w_pwr_r = np.var(w_real) + w_pwr_c = np.var(w_complex) + + # Making sure that power of the real part is not affected with + # complex conversion operation + err = w_pwr_r - np.real(w_pwr_c) + + assert(err < 1e-6) + + def test_linear_complex_at_zero(self): + w = waveforms.chirp(t=0, f0=-10.0, f1=1.0, t1=1.0, method='linear', + complex=True) + xp_assert_close(w, 1.0+0.0j) # dtype must match + def test_quadratic_at_zero(self): w = waveforms.chirp(t=0, f0=1.0, f1=2.0, t1=1.0, method='quadratic') assert_almost_equal(w, 1.0) @@ -82,6 +104,11 @@ def test_quadratic_at_zero2(self): vertex_zero=False) assert_almost_equal(w, 1.0) + def test_quadratic_complex_at_zero(self): + w = waveforms.chirp(t=0, f0=-1.0, f1=2.0, t1=1.0, method='quadratic', + complex=True) + xp_assert_close(w, 1.0+0j) + def test_quadratic_freq_01(self): method = 'quadratic' f0 = 1.0