diff --git a/.gitignore b/.gitignore index 704b1e1..ecba673 100644 --- a/.gitignore +++ b/.gitignore @@ -52,6 +52,7 @@ coverage.xml # Sphinx documentation doc/_build/ +doc/apiref/*.rst # PyBuilder target/ diff --git a/doc/apiref/modules.rst b/doc/apiref/modules.rst deleted file mode 100644 index 8ebe0b5..0000000 --- a/doc/apiref/modules.rst +++ /dev/null @@ -1,7 +0,0 @@ -tftb -==== - -.. toctree:: - :maxdepth: 4 - - tftb diff --git a/doc/apiref/tftb.generators.rst b/doc/apiref/tftb.generators.rst deleted file mode 100644 index 43bb0a8..0000000 --- a/doc/apiref/tftb.generators.rst +++ /dev/null @@ -1,69 +0,0 @@ -tftb.generators package -======================= - -Subpackages ------------ - -.. toctree:: - - tftb.generators.tests - -Submodules ----------- - -tftb.generators.amplitude_modulated module ------------------------------------------- - -.. automodule:: tftb.generators.amplitude_modulated - :members: - :undoc-members: - :show-inheritance: - -tftb.generators.analytic_signals module ---------------------------------------- - -.. automodule:: tftb.generators.analytic_signals - :members: - :undoc-members: - :show-inheritance: - -tftb.generators.frequency_modulated module ------------------------------------------- - -.. automodule:: tftb.generators.frequency_modulated - :members: - :undoc-members: - :show-inheritance: - -tftb.generators.misc module ---------------------------- - -.. automodule:: tftb.generators.misc - :members: - :undoc-members: - :show-inheritance: - -tftb.generators.noise module ----------------------------- - -.. automodule:: tftb.generators.noise - :members: - :undoc-members: - :show-inheritance: - -tftb.generators.utils module ----------------------------- - -.. automodule:: tftb.generators.utils - :members: - :undoc-members: - :show-inheritance: - - -Module contents ---------------- - -.. automodule:: tftb.generators - :members: - :undoc-members: - :show-inheritance: diff --git a/doc/apiref/tftb.generators.tests.rst b/doc/apiref/tftb.generators.tests.rst deleted file mode 100644 index 2412638..0000000 --- a/doc/apiref/tftb.generators.tests.rst +++ /dev/null @@ -1,62 +0,0 @@ -tftb.generators.tests package -============================= - -Submodules ----------- - -tftb.generators.tests.test_amplitude_modulations module -------------------------------------------------------- - -.. automodule:: tftb.generators.tests.test_amplitude_modulations - :members: - :undoc-members: - :show-inheritance: - -tftb.generators.tests.test_analytic_signals module --------------------------------------------------- - -.. automodule:: tftb.generators.tests.test_analytic_signals - :members: - :undoc-members: - :show-inheritance: - -tftb.generators.tests.test_frequency_modulations module -------------------------------------------------------- - -.. automodule:: tftb.generators.tests.test_frequency_modulations - :members: - :undoc-members: - :show-inheritance: - -tftb.generators.tests.test_misc module --------------------------------------- - -.. automodule:: tftb.generators.tests.test_misc - :members: - :undoc-members: - :show-inheritance: - -tftb.generators.tests.test_noise module ---------------------------------------- - -.. automodule:: tftb.generators.tests.test_noise - :members: - :undoc-members: - :show-inheritance: - -tftb.generators.tests.test_utils module ---------------------------------------- - -.. automodule:: tftb.generators.tests.test_utils - :members: - :undoc-members: - :show-inheritance: - - -Module contents ---------------- - -.. automodule:: tftb.generators.tests - :members: - :undoc-members: - :show-inheritance: diff --git a/doc/apiref/tftb.processing.rst b/doc/apiref/tftb.processing.rst deleted file mode 100644 index 5893603..0000000 --- a/doc/apiref/tftb.processing.rst +++ /dev/null @@ -1,102 +0,0 @@ -tftb.processing package -======================= - -Submodules ----------- - -tftb.processing.affine module ------------------------------ - -.. automodule:: tftb.processing.affine - :members: - :undoc-members: - :show-inheritance: - -tftb.processing.ambiguity module --------------------------------- - -.. automodule:: tftb.processing.ambiguity - :members: - :undoc-members: - :show-inheritance: - -tftb.processing.base module ---------------------------- - -.. automodule:: tftb.processing.base - :members: - :undoc-members: - :show-inheritance: - -tftb.processing.cohen module ----------------------------- - -.. automodule:: tftb.processing.cohen - :members: - :undoc-members: - :show-inheritance: - -tftb.processing.freq_domain module ----------------------------------- - -.. automodule:: tftb.processing.freq_domain - :members: - :undoc-members: - :show-inheritance: - -tftb.processing.linear module ------------------------------ - -.. automodule:: tftb.processing.linear - :members: - :undoc-members: - :show-inheritance: - -tftb.processing.plotifl module ------------------------------- - -.. automodule:: tftb.processing.plotifl - :members: - :undoc-members: - :show-inheritance: - -tftb.processing.postprocessing module -------------------------------------- - -.. automodule:: tftb.processing.postprocessing - :members: - :undoc-members: - :show-inheritance: - -tftb.processing.reassigned module ---------------------------------- - -.. automodule:: tftb.processing.reassigned - :members: - :undoc-members: - :show-inheritance: - -tftb.processing.time_domain module ----------------------------------- - -.. automodule:: tftb.processing.time_domain - :members: - :undoc-members: - :show-inheritance: - -tftb.processing.utils module ----------------------------- - -.. automodule:: tftb.processing.utils - :members: - :undoc-members: - :show-inheritance: - - -Module contents ---------------- - -.. automodule:: tftb.processing - :members: - :undoc-members: - :show-inheritance: diff --git a/doc/apiref/tftb.rst b/doc/apiref/tftb.rst deleted file mode 100644 index 9bfff10..0000000 --- a/doc/apiref/tftb.rst +++ /dev/null @@ -1,31 +0,0 @@ -tftb package -============ - -Subpackages ------------ - -.. toctree:: - - tftb.generators - tftb.processing - tftb.tests - -Submodules ----------- - -tftb.utils module ------------------ - -.. automodule:: tftb.utils - :members: - :undoc-members: - :show-inheritance: - - -Module contents ---------------- - -.. automodule:: tftb - :members: - :undoc-members: - :show-inheritance: diff --git a/doc/apiref/tftb.tests.rst b/doc/apiref/tftb.tests.rst deleted file mode 100644 index 0e5b7b7..0000000 --- a/doc/apiref/tftb.tests.rst +++ /dev/null @@ -1,30 +0,0 @@ -tftb.tests package -================== - -Submodules ----------- - -tftb.tests.base module ----------------------- - -.. automodule:: tftb.tests.base - :members: - :undoc-members: - :show-inheritance: - -tftb.tests.test_utils module ----------------------------- - -.. automodule:: tftb.tests.test_utils - :members: - :undoc-members: - :show-inheritance: - - -Module contents ---------------- - -.. automodule:: tftb.tests - :members: - :undoc-members: - :show-inheritance: diff --git a/doc/index.rst b/doc/index.rst index 2b04fb5..1c459e1 100644 --- a/doc/index.rst +++ b/doc/index.rst @@ -12,6 +12,7 @@ Contents: :maxdepth: 2 introduction + nonstationary_signals auto_gallery api diff --git a/doc/misc_plots/nonstationary_phase_plot.py b/doc/misc_plots/nonstationary_phase_plot.py new file mode 100644 index 0000000..2afaa33 --- /dev/null +++ b/doc/misc_plots/nonstationary_phase_plot.py @@ -0,0 +1,10 @@ +from tftb.generators import fmlin, amgauss +import numpy as np +import matplotlib.pyplot as plt + +y_nonstat, _ = fmlin(2048) # Already analytic, no need of Hilbert transorm +y_nonstat *= amgauss(2048) +plt.plot(np.real(y_nonstat), np.imag(y_nonstat)) +plt.xlabel("Real part") +plt.ylabel("Imaginary part") +plt.show() diff --git a/doc/misc_plots/stationary_phase_plot.py b/doc/misc_plots/stationary_phase_plot.py new file mode 100644 index 0000000..bddc198 --- /dev/null +++ b/doc/misc_plots/stationary_phase_plot.py @@ -0,0 +1,17 @@ +import numpy as np +import matplotlib.pyplot as plt +from scipy.signal import hilbert + +fs = 32768 +ts = np.linspace(0, 1, fs) +y1 = np.sin(2 * np.pi * 697 * ts) +y2 = np.sin(2 * np.pi * 1336 * ts) +y = (y1 + y2) / 2 + + +y = y[:(fs / 16)] +y_analytic = hilbert(y) +plt.plot(np.real(y_analytic), np.imag(y_analytic)) +plt.xlabel("Real part") +plt.ylabel("Imaginary part") +plt.show() diff --git a/doc/misc_plots/touchtone.py b/doc/misc_plots/touchtone.py new file mode 100644 index 0000000..4f8bbf1 --- /dev/null +++ b/doc/misc_plots/touchtone.py @@ -0,0 +1,25 @@ +#! /usr/bin/env python +# -*- coding: utf-8 -*- +# vim:fenc=utf-8 +# +# Cube26 product code +# +# (C) Copyright 2015 Cube26 Software Pvt Ltd +# All right reserved. +# +# This file is confidential and NOT open source. Do not distribute. +# + +""" + +""" +import numpy as np +import matplotlib.pyplot as plt +fs = 32768 +ts = np.linspace(0, 1, fs) +y1 = np.sin(2 * np.pi * 697 * ts) +y2 = np.sin(2 * np.pi * 1336 * ts) +y = (y1 + y2) / 2 +plt.plot(ts, y) +plt.xlim(0, 0.1) +plt.show() diff --git a/doc/misc_plots/touchtone_mean_convolve.py b/doc/misc_plots/touchtone_mean_convolve.py new file mode 100644 index 0000000..96bc26a --- /dev/null +++ b/doc/misc_plots/touchtone_mean_convolve.py @@ -0,0 +1,27 @@ +#! /usr/bin/env python +# -*- coding: utf-8 -*- +# vim:fenc=utf-8 +# +# Cube26 product code +# +# (C) Copyright 2015 Cube26 Software Pvt Ltd +# All right reserved. +# +# This file is confidential and NOT open source. Do not distribute. +# + +""" + +""" +import numpy as np +import matplotlib.pyplot as plt +fs = 32768 +ts = np.linspace(0, 1, fs) +y1 = np.sin(2 * np.pi * 697 * ts) +y2 = np.sin(2 * np.pi * 1336 * ts) +y = (y1 + y2) / 2 + +Y = np.fft.fft(y) +freqs = np.arange(fs) +plt.plot(freqs, np.abs(Y)) +plt.show() diff --git a/doc/misc_plots/uncertainty_example_plot.py b/doc/misc_plots/uncertainty_example_plot.py new file mode 100644 index 0000000..239ed99 --- /dev/null +++ b/doc/misc_plots/uncertainty_example_plot.py @@ -0,0 +1,12 @@ +import numpy as np +import matplotlib.pyplot as plt +f1, f2 = 500, 1000 +t1, t2 = 0.192, 0.196 +f_sample = 8000 +n_points = 2048 +ts = np.arange(n_points, dtype=float) / f_sample +signal = np.sin(2 * np.pi * f1 * ts) + np.sin(2 * np.pi * f2 * ts) +signal[int(t1 * f_sample) - 1] += 3 +signal[int(t2 * f_sample) - 1] += 3 +plt.plot(ts, signal) +plt.show() diff --git a/doc/misc_plots/uncertainty_stft.py b/doc/misc_plots/uncertainty_stft.py new file mode 100644 index 0000000..5d1c13a --- /dev/null +++ b/doc/misc_plots/uncertainty_stft.py @@ -0,0 +1,29 @@ +import numpy as np +import matplotlib.pyplot as plt +from tftb.processing import ShortTimeFourierTransform +f1, f2 = 500, 1000 +t1, t2 = 0.192, 0.196 +f_sample = 8000 +n_points = 2048 +ts = np.arange(n_points, dtype=float) / f_sample +signal = np.sin(2 * np.pi * f1 * ts) + np.sin(2 * np.pi * f2 * ts) +signal[int(t1 * f_sample) - 1] += 3 +signal[int(t2 * f_sample) - 1] += 3 + +wlengths = [2, 4, 8, 16] +nf = [(w * 0.001 * f_sample) + 1 for w in wlengths] +fig = plt.figure() +extent = [0, ts.max(), 0, 2000] +for i, wlen in enumerate(wlengths): + window = np.ones((nf[i],), dtype=float) + stft = ShortTimeFourierTransform(signal, fwindow=window) + stft.run() + ax = fig.add_subplot(4, 1, i + 1) + stft.plot(ax=ax, default_annotation=False, show=False, + extent=extent) + ax.set_yticklabels([]) + if i != 3: + ax.set_xticklabels([]) + ax.set_title("window duration = {} ms".format(wlen)) +plt.subplots_adjust(hspace=0.5) +plt.show() diff --git a/doc/nonstationary_signals.rst b/doc/nonstationary_signals.rst new file mode 100644 index 0000000..a15774d --- /dev/null +++ b/doc/nonstationary_signals.rst @@ -0,0 +1,382 @@ +====================== +Non Stationary Signals +====================== + +Frequency Domain Representations +-------------------------------- + +The most straightforward frequency domain representation of a signal is +obtained by the `Fourier transform `_ +of the signal, defined as: + + .. math:: + + X(\nu) = \int_{-\infty}^{\infty}x(t)e^{-j2\pi\nu t}dt + +The spectrum :math:`X(\nu)` is perfectly valid, but the Fourier transform is +essentially an integral over time. Thus, we lose all information that varies +with time. All we can tell from the spectrum is that the signal has two +distinct frequency components. In other words, we can comment on what happens +a signal, not when it happens. Consider a song as the signal under +consideration. If you were not interested in time, the whole point of +processing that signal would be lost. Rhythm and timing are the very heart of +good music, after all. In this case, we want to know when the drums kicked i +, as well as what notes were being played on the guitar. If we perform only +frequency analysis, all time information would be lost and the only +information we would have would be about what frequencies were played in the +song, and what their respective amplitudes were, averaged over the duration +of the entire song. So even if the drums stop playing after the second stanza, +the frequency spectrum would show them playing throughout the song. +Conversely, if we were only interested in the time information, we would be +hardly better off than simply listening to the song. + +The solution to this problem is essentially time-frequency analysis which is +a field that deals with signal processing in both time and frequency domain. +It consists of a collection of methods that allow us to make tradeoffs +between time and frequency processing of a signal, depending on what makes +more sense for a particular application, as we shall see through the rest of +this tutorial. + +The Heisenberg-Gabor Inequality +------------------------------- + +Before delving into joint time frequency representations, it is necessary to +understand that any signal is characterized in the time-frequncy space by two +quantities: + +1. The *mean* position of the signal, defined as pair of two figures: average + time (:math:`t_{m}`) and average frequency (:math:`\nu_{m}`) +2. The energy localization of the signal in the time-frequency space, whose + area is proportional to the *Time-Bandwidth product*. An important + constraint related to this quantity is called the Heisenberg-Gabor + inequality, which we shall explore later in this section. + +If a signal :math:`x(t)` has finite energy, i.e. + + .. math:: + + E_{x} = \int_{-\infty}^{\infty} \left|x(t)\right|^{2} dt < \infty + +then the time and frequency domain energies of the signal can be considered as +probability distributions, and their respective means and standard deviations +can be used to estimate the time and frequency localizations and dispersions of +the signal. + +* Average time: + :math:`t_{m} = \frac{1}{E_{x}}\int_{-\infty}^{\infty}t\left|x(t)\right|^{2}dt` +* Average frequency: + :math:`\nu_{m} = \frac{1}{E_{x}}\int_{-\infty}^{\infty}\nu\left|X(\nu)\right|^{2}d\nu` +* Time spreading: + :math:`T^{2} = \frac{4\pi}{E_{x}}\int_{-\infty}^{\infty}(t-t_{m})^{2}\left|x(t)\right|^{2}dt` +* Frequency spreading: + :math:`B^{2} = \frac{4\pi}{E_{x}}\int_{-\infty}^{\infty}(\nu-\nu_{m})^{2}\left|X(\nu)\right|^{2}d\nu` + +Let's play around with these values with some examples. + +Example: Time and Frequency Localizations +````````````````````````````````````````` + +Time and frequency localizations can be calculated with the functions +``tftb.processing.loctime`` and ``tftb.processing.locfreq``. Consider a linear +chirp with Gaussian amplitude modulation as an example, shown below: + + .. plot:: _gallery/plot_2_2_1_time_freq_localization.py + + >>> from tftb.generators import fmlin, amgauss + >>> from tftb.processing import loctime, locfreq + >>> sig = fmlin(256)[0] * amgauss(256) + >>> t_mean, time_spreading = loctime(sig) + >>> print t_mean, time_spreading + 127.0 32.0 + >>> f_mean, freq_spreading = locfreq(sig) + >>> print f_mean, freq_spreading + 0.249019607843 0.0700964323482 + +The time-bandwidth product of the signal can be obtained by multiplying the +``time_spreading`` and ``frequency_spreading`` variables in the snippet above. +An important inequality concering the time bandwidth product is called the +`uncertainty principle `_. +which states that a signal cannot be localized simultaneously in both time and +frequency with arbitrariy high resolution. The next example demonstrates this +concept. + +Example: The Uncertainty Principle +`````````````````````````````````` + +The uncertainty principle is a very manifest limitation of the Fourier +transform. Consider the signal shown here:: + + >>> import numpy as np + >>> import matplotlib.pyplot as plt + >>> f1, f2 = 500, 1000 + >>> t1, t2 = 0.192, 0.196 + >>> f_sample = 8000 + >>> n_points = 2048 + >>> ts = np.arange(n_points, dtype=float) / f_sample + >>> signal = np.sin(2 * np.pi * f1 * ts) + np.sin(2 * np.pi * f2 * ts) + >>> signal[int(t1 * f_sample) - 1] += 3 + >>> signal[int(t2 * f_sample) - 1] += 3 + >>> plt.plot(ts, signal) + >>> plt.show() + +.. plot:: misc_plots/uncertainty_example_plot.py + +It is a sum of two sinusiodal signals of frequencies 500 Hz and 1000 Hz. It has +two spikes at :math:`t` = 0.192s and :math:`t` = 0.196s. The purpose of a time frequency +distribution would be to clearly identify both the frequencies and both the spikes, +thus resolving events in both frequency and time. Let's check out the spectrograms of +of the signal with four different window lengths: + +.. plot:: misc_plots/uncertainty_stft.py + +As can be clearly seen, resolution in time and frequency +cannot be obtained simultaneously. In the last (bottom) image, where the +window length is high, the STFT manages to discriminate between frequencies +of 500 Hz and 1000 Hz very clearly, but the time resolution between the +events at t = 0.192 s and t = 0.196 s is ambiguous. As we reduce the length +of the window function, the resolution between the time events goes on +becoming better, but only at the cost of resolution in frequencies. + +Informally, the uncertainty principle states +that arbitrarily high resolution cannot be obtained in both time and frequency. +This is a consequence of the definition of the Fourier transform. The +definition insists that a signal be represented as a weighted sum of sinusoids, +and therefore identifies frequency information that is globally prevalent. As +a workaround to this interpretation, we use the STFT which performs the +Fourier transform on limited periods of the signals. Mathematically this +uncertainty can be quantified with the Heisenberg-Gabor Inequality (also +sometimes called the Gabor limit): + +.. topic:: Heisenberg - Gabor Inequality + + If :math:`T` and :math:`B` are standard deviations of the time + characteristics and the bandwidth respectively of a signal :math:`s(t)`, + then + + .. math:: + + TB ≥ 1 + +The expression states that the time-bandwidth product of a signal is lower +bounded by unity. Gaussian functions satisfy the equality condition in the +equation. This can be verified as follows:: + + >>> from tftb.generators import fmconst, amgauss + >>> x = gen.amgauss(128) * gen.fmconst(128)[0] + >>> plot(real(x)) + +.. plot:: + + from tftb.generators import fmconst, amgauss + import matplotlib.pyplot as plt + from numpy import real + x = amgauss(128) * fmconst(128)[0] + plt.plot(real(x)) + plt.grid() + plt.xlim(0, 128) + plt.title("Gaussian amplitude modulation") + plt.show() + +.. code-block:: python + + >>> from tftb.processing import loctime, locfreq + >>> time_mean, time_duration = loctime(x) + >>> freq_center, bandwidth = locfreq(x) + >>> time_duration * bandwidth + 1.0 + +A remarkably insightful commentary on the Uncertainty principle is provided +in [1]_, which states that the Uncertainty principle is a statement about two +variables whose associated operators do not mutually commute. This helps us +apply the Uncertainty principle in signal processing in the same way as in +quantum physics. + +Instantaneous Frequency +----------------------- + +An alternative way to localize a signal in time and frequency is its +instantaneous frequency. Instantaneous frequencies are defined for analytic +signals, which are defined as follows:: + + .. math:: + + x_{a}(t) = x(t) + jH(x(t)) + +where :math:`x(t)` is a real valued time domain signal and `H` denotes the +Hilbert transform (``scipy.signal.hilbert``). From this defition of the +analytic signal, the following quantities can be derived: + +* Instantaneous amplitude + :math:`a(t) = \left|x_{a}(t)\right|` +* Instantaneous frequency + :math:`f(t) = \frac{1}{2\pi}\frac{d}{dt}arg(x_{a}(t))` + +An implementation of these functions can be found in ``tftb.processing.instfreq`` + +Example: Instantaneous Frequency +```````````````````````````````` + + >>> from tftb.processign import inst_freq, plotifl + >>> signal, _ = fmlin(256) + >>> time_samples = np.arange(3, 257) + >>> ifr = inst_freq(signal)[0] + >>> plotifl(time_samples, ifr) + + .. plot:: _gallery/plot_2_3_instantaneous_frequency.py + + +Group Delay +----------- + +The frequency domain equivalent of instantaneous frequency is called group +delay, which localizes time characteristics of a signal as function of the +frequency. + + .. math:: + + t_{x}(\nu) = -\frac{1}{2\pi}\frac{d}{d\nu}arg(X_{a}(\nu)) + + +Example: Group Delay +```````````````````` + +The group delay of the signal in the previous example can be obtained as +follows + + >>> from tftb.processign import group_delay + >>> fnorm = np.linspace(0, .5, 10) + >>> gd = group_delay(signal, fnorm) + >>> plt.plot(gd, fnorm) + + .. plot:: _gallery/plot_2_4_group_delay.py + + +Example: Comparison of Instantaneous Frequency and Group Delay +`````````````````````````````````````````````````````````````` + +Ideally, for a signal localized perfectly in time and frequency, its +instantaneous frequency and group delay would be expected to be identical. +However, mathematically they are two different fuctions in the time-frequency +space, and only coincide for signals with high time-bandwidth products. This +makes sense, since a high time-bandwidth product implies that the +signal would be pushed away from the Heisenberg-Gabor inequality, thereby +leading to lesser ambiguity. Consider the following example, where we construct +two signals - one with a high time-bandwidth product, and one with a low one - +and then estimate their respective instantaneous frequencies and group delays. + + >>> # generate a signal with a high TB + >>> time_instants = np.arange(2, 256) + >>> sig1 = amgauss(256, 128, 90) * fmlin(256)[0] + >>> tm, T1 = loctime(sig1) + >>> fm, B1 = locfreq(sig1) + >>> print T1 * B1 + 15.9138 + >>> ifr1 = inst_freq(sig1, time_instants)[0] + >>> f1 = np.linspace(0, 0.5 - 1.0 / 256, 256) + >>> gd1 = group_delay(sig1, f1) + >>> + >>> plt.subplot(211) + >>> plt.plot(time_instants, ifr1, '*', label='inst_freq') + >>> plt.plot(gd1, f1, '-', label='group delay') + >>> + >>> # generate a signal with low TB + >>> sig2 = amgauss(256, 128, 30) * fmlin(256, 0.2, 0.4)[0] + >>> tm, T2 = loctime(sig2) + >>> fm, B2 = locfreq(sig2) + >>> print T2 * B2 + 1.224 + >>> ifr2 = inst_freq(sig2, time_instants)[0] + >>> f2 = np.linspace(0.02, 0.4, 256) + >>> gd2 = group_delay(sig2, f2) + >>> + >>> plt.subplot(212) + >>> plt.plot(time_instants, ifr2, '*', label='inst_freq') + >>> plt.plot(gd2, f2, '-', label='group delay') + + .. plot:: _gallery/plot_2_4_group_delay.py + + +A Note on Stationarity +---------------------- + +As per Wikipedia, a "stationary" process is one whose joint probability +distribution does not change with time (or space). Let's try and see what a +stationary process looks like. Consider a signal generated as follows: + + >>> import numpy as np + >>> import matplotlib.pyplot as plt + >>> fs = 32768 + >>> ts = np.linspace(0, 1, fs) + >>> y1 = np.sin(2 * np.pi * 697 * ts) + >>> y2 = np.sin(2 * np.pi * 1336 * ts) + >>> y = (y1 + y2) / 2 + >>> plt.plot(ts, y) + >>> plt.xlim(0, 0.1) + >>> plt.show() + + .. plot:: misc_plots/touchtone.py + +The plot shows a slice of a touchtone signal with the duration of a tenth of a +second. This is the signal you hear when you press the "2" key on a telephone +dialing pad. (You can save the generated signal as a WAV file as follows: + + >>> from scipy.io import wavfile + >>> wavfile.write("tone.wav", fs, y) + +and listen to the file ``tone.wav`` with your favourite music player.) + +Since the signal is composed of two sinusoids, ``y1`` and ``y2``, we would +expect it to be stationary. Let's try and assert this qualitatively. Let's try +to plot the signal in its phase space. In order to do so, we will first need to +construct an analytic representation of the signal. For simplicity, we shall only +consider a part of the original signal + + + >>> y = y[:(fs / 16)] + >>> y_analytic = hilbert(y) + >>> plt.plot(np.real(y_analytic), np.imag(y_analytic)) + >>> plt.xlabel("Real part") + >>> plt.ylabel("Imaginary part") + >>> plt.show() + + .. plot:: misc_plots/stationary_phase_plot.py + +This visualization can be interpreted as follows. Imagine that there is a +vector centered at the origin of this plot which traces out the signal as it +rotates about the origin. Then, at any time :math:`t`, the angle which the +vector makes with the real axis is the instantaneous phase of the signal, +:math:`\theta(t)`. The angular speed with which the phasor rotates is the +instantaneous frequency of the signal: + + .. math:: + + \omega(t) = \frac{d}{dt} \theta(t) + +Now, let's compare this phase plot with that of a known nonstationary signal. + + + >>> from tftb.generators import fmlin, amgauss + >>> y_ns, _ = fmlin(2048) # Already analytic, no need of Hilbert transorm + >>> y_nonstat = y_ns * amgauss(2048) + >>> plt.plot(np.real(y_nonstat), np.imag(y_nonstat)) + >>> plt.xlabel("Real part") + >>> plt.ylabel("Imaginary part") + >>> plt.show() + + .. plot:: misc_plots/nonstationary_phase_plot.py + +Notice that the second plot has a lot of rough edges and sharp angles. This +means that when the signal vector rotates through the phase space, it will have +to make sharp jumps in order to trace the signal. Moreover, by the definition +of instantaneous frequency, when the instantaneous phase is not differentiable, +the instantaneous frequency will be indeterminate. By contrast, the phase plot +for the stationary signal is a lot smoother, and we can expect that the +instantaneous frequency will be finite at all times. Physically, this means +that in the nonstationary signal, the variation in frequency components has no +structure, and these components can change arbitrarily. + +This phenomenon of arbitrary, unstructured changes in frequency over time is a +symptom of nonstationarity, and will become increasingly relevant as we +proceed. + +.. [1] http://www.amazon.com/Time-Frequency-Analysis-Theory-Applications/dp/0135945321 diff --git a/doc/quickstart/intro_examples_1.rst b/doc/quickstart/intro_examples_1.rst index caf78e5..2e3f367 100644 --- a/doc/quickstart/intro_examples_1.rst +++ b/doc/quickstart/intro_examples_1.rst @@ -34,7 +34,7 @@ If we now consider the energy spectrum of the signal ``z`` by squaring the modul >>> plt.grid() >>> plt.show() - .. plot:: auto_examples/plot_1_3_1_chirp_spectrum.py + .. plot:: _gallery/plot_1_3_1_chirp_spectrum.py we still can not say, from this plot, anything about the evolution in time of the frequency content. This is due to the fact that the Fourier transform is a decomposition on complex exponentials, which are of infinite duration and @@ -50,7 +50,7 @@ Wigner-Ville distribution on this signal. >>> wvd.run() >>> wvd.plot(kind='contour', extent=[0, n_points, fmin, fmax]) - .. plot:: auto_examples/plot_1_3_1_chirp_wv.py + .. plot:: _gallery/plot_1_3_1_chirp_wv.py we can see that the linear progression of the frequency with time, from 0 to 0.5, is clearly shown. @@ -66,7 +66,7 @@ If we now add some complex white gaussian noise on this signal, >>> plt.grid() >>> plt.show() - .. plot:: auto_examples/plot_1_3_1_noisy_chirp.py + .. plot:: _gallery/plot_1_3_1_noisy_chirp.py and consider the spectrum of it, @@ -79,7 +79,7 @@ and consider the spectrum of it, >>> plt.grid() >>> plt.show() - .. plot:: auto_examples/plot_1_3_1_noisy_chirp_spectrum.py + .. plot:: _gallery/plot_1_3_1_noisy_chirp_spectrum.py it is worse than before to interpret these plots. On the other hand, the Wigner-Ville distribution still show quite clearly the linear progression of the frequency with time. @@ -88,4 +88,4 @@ clearly the linear progression of the frequency with time. >>> wvd.run() >>> wvd.plot(kind='contour') - .. plot:: auto_examples/plot_1_3_1_noisy_chirp_wv.py \ No newline at end of file + .. plot:: _gallery/plot_1_3_1_noisy_chirp_wv.py diff --git a/doc/quickstart/intro_examples_2.rst b/doc/quickstart/intro_examples_2.rst index d0d1981..695ad75 100644 --- a/doc/quickstart/intro_examples_2.rst +++ b/doc/quickstart/intro_examples_2.rst @@ -35,7 +35,7 @@ follows: >>> >>> plt.show() - .. plot:: auto_examples/plot_1_3_3_transient.py + .. plot:: _gallery/plot_1_3_3_transient.py From these representations, it is difficult to localize precisely the signal in the time-domain as well as in the frequency domain. Now let us have a look at the spectrogram of this signal: @@ -47,6 +47,6 @@ frequency domain. Now let us have a look at the spectrogram of this signal: >>> spec.run() >>> spec.plot(kind="contour", threshold=0.1, show_tf=False) - .. plot:: auto_examples/plot_1_3_3_transient_spectrogram.py + .. plot:: _gallery/plot_1_3_3_transient_spectrogram.py -the transient signal appears distinctly around the normalized frequency 0.25, and between time points 125 and 160. \ No newline at end of file +the transient signal appears distinctly around the normalized frequency 0.25, and between time points 125 and 160.