Skip to content

Commit ec30b43

Browse files
DietBruDietBru
authored
Merge pull request scipy#17318 from anilgurses/complex_chirp
ENH: signal: Add functionality of Complex Chirp waveform
2 parents 4d2a0f3 + 96d56d9 commit ec30b43

File tree

2 files changed

+129
-90
lines changed

2 files changed

+129
-90
lines changed

‎scipy/signal/_waveforms.py

Lines changed: 102 additions & 90 deletions
Original file line numberDiff line numberDiff line change
@@ -257,8 +257,9 @@ def gausspulse(t, fc=1000, bw=0.5, bwr=-6, tpr=-60, retquad=False,
257257
return yI, yQ, yenv
258258

259259

260-
def chirp(t, f0, t1, f1, method='linear', phi=0, vertex_zero=True):
261-
"""Frequency-swept cosine generator.
260+
def chirp(t, f0, t1, f1, method='linear', phi=0, vertex_zero=True, *,
261+
complex=False):
262+
r"""Frequency-swept cosine generator.
262263
263264
In the following, 'Hz' should be interpreted as 'cycles per unit';
264265
there is no requirement here that the unit is one second. The
@@ -284,135 +285,146 @@ def chirp(t, f0, t1, f1, method='linear', phi=0, vertex_zero=True):
284285
This parameter is only used when `method` is 'quadratic'.
285286
It determines whether the vertex of the parabola that is the graph
286287
of the frequency is at t=0 or t=t1.
288+
complex : bool, optional
289+
This parameter creates a complex-valued analytic signal instead of a
290+
real-valued signal. It allows the use of complex baseband (in communications
291+
domain). Default is False.
292+
293+
.. versionadded:: 1.15.0
287294
288295
Returns
289296
-------
290297
y : ndarray
291-
A numpy array containing the signal evaluated at `t` with the
292-
requested time-varying frequency. More precisely, the function
293-
returns ``cos(phase + (pi/180)*phi)`` where `phase` is the integral
294-
(from 0 to `t`) of ``2*pi*f(t)``. ``f(t)`` is defined below.
298+
A numpy array containing the signal evaluated at `t` with the requested
299+
time-varying frequency. More precisely, the function returns
300+
``exp(1j*phase + 1j*(pi/180)*phi) if complex else cos(phase + (pi/180)*phi)``
301+
where `phase` is the integral (from 0 to `t`) of ``2*pi*f(t)``.
302+
The instantaneous frequency ``f(t)`` is defined below.
295303
296304
See Also
297305
--------
298306
sweep_poly
299307
300308
Notes
301309
-----
302-
There are four options for the `method`. The following formulas give
303-
the instantaneous frequency (in Hz) of the signal generated by
304-
`chirp()`. For convenience, the shorter names shown below may also be
305-
used.
306-
307-
linear, lin, li:
308-
309-
``f(t) = f0 + (f1 - f0) * t / t1``
310-
311-
quadratic, quad, q:
310+
There are four possible options for the parameter `method`, which have a (long)
311+
standard form and some allowed abbreviations. The formulas for the instantaneous
312+
frequency :math:`f(t)` of the generated signal are as follows:
312313
313-
The graph of the frequency f(t) is a parabola through (0, f0) and
314-
(t1, f1). By default, the vertex of the parabola is at (0, f0).
315-
If `vertex_zero` is False, then the vertex is at (t1, f1). The
316-
formula is:
314+
1. Parameter `method` in ``('linear', 'lin', 'li')``:
317315
318-
if vertex_zero is True:
316+
.. math::
317+
f(t) = f_0 + \beta\, t \quad\text{with}\quad
318+
\beta = \frac{f_1 - f_0}{t_1}
319319
320-
``f(t) = f0 + (f1 - f0) * t**2 / t1**2``
320+
Frequency :math:`f(t)` varies linearly over time with a constant rate
321+
:math:`\beta`.
321322
322-
else:
323+
2. Parameter `method` in ``('quadratic', 'quad', 'q')``:
323324
324-
``f(t) = f1 - (f1 - f0) * (t1 - t)**2 / t1**2``
325+
.. math::
326+
f(t) =
327+
\begin{cases}
328+
f_0 + \beta\, t^2 & \text{if vertex_zero is True,}\\
329+
f_1 + \beta\, (t_1 - t)^2 & \text{otherwise,}
330+
\end{cases}
331+
\quad\text{with}\quad
332+
\beta = \frac{f_1 - f_0}{t_1^2}
325333
326-
To use a more general quadratic function, or an arbitrary
327-
polynomial, use the function `scipy.signal.sweep_poly`.
334+
The graph of the frequency f(t) is a parabola through :math:`(0, f_0)` and
335+
:math:`(t_1, f_1)`. By default, the vertex of the parabola is at
336+
:math:`(0, f_0)`. If `vertex_zero` is ``False``, then the vertex is at
337+
:math:`(t_1, f_1)`.
338+
To use a more general quadratic function, or an arbitrary
339+
polynomial, use the function `scipy.signal.sweep_poly`.
328340
329-
logarithmic, log, lo:
341+
3. Parameter `method` in ``('logarithmic', 'log', 'lo')``:
330342
331-
``f(t) = f0 * (f1/f0)**(t/t1)``
343+
.. math::
344+
f(t) = f_0 \left(\frac{f_1}{f_0}\right)^{t/t_1}
332345
333-
f0 and f1 must be nonzero and have the same sign.
346+
:math:`f_0` and :math:`f_1` must be nonzero and have the same sign.
347+
This signal is also known as a geometric or exponential chirp.
334348
335-
This signal is also known as a geometric or exponential chirp.
349+
4. Parameter `method` in ``('hyperbolic', 'hyp')``:
336350
337-
hyperbolic, hyp:
351+
.. math::
352+
f(t) = \frac{\alpha}{\beta\, t + \gamma} \quad\text{with}\quad
353+
\alpha = f_0 f_1 t_1, \ \beta = f_0 - f_1, \ \gamma = f_1 t_1
338354
339-
``f(t) = f0*f1*t1 / ((f0 - f1)*t + f1*t1)``
355+
:math:`f_0` and :math:`f_1` must be nonzero.
340356
341-
f0 and f1 must be nonzero.
342357
343358
Examples
344359
--------
345-
The following will be used in the examples:
360+
For the first example, a linear chirp ranging from 6 Hz to 1 Hz over 10 seconds is
361+
plotted:
346362
347363
>>> import numpy as np
348-
>>> from scipy.signal import chirp, spectrogram
364+
>>> from matplotlib.pyplot import tight_layout
365+
>>> from scipy.signal import chirp, square, ShortTimeFFT
366+
>>> from scipy.signal.windows import gaussian
349367
>>> import matplotlib.pyplot as plt
350-
351-
For the first example, we'll plot the waveform for a linear chirp
352-
from 6 Hz to 1 Hz over 10 seconds:
353-
354-
>>> t = np.linspace(0, 10, 1500)
355-
>>> w = chirp(t, f0=6, f1=1, t1=10, method='linear')
356-
>>> plt.plot(t, w)
357-
>>> plt.title("Linear Chirp, f(0)=6, f(10)=1")
358-
>>> plt.xlabel('t (sec)')
359-
>>> plt.show()
360-
361-
For the remaining examples, we'll use higher frequency ranges,
362-
and demonstrate the result using `scipy.signal.spectrogram`.
363-
We'll use a 4 second interval sampled at 7200 Hz.
364-
365-
>>> fs = 7200
366-
>>> T = 4
367-
>>> t = np.arange(0, int(T*fs)) / fs
368-
369-
We'll use this function to plot the spectrogram in each example.
370-
371-
>>> def plot_spectrogram(title, w, fs):
372-
... ff, tt, Sxx = spectrogram(w, fs=fs, nperseg=256, nfft=576)
373-
... fig, ax = plt.subplots()
374-
... ax.pcolormesh(tt, ff[:145], Sxx[:145], cmap='gray_r',
375-
... shading='gouraud')
376-
... ax.set_title(title)
377-
... ax.set_xlabel('t (sec)')
378-
... ax.set_ylabel('Frequency (Hz)')
379-
... ax.grid(True)
380368
...
381-
382-
Quadratic chirp from 1500 Hz to 250 Hz
383-
(vertex of the parabolic curve of the frequency is at t=0):
384-
385-
>>> w = chirp(t, f0=1500, f1=250, t1=T, method='quadratic')
386-
>>> plot_spectrogram(f'Quadratic Chirp, f(0)=1500, f({T})=250', w, fs)
387-
>>> plt.show()
388-
389-
Quadratic chirp from 1500 Hz to 250 Hz
390-
(vertex of the parabolic curve of the frequency is at t=T):
391-
392-
>>> w = chirp(t, f0=1500, f1=250, t1=T, method='quadratic',
393-
... vertex_zero=False)
394-
>>> plot_spectrogram(f'Quadratic Chirp, f(0)=1500, f({T})=250\\n' +
395-
... '(vertex_zero=False)', w, fs)
369+
>>> N, T = 1000, 0.01 # number of samples and sampling interval for 10 s signal
370+
>>> t = np.arange(N) * T # timestamps
371+
...
372+
>>> x_lin = chirp(t, f0=6, f1=1, t1=10, method='linear')
373+
...
374+
>>> fg0, ax0 = plt.subplots()
375+
>>> ax0.set_title(r"Linear Chirp from $f(0)=6\,$Hz to $f(10)=1\,$Hz")
376+
>>> ax0.set(xlabel="Time $t$ in Seconds", ylabel=r"Amplitude $x_\text{lin}(t)$")
377+
>>> ax0.plot(t, x_lin)
396378
>>> plt.show()
397379
398-
Logarithmic chirp from 1500 Hz to 250 Hz:
380+
The following four plots each show the short-time Fourier transform of a chirp
381+
ranging from 45 Hz to 5 Hz with different values for the parameter `method`
382+
(and `vertex_zero`):
399383
400-
>>> w = chirp(t, f0=1500, f1=250, t1=T, method='logarithmic')
401-
>>> plot_spectrogram(f'Logarithmic Chirp, f(0)=1500, f({T})=250', w, fs)
384+
>>> x_qu0 = chirp(t, f0=45, f1=5, t1=N*T, method='quadratic', vertex_zero=True)
385+
>>> x_qu1 = chirp(t, f0=45, f1=5, t1=N*T, method='quadratic', vertex_zero=False)
386+
>>> x_log = chirp(t, f0=45, f1=5, t1=N*T, method='logarithmic')
387+
>>> x_hyp = chirp(t, f0=45, f1=5, t1=N*T, method='hyperbolic')
388+
...
389+
>>> win = gaussian(50, std=12, sym=True)
390+
>>> SFT = ShortTimeFFT(win, hop=2, fs=1/T, mfft=800, scale_to='magnitude')
391+
>>> ts = ("'quadratic', vertex_zero=True", "'quadratic', vertex_zero=False",
392+
... "'logarithmic'", "'hyperbolic'")
393+
>>> fg1, ax1s = plt.subplots(2, 2, sharex='all', sharey='all',
394+
... figsize=(6, 5), layout="constrained")
395+
>>> for x_, ax_, t_ in zip([x_qu0, x_qu1, x_log, x_hyp], ax1s.ravel(), ts):
396+
... aSx = abs(SFT.stft(x_))
397+
... im_ = ax_.imshow(aSx, origin='lower', aspect='auto', extent=SFT.extent(N),
398+
... cmap='plasma')
399+
... ax_.set_title(t_)
400+
... if t_ == "'hyperbolic'":
401+
... fg1.colorbar(im_, ax=ax1s, label='Magnitude $|S_z(t,f)|$')
402+
>>> _ = fg1.supxlabel("Time $t$ in Seconds") # `_ =` is needed to pass doctests
403+
>>> _ = fg1.supylabel("Frequency $f$ in Hertz")
402404
>>> plt.show()
403405
404-
Hyperbolic chirp from 1500 Hz to 250 Hz:
406+
Finally, the short-time Fourier transform of a complex-valued linear chirp
407+
ranging from -30 Hz to 30 Hz is depicted:
405408
406-
>>> w = chirp(t, f0=1500, f1=250, t1=T, method='hyperbolic')
407-
>>> plot_spectrogram(f'Hyperbolic Chirp, f(0)=1500, f({T})=250', w, fs)
409+
>>> z_lin = chirp(t, f0=-30, f1=30, t1=N*T, method="linear", complex=True)
410+
>>> SFT.fft_mode = 'centered' # needed to work with complex signals
411+
>>> aSz = abs(SFT.stft(z_lin))
412+
...
413+
>>> fg2, ax2 = plt.subplots()
414+
>>> ax2.set_title(r"Linear Chirp from $-30\,$Hz to $30\,$Hz")
415+
>>> ax2.set(xlabel="Time $t$ in Seconds", ylabel="Frequency $f$ in Hertz")
416+
>>> im2 = ax2.imshow(aSz, origin='lower', aspect='auto',
417+
... extent=SFT.extent(N), cmap='viridis')
418+
>>> fg2.colorbar(im2, label='Magnitude $|S_z(t,f)|$')
408419
>>> plt.show()
409420
421+
Note that using negative frequencies makes only sense with complex-valued signals.
422+
Furthermore, the magnitude of the complex exponential function is one whereas the
423+
magnitude of the real-valued cosine function is only 1/2.
410424
"""
411425
# 'phase' is computed in _chirp_phase, to make testing easier.
412-
phase = _chirp_phase(t, f0, t1, f1, method, vertex_zero)
413-
# Convert phi to radians.
414-
phi *= pi / 180
415-
return cos(phase + phi)
426+
phase = _chirp_phase(t, f0, t1, f1, method, vertex_zero) + np.deg2rad(phi)
427+
return np.exp(1j*phase) if complex else np.cos(phase)
416428

417429

418430
def _chirp_phase(t, f0, t1, f1, method='linear', vertex_zero=True):

‎scipy/signal/tests/test_waveforms.py

Lines changed: 27 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -73,6 +73,28 @@ def test_linear_freq_02(self):
7373
abserr = np.max(np.abs(f - chirp_linear(tf, f0, f1, t1)))
7474
assert abserr < 1e-6
7575

76+
def test_linear_complex_power(self):
77+
method = 'linear'
78+
f0 = 1.0
79+
f1 = 2.0
80+
t1 = 1.0
81+
t = np.linspace(0, t1, 100)
82+
w_real = waveforms.chirp(t, f0, t1, f1, method, complex=False)
83+
w_complex = waveforms.chirp(t, f0, t1, f1, method, complex=True)
84+
w_pwr_r = np.var(w_real)
85+
w_pwr_c = np.var(w_complex)
86+
87+
# Making sure that power of the real part is not affected with
88+
# complex conversion operation
89+
err = w_pwr_r - np.real(w_pwr_c)
90+
91+
assert(err < 1e-6)
92+
93+
def test_linear_complex_at_zero(self):
94+
w = waveforms.chirp(t=0, f0=-10.0, f1=1.0, t1=1.0, method='linear',
95+
complex=True)
96+
xp_assert_close(w, 1.0+0.0j) # dtype must match
97+
7698
def test_quadratic_at_zero(self):
7799
w = waveforms.chirp(t=0, f0=1.0, f1=2.0, t1=1.0, method='quadratic')
78100
assert_almost_equal(w, 1.0)
@@ -82,6 +104,11 @@ def test_quadratic_at_zero2(self):
82104
vertex_zero=False)
83105
assert_almost_equal(w, 1.0)
84106

107+
def test_quadratic_complex_at_zero(self):
108+
w = waveforms.chirp(t=0, f0=-1.0, f1=2.0, t1=1.0, method='quadratic',
109+
complex=True)
110+
xp_assert_close(w, 1.0+0j)
111+
85112
def test_quadratic_freq_01(self):
86113
method = 'quadratic'
87114
f0 = 1.0

0 commit comments

Comments
 (0)