signal_filters package

Submodules

signal_filters.filters module

Collection of digital filters and filter front ends for filters defined in the scipy package

In the filtering module of signal_processing tool box three base filters have been implemented:

  1. Ideal block filter (bandpass_block_filter)

  2. Kaiser filter (kaiser_bandpass_filter)

  3. Butterworth filter (butter_bandpass_filter).

Notes

  • All filters need at least a low or a high cut-off frequency as an input argument.

  • In case only a low cut-off frequency is passed, the filters behave as a high-pass filter, as only the frequencies below f_cut_low are suppress

  • In case only a high cut-off frequency is passed, the filters behave as a low pass-filter, as only the frequencies above the f_cut_high are suppressed.

  • In case both the low and high cut-off frequencies are given, the filters are band-pass filters.

  • The cut-off frequencies and the sample frequency are given in Hz.

  • The kaiser_bandpass_filter and butter_bandpass_filter are in fact front-ends to the scipy filters. To use the scipy filter, quite some steps are required to calculate the filter coefficients. The implementations in this signal_processing package take care of that and then

    call the scipy filters.

  • The Ideal block filter bandpass_block_filter is a literal copy of the Matlab band_pass2 filter

  • On top of the base filters an extra front end is defined, filter_signal, which receives the same input arguments as all three filtered mentioned above, plus an argument filter_type. With this function, we can pick which filter we want to use: kaiser, block or butterworth. The default choice is block.

signal_filters.filters.band_pass_block(omega, fs=1, lowcut=None, highcut=None)[source]

A perfect block filter with unity response

Parameters:

  • omega (class:numpy.ndarray) – Frequencies in rad/s

  • fs (float) – Sample frequency in Hz

  • lowcut (float) – Cut-off low frequency in Hz

  • highcut (float) – Cut-off high frequency in Hz

Returns:

1-D array with 1 in between lowcut, highcut and 0 elsewhere

Return type:

numpy.ndarray

Examples

Create a array with some radial frequencies running from 0 to 2.5 rad/s

>>> frequencies = np.linspace(0, 2.5, 10)

Calculate the block response

>>> band_pass_block(omega=frequencies, fs=1.0, lowcut=0.1, highcut=0.2)
array([ 0.,  0.,  0.,  1.,  1.,  0.,  0.,  0.,  0.,  0.])
signal_filters.filters.bandpass_block_filter(T, X, wfiltlo=None, wfiltup=None)[source]

Ideal band pass filter

Parameters:

  • T (arry_like) – time base ‘T’ (seconds !!) of size N OR delta time already.

  • X (array_like) – signal ‘X’ of size N. Note that I took out the M dimensional signals

  • wfiltlo (float, optional) – Filter frequency ‘wfiltlo’: lower boundary (high pass filter) in rad/s. In case wfillo is not given, zero is imposed, i.e. the band pass filter turns into a low pass filter. The high cut-off frequency must be given in that case

  • wfiltup (float, optional) – filter frequency ‘wfiltup’: upper boundary (low pass filter) in rad/s. In case wfilup is not given, the Nyquist frequency is imposed , i.e. the band pass filter turns into a high pass filter. The low cut-off frequency must be given in that case

Returns:

1-D array with the band passed filtered signal

Return type:

ndarray

Raises:

ValueError: –

  • In case no cut-off frequencies are defined * In case that wfiltlo >= wfiltup

Examples

Define as signal with some noise

>>> x = np.linspace(0, 10, num=500, endpoint=True)
>>> fs = 1 / (x[1] - x[0])
>>> y_orig = np.sin(2 * pi * x / 2.5)
>>> y_noise = y_orig + np.random.random_sample(x.size)

Filter the signal with the ideal block filter to recover the peak period at 2.5 s. Set the low and high cut-off frequencies at 3 s and 1 s, respectively

>>> f_low = 2 * pi / 3.0  # T =  3.0 s. f_low is the frequency in rad/s
>>> f_hig = 2 * pi / 1.0  # T =  1.0 s. f_hig is the frequency in rad/s
>>> y_recov = bandpass_block_filter(x, y_noise, wfiltlo=f_low, wfiltup=f_hig)

