Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

ENH: signal: Add functionality of Complex Chirp waveform #17318

Merged
merged 1 commit into from
Nov 18, 2024
Merged

ENH: signal: Add functionality of Complex Chirp waveform #17318

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
192 changes: 102 additions & 90 deletions scipy/signal/_waveforms.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -284,135 +285,146 @@ 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
--------
sweep_poly

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):
Expand Down
27 changes: 27 additions & 0 deletions scipy/signal/tests/test_waveforms.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand All @@ -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):
Copy link
Member

Choose a reason for hiding this comment

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

We would need to test more known values to ensure the function behaves as expected and there are no corner cases it cannot handle.

Copy link
Contributor Author

Choose a reason for hiding this comment

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

Do you have any suggestions?

PS. Sorry, this fell through the cracks. If you can guide me through some example scenarios, I can write them.

Copy link
Contributor Author

Choose a reason for hiding this comment

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

Added one more test. I'll think about what else I can add.

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
Expand Down