Source code for skrf.time

.. module:: skrf.time
time (:mod:`skrf.time`)

Time domain functions

.. autosummary::
   :toctree: generated/


from __future__ import annotations

import warnings

# imports for type hinting
from typing import TYPE_CHECKING, Callable

import numpy as np
from numpy.fft import fft, fftshift, ifft, ifftshift, irfft, rfft
from scipy import signal
from scipy.ndimage import convolve1d

from .util import find_nearest_index

    from .network import Network

[docs] def indexes(y: np.ndarray, thres: float = 0.3, min_dist: int = 1) -> np.ndarray: """ Peak detection routine. Finds the numeric index of the peaks in *y* by taking its first order difference. By using *thres* and *min_dist* parameters, it is possible to reduce the number of detected peaks. *y* must be signed. Parameters ---------- y : ndarray (signed) 1D amplitude data to search for peaks. thres : float between [0., 1.], optional Normalized threshold. Only the peaks with amplitude higher than the threshold will be detected. Default is 0.3 min_dist : int, optional Minimum distance between each detected peak. The peak with the highest amplitude is preferred to satisfy this constraint. Default is 1 Returns ------- ndarray Array containing the numeric indexes of the peaks that were detected Notes ----- This function was taken from peakutils-1.1.0 """ #This function was taken from peakutils, and is covered # by the MIT license, included below: #The MIT License (MIT) #Copyright (c) 2014 Lucas Hermann Negri #Permission is hereby granted, free of charge, to any person obtaining a copy #of this software and associated documentation files (the "Software"), to deal #in the Software without restriction, including without limitation the rights #to use, copy, modify, merge, publish, distribute, sublicense, and/or sell #copies of the Software, and to permit persons to whom the Software is #furnished to do so, subject to the following conditions: #The above copyright notice and this permission notice shall be included in #all copies or substantial portions of the Software. #THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR #IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, #FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE #AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER #LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, #OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN #THE SOFTWARE. if isinstance(y, np.ndarray) and np.issubdtype(y.dtype, np.unsignedinteger): raise ValueError("y must be signed") thres = thres * (np.max(y) - np.min(y)) + np.min(y) min_dist = int(min_dist) # compute first order difference dy = np.diff(y) # propagate left and right values successively to fill all plateau pixels (0-value) zeros, = np.where(dy == 0) # check if the singal is totally flat if len(zeros) == len(y) - 1: return np.array([]) while len(zeros): # add pixels 2 by 2 to propagate left and right value onto the zero-value pixel zerosr = np.hstack([dy[1:], 0.]) zerosl = np.hstack([0., dy[:-1]]) # replace 0 with right value if non zero dy[zeros]=zerosr[zeros] zeros, = np.where(dy == 0) # replace 0 with left value if non zero dy[zeros] = zerosl[zeros] zeros, = np.where(dy == 0) # find the peaks by using the first order difference peaks = np.where((np.hstack([dy, 0.]) < 0.) & (np.hstack([0., dy]) > 0.) & (y > thres))[0] # handle multiple peaks, respecting the minimum distance if peaks.size > 1 and min_dist > 1: highest = peaks[np.argsort(y[peaks])][::-1] rem = np.ones(y.size, dtype=bool) rem[peaks] = False for peak in highest: if not rem[peak]: sl = slice(max(0, peak - min_dist), peak + min_dist + 1) rem[sl] = True rem[peak] = False peaks = np.arange(y.size)[~rem] return peaks
[docs] def find_n_peaks(x: np.ndarray, n: int, thres: float = 0.9, **kwargs) -> list[int]: """ Find a given number of peaks in a signal. Parameters ---------- x : np.ndarray signal n : int number of peaks to search for thres : float, optional threshold, default is 0.9 **kwargs : optional keyword arguments passed to :func:`indexes` Returns ------- peak_idxs : list of int List containing the numeric indexes of the peaks that were detected Raises ------ ValueError If no peaks are found. """ for _dummy in range(10): idx = indexes(x, **kwargs) if len(idx) < n: thres *= .5 else: peak_vals = sorted(x[idx], reverse=True)[:n] peak_idxs = [x.tolist().index(k) for k in peak_vals] return peak_idxs raise ValueError('Couldnt find %i peaks' % n)
time_lookup_dict = { "s": 1, "ms": 1e-3, "us": 1e-6, "µs": 1e-6, "ns": 1e-9, "ps": 1e-12 }
[docs] def detect_span(ntwk: Network, t_unit: str = "") -> float: """ Detect the correct time-span between two largest peaks. Parameters ---------- ntwk : :class:`` network to get data from t_unit : str Time unit for start, stop, center and span arguments, defaults to seconds (s). Possible values: * 's': seconds * 'ms': milliseconds * 'µs' or 'us': microseconds * 'ns': nanoseconds (default) * 'ps': picoseconds Returns ------- span : float in unit t_unit """ if t_unit == "": # Do not raise in autogate mode, where all parameters are None warnings.warn(''' Time unit not passed: uses 's' per default. ''', DeprecationWarning, stacklevel=2) t_unit = 's' x = ntwk.s_time_db.flatten() p1, p2 = find_n_peaks(x, n=2) # distance to nearest neighbor peak span = abs(ntwk.frequency.t[p1]-ntwk.frequency.t[p2]) return span / time_lookup_dict[t_unit]
[docs] def get_window(window: str | tuple | Callable, Nx: int, **kwargs) -> np.ndarray: """Calls a custom window function or `scipy.signal.get_window()` depending on the window argument. Parameters ---------- window : str, tuple or Callable The callable window function or a valid window string for `scipy.signal.get_window()`. Nx : int Number of samples Returns ------- window : np.ndarray Window samples. """ if callable(window): return window(Nx, **kwargs) else: return signal.get_window(window, Nx=Nx, **kwargs)
[docs] def time_gate(ntwk: Network, start: float = None, stop: float = None, center: float = None, span: float = None, mode: str = 'bandpass', window=('kaiser', 6), method: str ='fft', fft_window: str='cosine', conv_mode: str='wrap', t_unit: str = "") -> Network: """ Time-domain gating of one-port s-parameters with a window function from The gate can be defined with start/stop times, or by center/span. All times are in units of seconds but can be changed using the `t_unit` parameter. Common windows are: * ('kaiser', 6) * 6 # integers are interpreted as kaiser beta-values * 'hamming' * 'boxcar' # a straight up rect * Callable function which accepts integer with window size as argument If no parameters are passed this will try to auto-gate the largest peak. Parameters ---------- ntwk : :class:`` network to operate on start : number, or None start of time gate in t_unit. stop : number, or None stop of time gate in t_unit. center : number, or None center of time gate, in t_unit. If None, and span is given, the gate will be centered on the peak. span : number, or None span of time gate, in t_unit. If None span will be half of the distance to the second tallest peak mode : ['bandpass', 'bandstop'] mode of gate window : string, float, or tuple passed to `window` arg of `scipy.signal.get_window()` or callable function method : str Gating method. There are 3 option: 'convolution', 'fft', 'rfft'. With *'convolution'*, the time-domain gate gets transformed into frequency-domain using inverse FFT and the gating is then achieved by convolution with the frequency-domain data. With *'fft'* (default), the data gets transformed into time-domain using inverse FFT and the gating is achieved by multiplication with the time-domain gate. The gated time-domain signal is then transformed back into frequency-domain using FFT. As only positive signal frequencies are considered for the inverse FFT (with or without a dc component), the resulting time-domain signal has the same number of samples as in the frequency-domain, but is complex-valued. This method is also know as *time-domain band-pass mode*. With *'rfft'*, the procedure is the same as with *'fft'*, but the inverse FFT uses a complex-conjugate copy of the positive signal frequencies for the negative frequencies (Hermitian frequency response). A dc sample is also required. The resulting time-domain signal is real-valued and has twice the number of samples, which gives an improved time resolution. This method is also known as *time-domain low-pass mode*. fft_window : str or tuple or None Frequency-domain window applied before the inverse FFT in case of the (R)FFT method. This parameter takes the same values as the `window` parameter. Example: `window='hann` (default), or `window=('kaiser', 5)`, or `window=None`. The window helps to remove artefacts such as time-domain sidelobes of the pulses, but it is a trade-off with the achievable pulse width. The window is removed when the gated time-domain signals is transformed back into frequency-domain. conv_mode : str Extension mode for the convolution (if selected) determining how the frequency-domain gate is extended beyond the boundaries. This has a large effect on the generation of gating artefacts due to boundary effects. The optimal mode depends on the data. See the parameter description of `scipy.ndimage.convolve1d` for the available options. t_unit : str Time unit for start, stop, center and span arguments, defaults to seconds (s). Possible values: * 's': seconds * 'ms': milliseconds * 'µs' or 'us': microseconds * 'ns': nanoseconds (default) * 'ps': picoseconds Note ---- You cant gate things that are 'behind' strong reflections. This is due to the multiple reflections that occur. If you need to time-gate an N-port network, then you should gate each s-parameter independently. Returns ------- ntwk : Network copy of ntwk with time-gated s-parameters .. warning:: Depending on sharpness of the gate, the band edges may be inaccurate, due to properties of FFT. We do not re-normalize anything. """ if ntwk.nports >1: raise ValueError('Time-gating only works on one-ports. Try passing `ntwk.s11` or `ntwk.s21`.') if t_unit == "": if not all([e is None for e in [start, stop, center, span]]): # Do not raise in autogate mode, where all parameters are None warnings.warn(''' Time unit not passed: uses 's' per default. ''', DeprecationWarning, stacklevel=2) t_unit = 's' t_mult = time_lookup_dict[t_unit] if start is not None and stop is not None: start *= t_mult stop *= t_mult span = abs(stop-start) center = (stop+start)/2. else: if center is None: # they didnt provide center, so find the peak n = ntwk.s_time_mag.argmax() center = ntwk.frequency.t_ns[n] if span is None: span = detect_span(ntwk, t_unit='ns') center *= t_mult span *= t_mult start = center - span / 2. stop = center + span / 2. ntwk_gated = ntwk.copy() method = method.lower() n_fd = ntwk.frequency.npoints df = ntwk.frequency.step if method == 'convolution': # frequency-domain gating n_td = n_fd # create dummy-window window_fd = np.ones(n_fd) elif method == 'fft': # time-domain band-pass mode n_td = n_fd if fft_window is not None: # create band-pass window (zero on both lower and upper limit, one at center) window_fd = get_window(fft_window, n_fd) else: # create dummy-window window_fd = np.ones(n_fd) elif method == 'rfft': # time-domain low-pass mode if ntwk.f[0] > 0.0: # no dc point included warnings.warn('The network data to be gated does not contain the dc point (0 Hz). This is required for the ' 'selected low-pass gating mode. Please consider to include the dc point if the results are ' 'inaccurate, either by direct measurement of by extrapolation using ' 'skrf.Network.extrapolate_to_dc().', UserWarning, stacklevel=2) n_td = 2 * n_fd - 1 if fft_window is not None: # create low-pass window (one at lower limit at f=0, zero on upper limit) window_fd = get_window(fft_window, 2 * n_fd) window_fd = window_fd[n_fd:] else: # create dummy-window window_fd = np.ones(n_fd) else: raise ValueError(f'Invalid parameter method=`{method}`') # apply frequency-domain window ntwk_gated.s[:, 0, 0] = ntwk_gated.s[:, 0, 0] * window_fd # create time vector t = np.fft.fftshift(np.fft.fftfreq(n_td, df)) # find start/stop gate indices start_idx = find_nearest_index(t, start) stop_idx = find_nearest_index(t, stop) # create gating window window_width = abs(stop_idx - start_idx) + 1 window = get_window(window, window_width) # create the gate by padding the window with zeros gate = np.zeros_like(t) gate[start_idx:stop_idx+1] = window if method == 'convolution': # frequency-domain gating kernel = fftshift(fft(ifftshift(gate), norm='forward')) ntwk_gated.s[:, 0, 0] = convolve1d(ntwk_gated.s[:, 0, 0], kernel, mode=conv_mode) elif method == 'fft': # time-domain band-pass mode s_td = fftshift(ifft(ntwk_gated.s[:, 0, 0])) s_td_g = s_td * gate ntwk_gated.s[:, 0, 0] = fft(ifftshift(s_td_g)) elif method == 'rfft': # time-domain low-pass mode s_td = fftshift(irfft(ntwk_gated.s[:, 0, 0], n=len(t))) s_td_g = s_td * gate ntwk_gated.s[:, 0, 0] = rfft(ifftshift(s_td_g)) # remove frequency-domain window ntwk_gated.s[:, 0, 0] = ntwk_gated.s[:, 0, 0] / window_fd if mode == 'bandstop': ntwk_gated = ntwk - ntwk_gated elif mode == 'bandpass': pass else: raise ValueError('mode should be \'bandpass\' or \'bandstop\'') return ntwk_gated