The recovered signal y_recov now contains the filter signal.

Notes

  • Explicity extend the input signal with one value in case we have an odd signal. This is to avoid a bug in the fft making the function extremely slow for an odd number of points.

  • For the rest the original matlab code is followed as close as possible, which means that the variable naming does not follow the PIP standards

signal_filters.filters.butter_bandpass_coefficients(f_cut_low, f_cut_high, fs, order=2)[source]

Return the filter coefficients of a Butterworth bandpass filter

Parameters:

  • f_cut_low – Frequency low in Hz: suppress all frequencies below this value

  • f_cut_high – Frequency high in Hz: suppress all frequencies above this value

  • fs – the sample frequencies in Hz

  • order – order of the filter: higher order is sharper drop off edge at the cutoff frequencies (order n has 20n dB/decaded drop off) (Default value = 2)

Returns:

(b, a) with the band pass filter coefficients

Return type:

tuple

Notes

  • This function is a front-end to the butter function from the scipy.signal module.

  • The butter function requires the low and/or high cut-off frequency to be expressed as a fraction of the Nyquist frequency defined as 0.5 * f_sample

  • This function allows to get the Butterworth filter coefficients using the absolute high and low cut frequency as an input

  • This function is used from the butterworth_filter function, which is a front end to the digital Butterworth filter

signal_filters.filters.butter_highpass_coefficients(f_cut_low, fs, order=2)[source]

Return the filter coefficients of a Butterworth high-pass filter

Parameters:

  • f_cut_low – suppress all frequencies below this value

  • fs – the sample frequencies

  • order – order of the filter: higher order is sharper drop off edge at the cutoff frequencies (order n has 20n dB/decade drop off) (Default value = 2)

Returns:

(b, a) with the high pass filter coefficients

Return type:

tuple

Notes

  • The high-pass filter requires a low cut-off frequency: all frequencies below f_cut_low will be suppressed

signal_filters.filters.butter_lowpass_coefficients(f_cut_high, fs, order=2)[source]

Return the filter coefficients of a Butterworth low-pass filter

Parameters:

  • f_cut_high – suppress all frequencies above this value

  • fs – the sample frequencies

  • order – order of the filter: higher order is sharper drop off edge at the cutoff frequencies (order n has 20n dB/decade drop off) (Default value = 2)

Returns:

(b, a) with the low pass filter coefficients

Return type:

tuple

Notes

  • The low-pass filter requires a high cut-off frequency: all frequencies above f_cut_high will be suppressed

signal_filters.filters.butterworth_filter(data, f_cut_low=None, f_cut_high=None, fs=1, order=2, use_filtfilt=True)[source]

Apply the Butterworth filter (high-, low, or band-pass) to the signal ‘data’.

Parameters:

  • data (ndarray) – Array with the signal to be filtered

  • f_cut_low (float or None, optional) – low frequency below which to suppress all signals (Default value = None)

  • f_cut_high (float or None, optional) – high frequency above which to suppress all signals (Default value =None)

  • fs (float, optional) – sample frequency (Default value = 1)

  • order (int) – order of the filter (Default value = 2)

  • use_filtfilt (bool, optional) – If true use the filtfilt filter from scipy to take out the phase shift. Default = True. If False, the lfilter is used instead, causing a phase shift.

Returns:

Filtered signal

Return type:

ndarray

Raises:

ValueError: –

  • In case no cut-off frequencies are defined * In case that wfiltlo >= wfiltup

Examples

First create a noise signal with 2 harmonics with a period of 2.5 s and 1.0 s

>>> x = np.linspace(0, 10, num=500)
>>> fs = 1 / (x[1] - x[0])
>>> y_orig = np.sin(2 * pi * x / 2.5) + np.sin(2 * pi * x / 1.0)
>>> y_noise = y_orig + np.random.random_sample(x.size)

Filter the signal with the high pass butt block filter to recover the peak with a period of 1.0 s. Note that for the high-pass filter, we have to define the low-cut frequency f_low (as all the frequencies below f_low are suppressed)

