Shulou(Shulou.com)06/01 Report--

This article mainly explains "how to implement EMD algorithm in Python". Interested friends may wish to have a look. The method introduced in this paper is simple, fast and practical. Let's let the editor take you to learn "how Python implements the EMD algorithm".

SSVEP signal contains spontaneous EEG and a large number of external interference signals, which is a typical nonlinear and non-stationary signal. Traditional filtering methods usually do not meet the conditions for nonlinear and non-stationary analysis. In 1998, Huang E proposed Hilbert-Huang transform (HHT), which includes empirical mode decomposition (EMD) and Hilbert transform (HT). EMD can decompose the original signal into a series of intrinsic mode functions (IMF) [1]. The IMF component is an oscillation function with time-varying frequency, which can reflect the local characteristics of non-stationary signals, so it is more appropriate to use it to decompose nonlinear and non-stationary SSVEP signals.

Netizen Aeo [2] provides the following algorithm process analysis.

Algorithm process analysis

Screening (Sifting) to find the extreme points all the maximum and minimum fitting envelope curves of the signal sequence through the Find Peaks algorithm through the maximum and minimum groups of the signal sequence, two smooth peak / trough fitting curves are obtained by cubic spline interpolation. That is, the mean envelope of the upper envelope and the lower envelope of the signal averages the two extreme curves to obtain the average envelope of the original signal of the intermediate signal minus the mean envelope, and the intermediate signal to judge the eigenmode function (IMF) IMF needs to meet two conditions: 1) in the whole data segment, the number of extreme points and zero-crossing points must be equal or the difference can not exceed one at most. 2) at any time, the average value of the upper envelope formed by the local maximum point and the lower envelope formed by the local minimum point is zero, that is, the upper and lower envelopes are locally symmetrical relative to the time axis. The first intermediate signal obtained by IMF 1 which satisfies the IMF condition is the first eigenmode function component IMF 1 of the original signal (the new data after subtracting the envelope average from the original data. If there are negative local Maxima and positive local minima, it shows that it is not an eigenmode function and needs to be screened. ) after the first IMF is obtained by using the above method, the original signal minus IMF1 is used as the new original signal, and then through the above screening analysis, the IMF2 can be obtained, and so on, the EMD decomposition can be completed.

The above analysis process is explained by the formula below.

EMD algorithm steps

Any complex signal can be regarded as the sum of many different intrinsic modal functions, any modal function can be linear or nonlinear, and any two modes are independent of each other. Based on this assumption, the EMD decomposition steps of complex signals are as follows:

Step 1:

Find all the extreme points of the signal, connect the local maximum points into the upper envelope through the cubic spline curve, and connect the local minimum points into the lower envelope. The upper and lower envelopes contain all the data points.

Step 2:

From the average value of the upper envelope and the lower envelope, it is concluded that

If the condition of IMF is satisfied, it can be considered as the first IMF component of.

Step 3:

If it does not meet the IMF condition, it will be used as the original data and repeat steps 1 and 2 to get the average value of the upper and lower envelopes. By calculating whether it is suitable for the necessary conditions of the IMF component, if it does not meet, repeat as the previous two steps until the premise is met. The first IMF is expressed as follows:

Step 4:

Will be separated from the signal:

The above three steps will be repeated as the original signal, and the second IMF component will be obtained until the first IMF component, and the following will be obtained:

Step 5:

When it becomes a monotone function, the rest becomes a residual component. The sum of all IMF components and residual components is the original signal:

Case 1---Python implements EMD case

Combined with the above algorithm analysis process, look at the algorithm from a code point of view.

1. Find the maximum point and minimum point

From scipy.signal import argrelextrema

"" obtain the extreme point of the signal sequence through the argrelextrema function of Scipy "" # construct 100 random numbers data = np.random.random (100) # get the maximum max_peaks = argrelextrema (data, np.greater) # get the minimum min_peaks = argrelextrema (data, np.less)

# draw the extreme point image plt.figure (figsize = (18pj6)) plt.plot (data) plt.scatter (max_peaks, data [max _ peaks], cymbals, label='Max Peaks') plt.scatter (min_peaks, data [min _ peaks], cymbals, label='Max Peaks') plt.legend () plt.xlabel ('time (s)') plt.ylabel ('Amplitude') plt.title ("Find Peaks")

two。 Fitting envelope function

This step is the core step of EMD and the premise of decomposing the eigenmode function IMFs.

From scipy.signal import argrelextrema

# perform spline difference import scipy.interpolate as spi