>>> f_low = 1 / 2.0 # T = 2.0 s. f_low is the frequency in Hz
>>> y_recover = butterworth_filter(y_noise, f_cut_low=f_low, order=4, fs=fs)

The recovered signal y_recov now contains the filtered signal.

We can also apply a low-pass filter to remove the peak with a frequency of 1 Hz

>>> f_high = 1 / 2.0  # T = 2.0 s. f_low is the frequency in Hz
>>> y_recover = butterworth_filter(y_noise, f_cut_high=f_high, order=4, fs=fs)

Perhaps better would be to use a band pass filter by both defining a high and low cut-off frequency. If we want to recover the peak with a period of 1.0 seconds we can do

>>> f_high = 1 / 0.5  # T = 0.5 s. Supress all frequencies above 2 Hz
>>> f_low = 1 / 2.0  # T = 2.0 s. Suppress all frequency below 0.5 Hz
>>> y_recover = butterworth_filter(y_noise, f_cut_high=f_high, f_cut_low=f_low,
...                                fs=fs)

Notes

  • The order argument is used to control the sharpness of the filter a high order will result in a stronger suppression of frequencies outside the cut-off limits.

  • An order which is too high also leads to an overflow of the filter. Use trial-and-error to find an optimal order. Typically, order should be between 2 and 5

  • In case the filter type is set to high, only the low cut-off frequence f_cut_low needs to be defined: all the frequencies below f_cut_low will be suppressed (hence: high-pass)

  • In case the filter type is set to low, only the high cut-off frequence f_cut_high needs to be defined: all the frequencies above f_cut_high will be suppressed (hence: low-pass)

signal_filters.filters.filter_signal(signal, filter_type='block', f_cut_low=None, f_cut_high=None, f_sampling=1, f_width_edge=0.05, ripple_db=100, order=2, constant_value=None, use_filtfilt=True)[source]

Filter the signal of input signal base on the input argument

Parameters:

  • signal (ndarray) – Input 1D signal

  • filter_type ({"mean", "butterworth", "kaiser", "block", "none"}) – Type of the filter to use. Default = block.

  • f_cut_low (Quantity or float or None) – Low cut-off frequency in Hz. Default = None, which means that f_cut_high must be given and that we are dealing with a low-pass filter

  • f_cut_high (Quantity or float or None) – High cut-off frequency in Hz. Default = None, which means that f_cut_low must be given and that we are dealing with a high-pass filter

  • f_sampling (Quantity or float) – Sampling frequency. Default = 1 Hz

  • f_width_edge (Quantity or float) – Sharpness of the stop band of the Kaiser filter. Only used when filter_type == “kaiser”. A smaller f_width_edge approaches the ideal block filter better. Default = 0.05 Hz

  • ripple_db (float, optional) – The desired attenuation in the stop band in dB. Only used when filter_type == “kaiser”. Default = 100 dB

  • order (int) – Order of the Butterworth filter. Only used when filter_type==”butterworth”. Default = 2. Higher values causes to approach the ideal block band filter, but also a more unstable filter.

  • constant_value (float) – Constant to use for the Kaiser filter when a phase shift correction is applied. Default = None, which means that the phase shift is not correct. Since we use filtfilt by default, this is not required. Only used when filter_type == “kaiser”

  • use_filtfilt (bool) – If true use the filtfilt filter from scipy to take out the phase shift. Default = True.

Returns:

Filtered signal

Return type:

ndarray

Raises:

ValueError: –

  • In case no cut-off frequencies are defined * In case that wfiltlo >= wfiltup

Notes

  • This function is a front end to the other filters in this module

  • The Nyquist frequency follows from the sample frequency as f_nyq = f_s / 2

  • At least one of the cut-off frequencies must be properly defined, i.e. in the range 0 ~ f_nyq

  • For f_cut_low == None (or 0), the filter is acting as a low-pass filter with its cut-off frequency given by f_cut_high. All frequencies above f_cut_high will be suppressed

  • For f_cut_high == None or f_cut_high >= f_nyquist, the filter is acting as a high-pass filter with its cut-off frequency given by f_cut_low. All frequencies below f_cut_low will be suppressed.

  • If both f_cut_low and f_cut_high are in the range 0 ~ f_nyquist, the filter is acting as a band pass filter

signal_filters.filters.kaiser_bandpass_coefficients(lowcut, highcut, fs, f_width_edge=None, ripple_db=200)[source]

Calculates the Kaiser band pass filter coefficients

Parameters:

  • lowcut (float) – The low cut-off frequency [Hz]

  • highcut (float) – The high cut-off frequency [Hz]

  • fs (float) – Sample frequency [Hz]

  • f_width_edge (float, optional) – Give the distance in Hz for the edges of the filter. If None, take 1 Hz. Default = None

  • ripple_db (float, option.) – The desired attenuation in the stop band, in dB. Default = 200 db

Returns:

(band_pass_coefficients, n_delay)

  • band_pass_coefficients: list with the band pass coefficients

  • n_delay: the number of delay samples

Return type:

tuple

signal_filters.filters.kaiser_bandpass_filter(data, lowcut=None, highcut=None, fs=1, f_width_edge=0.01, ripple_db=100.0, cval=None, replace_cval_with_unfiltered_signal=False, replace_cval_with_nan=False, use_filtfilt=True)[source]

Filter the data with the FIR kaiser

Parameters:

  • data (ndarray) – Input signal 1D

  • lowcut (float or None, optional) – Lower frequency in Hz. If the low-cut frequency is not given, the filter turns into a low-pass filter with only a high cut-off frequency defined. The high-cut must be defined in that case

  • highcut (float or None, optional) – higher frequency in Hz. If the high-cut frequency is not given the filter turns into a high-pass filter with only the low cut-off frequency defined. The low-cut frequency must be defined in that case

  • fs (float, optional) – Sample frequency Hz. Default = 1.0 Hz

  • f_width_edge (float, optional) – width of the edges. Default = 0.01

  • ripple_db (float, optional) – Suppression of the stop band in dB. Default = 100

  • cval (float, optional) – Constant value just for the shift. Default = 0

  • replace_cval_with_unfiltered_signal (bool, optional) – You may choose to replace the zero introduced due to the shift. Default = False

  • replace_cval_with_nan (bool, optional) – Replace the cval introduced by the shift correction with nans

  • use_filtfilt (bool, optional) – If true use the filtfilt which takes out the phase shift. Default = True

Returns:

Filtered signal

Return type:

ndarray

Raises:

ValueError: –

  • In case no cut off frequency is defined * In case the low cut-off frequency is higher than or equal to the high cut-off frequency

Notes

  • If both the low and high cut-off frequencies lowcut and highcut are defined, the kaiser filter is a band pass filter

  • In case only the low cut-off frequency lowcut is defined, all frequencies higher than lowcut will pass (i.e. the filter is a high-pass filter)

  • In case only the high cut-off frequency highcut is defined, all frequencies lower than highcut will pass (i.e. the filter is a low-pass filter)

signal_filters.filters.matrix_fft(om_spec=None, spec=None, eps=None, om_DFT=None, DFT=None, om_RAO=None, RAO=None, t=None, TT=None, s=4)[source]

Matrix FFT routine obtained from matlab

Parameters:

  • om_spec (array_like) – frequencies of spectrum

  • spec (array_like) – spectrum density

  • eps (float) – eps must always be in degrees

  • om_DFT (array_like) – frequencies of DFT

  • DFT (int) – DFT components

  • om_RAO (array_like) – array with frequencies

  • RAO (float) – rao array

  • t (array_like) – array with time values

  • TT (array_like) – array with values

  • s (int) – Switch to determine the mode of operation, must be integer between 1 and 5

Returns:

List with arrays. The amount of arrays returned depends on the mode used

Return type:

class:list

Notes

  • FFT routine copied from HMS matfft by RvD. Modified input argument treatment: named input arguments are used Major difference is that this routine is only for 1-D signals, in the matlab code the array may contain more than one 1-D signals.

  • In case more than one signal needs to be filtered, in the python version you need to loop over the signals. Since the filter is linear it does not influence the calculation time

  • Depending on the option s, the following input argument are used an may not be None:

    1.   spec 2 DFT    use: [om_DFT,DFT]       = matfft(om_spec,spec,eps,1)