Data = np.random.random (100)-0.5index = list (range (len (data)

# get the extreme point max_peaks = list (argrelextrema (data, np.greater) [0]) min_peaks = list (argrelextrema (data, np.less) [0])

# fit the extreme points into the curve ipo3_max = spi.splrep (max_peaks, data [max _ peaks], kappa 3) # import the sample points and generate the parameter iy3_max = spi.splev (index, ipo3_max) # generate interpolation according to the observation points and spline parameters

Ipo3_min = spi.splrep (min_peaks, data [min _ peaks], kappa 3) # sample points are imported to generate parameters iy3_min = spi.splev (index, ipo3_min) # interpolation is generated based on observation points and spline parameters

# calculate the average envelope iy3_mean = (iy3_max+iy3_min) / 2

# drawing image plt.figure (figsize = (18Power6)) plt.plot (data, label='Original') plt.plot (iy3_max, label='Maximun Peaks') plt.plot (iy3_min, label='Minimun Peaks') plt.plot (iy3_mean, label='Mean') plt.legend () plt.xlabel ('time (s)') plt.ylabel ('microvolts (uV)') plt.title ("Cubic Spline Interpolation")

The new signal is obtained by subtracting the average envelope from the original signal. if there are negative local Maxima and positive local minima in the new signal, it shows that it is not an eigenmode function and needs to be screened.

Get the eigenmode function (IMF) def sifting (data): index = list (range (len (data)

Max_peaks = list (argrelextrema (data, np.greater) [0]) min_peaks = list (argrelextrema (data, np.less) [0])

Ipo3_max = spi.splrep (max_peaks, data [max _ peaks], ipo3_max 3) # sample points are imported, generating parameters iy3_max = spi.splev (index, ipo3_max) # generate interpolation based on observation points and spline parameters

Ipo3_min = spi.splrep (min_peaks, data [min _ peaks], kappa 3) # sample points are imported to generate parameters iy3_min = spi.splev (index, ipo3_min) # interpolation is generated based on observation points and spline parameters

Iy3_mean = (iy3_max+iy3_min) / 2 return data-iy3_mean

Def hasPeaks (data): max_peaks = list (argrelextrema (data, np.greater) [0]) min_peaks = list (argrelextrema (data, np.less) [0])

If len (max_peaks) > 3 and len (min_peaks) > 3: return True else: return False

# judge IMFsdef isIMFs (data): max_peaks = list (argrelextrema (data, np.greater) [0]) min_peaks = list (argrelextrema (data, np.less) [0])

If min (data [max _ peaks])

0: return False else: return True

Def getIMFs (data): while (not isIMFs (data)): data = sifting (data) return data

# EMD function def EMD (data): IMFs = [] while hasPeaks (data): data_imf = getIMFs (data) data = data-data_imf IMFs.append (data_imf) return IMFs

# draw a comparison diagram data = np.random.random (1000)-0.5IMFs = EMD (data) n = len (IMFs) + 1

# original signal plt.figure (figsize = (18Power15)) plt.subplot (n, 1,1) plt.plot (data, label='Origin') plt.title ("Origin")

# several IMFs curves for i in range (0Len (IMFs)): plt.subplot (n, 1, iTun2) plt.plot (IMFS [I]) plt.ylabel ('Amplitude') plt.title ("IMFs" + str (iTun1))

Plt.legend () plt.xlabel ('time (s)') plt.ylabel ('Amplitude')

Case 2 Murray-using PyEMD tools to implement EMD

# Import tool library import numpy as npfrom PyEMD import EMD, Visualisation

Build signal

Time t: 0 to 1s, sampling frequency 100Hz, S is a composite signal.

# Construction signal t = np.arange (0jue 1,0.01) S = 2*np.sin (2*np.pi*15*t) + 4*np.sin (2*np.pi*10*t) * np.sin (2*np.pi*t*0.1) + np.sin (2*np.pi*5*t) # extract imfs and remaining emd = EMD () emd.emd (S) imfs Res = emd.get_imfs_and_residue () # draw IMFvis = Visualisation () vis.plot_imfs (imfs=imfs, residue=res, tweet, include_residue=True) # draw and display the instantaneous frequency vis.plot_instant_freq (t, imfs=imfs) vis.show () of all provided IMF

At this point, I believe you have a deeper understanding of "Python how to achieve EMD algorithm", might as well come to the actual operation of it! Here is the website, more related content can enter the relevant channels to inquire, follow us, continue to learn!

Welcome to subscribe "Shulou Technology Information " to get latest news, interesting things and hot topics in the IT industry, and controls the hottest and latest Internet news, technology news and IT industry trends.