2.   apply RAO     use: [om_DFT,DFT]       = matfft(om_DFT, DFT,om_RAO,
                                                                RAO,2)
3.   DFT 2 TT      use: [t     ,TT]        = matfft(om_DFT, DFT,3)
4.   TT 2 DFT      use: [om_DFT,DFT]       = matfft(t,TT,4)
5.   DFT 2 spec    use: [om_spec,spec,eps] = matfft(om_DFT, DFT,5)
     eps always in degrees!!!

  • The original code was followed as close as possible which means the variable naming does not follow the PIP standards

signal_filters.filters.remove_phase_shift(signal, degrees=True)[source]

Apply a phase shift to the angles of signal to get rid of phase-jumps at zero degrees

Parameters:

  • signal – the input array with the yaw angle in degrees or rad

  • degrees – flag to indicate that the yaw is in degrees. Default is True

Returns:

array with shifted phase such that there are no phase-jumps at zero degrees

Return type:

ndarray

Notes

  • The yaw angle sometimes has discontinuities when the phase angle goes from 0 to 360. To avoid this, in this routine the range of the yaw is put as good as possible in the range between 0 and 360.

  • In case the MRU 6-DOF motions need to be filtered in order to obtain the accelerations, the yaw needs to be put into the range 0~360, otherwise the discontinuity at 0 degrees will lead to an artificial spike in the acceleration.

  • In case the signal passes the zero degree two times, this routine will fail.

Examples

Lets assume we have a yaw signal which passed the 0 degrees. If we want to get the derivative of the yaw, we need to put it into a range 0~360 such that we get rid of the discontinuity at 0.

>>> yaw = np.hstack([np.linspace(350, 360, num=5, endpoint=False),
...                  np.linspace(0, 10, num=5)])
>>> print(yaw)
[ 350.   352.   354.   356.   358.     0.     2.5    5.     7.5   10. ]

In case we want to calculate the first derivative for the yaw, we get a spike at 0

>>> dyawdt = np.diff(yaw)
>>> print(dyawdt)
[   2.     2.     2.     2.  -358.     2.5    2.5    2.5    2.5]

In order to get rid of this spike, first remove the phase shift

>>> yaw2 = remove_phase_shift(yaw)
>>> print(yaw2)
[  0.    2.    4.    6.    8.   10.   12.5  15.   17.5  20. ]

Now the derivative yields proper values with the spike of -358 degrees

>>> dyaw2dt = np.diff(yaw2)
>>> print(dyaw2dt)
[ 2.   2.   2.   2.   2.   2.5  2.5  2.5  2.5]

signal_filters.utils module

Some utilities to support signal processing.

class signal_filters.utils.SignalComponent(key, type='cos', amplitude=1.0, frequency=1.0, phase_shift=0.0)[source]

Bases: object

Class to store the properties and data a signal signal component

Parameters:

  • key (str) – Name of the signal component

  • type ({"sin", "cos", "random_normal"}) – Type of the the signal. Can be ‘cos’, ‘sin’ and ‘random_normal’ for now. Default = “cos”

  • amplitude (float, optional) – Amplitude of the signal. Default = 1.0

  • frequency (float, optional) – Frequency of the periodic signals. Default = 1 Hz

  • phase_shift (float, optional) – Phase shift of the periodic components. Default = 0

report()[source]

Create a report of the current signal component properties

class signal_filters.utils.SignalGenerator(time_length=100, sample_frequency=1.0, create_data_frame=True)[source]

Bases: object

Create a signal and generate the data consisting of pure sin components and noise

Parameters:

  • time_length (float, optional) – Length of the signal in seconds. Default = 100 s

  • sample_frequency (float, optional) – Sample frequency in Hz. Default = 1 Hz

  • create_data_frame (bool, optional) – if true, also create a pandas data frame to keep all the components

Notes

  • A simple signal generator to generate a signal with given components

Examples

A signal consisting of 2 harmonic components and some noise can be created as follows

>>> logger = create_logger(console_log_format_clean=True)
>>> np.random.seed(0)  # makes sure to generate the same random numbers
>>> signal = SignalGenerator(time_length=5, sample_frequency=10)

The SignalGenerator object has been initialized but no signal components have been defined yet. We can add components by using the add_components method.

>>> signal.add_component(key="f_cos", amplitude=1.5, period=2.5, phase_shift=0.0)

The add_component method can be called several time to add more components. To add a sin component with phase shift 0.25 pi we can do

>>> signal.add_component(key="f_sin", signal_type="sin", amplitude=0.3, period=0.5,
...                     phase_shift=0.25 * pi)

And to add some extra noise do

>>> signal.add_component(key="noise", signal_type="random_normal", amplitude=0.2)

To invest which components have been adedd to the report_components method:

>>> signal.report_components()
# component: f_cos      --------------------------------------------------
Key                  : f_cos
Type                 : cos
Amplitude            : 1.5
Frequency            : 0.4
Phase                : 0
# component: f_sin      --------------------------------------------------
Key                  : f_sin
Type                 : sin
Amplitude            : 0.3
Frequency            : 2
Phase                : 0.79
# component: noise      --------------------------------------------------
Key                  : noise
Type                 : random_normal
Amplitude            : 0.2
Frequency            : None
Phase                : 0

Adding components does not generate the signal yet. For this use the generate method

>>> signal.generate()

After generate has been called the data given by the signal components has been generated and stored in the data attribute.

>>> print(signal.data)
[ 2.0649425   1.80020814  1.46327728  1.24532508  1.04105464  0.48020195
  0.55150542 -0.35827375 -0.95561921 -1.01021343 -0.97258474 -0.83650807
 -1.38289485 -1.76014355 -1.44208923 -0.93465859 -0.39001821 -0.72663093
 -0.51476493 -0.21283052  0.16505956  1.20176587  1.21940984  0.86972051
  1.77062852  1.4212589   1.7293284   1.23009291  1.10370228  0.9614148
  0.70664701  0.43712024 -0.50555946 -1.33113473 -1.16191556 -0.97012366
 -0.88130464 -1.29462642 -1.86194392 -1.59132243 -1.21110405 -0.97283761
 -1.02685331 -0.1872234  -0.14394181  0.58804267  0.82048308  1.20202067
  0.69537395  1.27412954  1.53303872]

The data attribute is jut a numpy array which contains the summation of all individual components.

The data is by default also stored as a pandas DataFrame in the data_frame attribute.

>>> from tabulate import tabulate
>>> print(tabulate(signal.data_frame.head(5), headers="keys", tablefmt="psql"))
+------------+---------+---------+------------+-----------+
|   Time [s] |   Total |   f_cos |      f_sin |     noise |
|------------+---------+---------+------------+-----------|
|        0   | 2.06494 | 1.5     |  0.212132  | 0.35281   |
|        0.1 | 1.80021 | 1.45287 |  0.267302  | 0.0800314 |
|        0.2 | 1.46328 | 1.31446 | -0.0469303 | 0.195748  |
|        0.3 | 1.24533 | 1.09345 | -0.296307  | 0.448179  |
|        0.4 | 1.04105 | 0.80374 | -0.136197  | 0.373512  |
+------------+---------+---------+------------+-----------+

The Total column of the data frame contains the same total signal as we have seen in the data attribute. However, all the sub component can be individually accessed in the data frame by the columns corresponding to the key name we have given to the components.

add_component(key, signal_type='cos', amplitude=1.0, frequency=None, period=None, phase_shift=0.0)[source]

Add a new signal component to the list of signals

Parameters:

  • key (str,) – Name of this component

  • signal_type ({"sin", "cos", "random_signal}) – Type of the signal. Default = cos

  • amplitude (float, optional) – Amplitude of the signal

  • frequency (float, optional) – Frequency of the current signal. Default = None

  • period (float, optional) – Period of the current signal. Only frequency -or- period can be given. Default = None

  • phase_shift (float, optional) – Phase shift of the current signal

generate()[source]

Generate the signal data based on the components added to this signal

generate_data_frame(index_name='Time [s]')[source]

Turn the data into a pandas data frame. Need to create the data first

Parameters:

index_name (str, optional) – The name given to the index. Default = “Time [s]”

info()[source]

Create some output with information of the signal. Just use pandas for that

report_components()[source]

Make a report of all the signal components added to the SignalGenerator

signal_filters.utils.get_peaks(freq, psd, thres=0.01, min_dist=3, max_number_of_peaks=3, sort_peaks=True)[source]

Get the peaks in the power spectral density ‘psd’

Parameters:

  • freq (ndarray) – frequencies belong to psd

  • psd (ndarray) – array with power spectral densities

  • thres (float between [0., 1.]) – Normalized threshold. Only the peaks with amplitude higher than the threshold will be detected.

  • max_number_of_peaks (int, optional) – maximum number of peaks return. Default = 3

  • min_dist (int) – Minimum distance between each detected peak. The peak with the highest amplitude is preferred to satisfy this constraint.

  • sort_peaks (bool) – If True, sort the peaks in the order of their height. Default = True

Returns:

(f_peak, psd_value_at_peak) Two ndarrays with peak value and psd values

Return type:

tuple

Examples

First create a signal with three harmonical components

>>> time_length = 1800 # 30 minutes sample
>>> f_sample = 25 # 25 Hz
>>> time = np.linspace(0, stop=time_length, num=time_length * f_sample)
>>> amplitudes = np.array([1, 2, 0.5])  # three amplitudes of the peaks
>>> f_peaks = np.array([0.03, 0.08, 1])  # three peak frequencies in Hz
>>> omega_p = 2 * np.pi * f_peaks   # convert frequencies to rad/s
>>> signal = np.zeros(time.shape)
>>> for ii, om_p in enumerate(omega_p):
...    signal += amplitudes[ii] * np.sin(om_p * time)
>>> signal += np.random.normal(omega_p.size, scale=0.4)

The spectrum of this signal should contain three peaks. First calculate the spectrum

>>> from scipy.signal import welch
>>> freq, psd_signal = welch(signal, fs=f_sample, nperseg=4096)

See if we can find the peaks in the spectrum with the get_peaks function

>>> f_peaks_found, v_p = get_peaks(freq=freq, psd=psd_signal, sort_peaks=False)
>>> for ii, fp in enumerate(f_peaks_found):
...    print("Peak {} at {:4.2f} Hz. Found peak {:4.2f} Hz".format(
                ii, f_peaks[ii], fp))
Peak 0 at 0.03 Hz. Found peak 0.03 Hz
Peak 1 at 0.08 Hz. Found peak 0.08 Hz
Peak 2 at 1.00 Hz. Found peak 1.00 Hz

Indeed it can be seen that the peaks are correctly found in the spectrum.

We can also sort the peaks in order of their magnitude (which is the default behaviour

>>> f_peaks_found, v_p = get_peaks(freq=freq, psd=psd_signal)
>>> for ii, fp in enumerate(f_peaks_found):
...    print("Peak {} at {:4.2f} Hz with A ={:7.2f}".format(ii, fp, v_p[ii]))
Peak 0 at 0.08 Hz with A = 215.22
Peak 1 at 0.03 Hz with A =  54.12
Peak 2 at 1.00 Hz with A =  13.23

The highest peaks corresponds to the component with the largest amplitude indeed. This makes it easy to get the most dominant peak frequency as it always the first components

>>> print("Most dominant frequency in spectrum is : {:.2f} Hz".format(
            f_peaks_found[0]))
Most dominant frequency in spectrum is : 0.08 Hz

Notes

  • This function is a front end to the peak function from the peakutils module.

  • The peakutil.peak function only returns the indices of the peaks in the array, whereas the get_peaks function returns the actual frequencies belonging to the peaks and also the peak values belonging to the peaks

  • Also, the peaks are sorted in height, such that the first peak is always the highest peak

Module contents