Electrooculogram processing method based on improved MIBSS-CEEMDAN

Abstract

The application relates to an electrooculogram processing method based on improved MIBSS-CEEMDAN, and belongs to the technical field of biomedical signal processing. The method proposes a fuzzy adaptive boundary information fitting method, solves the boundary defect phenomenon existing in EMD decomposition, designs a denoising method fusing double coefficients and double categories, solves the mode aliasing problem, proposes an ICA intelligent optimization criterion based on mutual information, further decomposes the IMF component, separates different types of electrooculogram signal components and noise components, and improves the accuracy and stability of source signal separation. The method can effectively realize decomposition of the electrooculogram signal and removal of the noise, and improves the accuracy and signal integrity of electrooculogram signal processing.

Description

Technical Field

[0001] This invention relates to an electrooculogram (EOG) signal processing method, specifically an improved fully adaptive noise-assisted ensemble empirical mode decomposition (CEEMDAN) combined with mutual information-based independent component analysis (MIBSS) for EOG signal processing, belonging to the field of biomedical signal processing technology. Background Technology

[0002] Electrooculography (EOG) is a bioelectrical signal that reflects the state of human eye movement by measuring the potential difference caused by eye movements. It has wide applications in fields such as human-computer interaction, virtual reality, and psychological research.

[0003] During the acquisition of electrooculogram (EOG) signals, various noises (such as electromyographic noise and power frequency interference) often interfere. These noises severely affect the accuracy and reliability of EOG signals. Therefore, how to effectively remove noise and extract pure EOG signals is a research hotspot and challenge in the field of EOG signal processing technology.

[0004] Traditional methods for denoising electrooculogram (EOG) signals (such as Fourier transform and wavelet transform) can effectively remove noise to some extent, but they still have some significant problems. For example, Fourier transform cannot reflect the frequency change of the signal over time when processing EOG signals, the choice of wavelet basis in wavelet transform has a great influence on the processing results of EOG signals, and it may produce large oscillations when processing EOG signals with abrupt changes, affecting the accuracy of denoising.

[0005] Fully Adaptive Noise-Assisted Ensemble Empirical Mode Decomposition (CEEMDAN) addresses mode aliasing by integrating the decomposition results of signals with added white noise. In Independent Component Analysis (ICA), source signals are assumed to be independent, and optimization algorithms separate the mixed signal into independent source signals. However, when processing electrooculogram (EOG) signals, CEEMDAN is still affected by noise, resulting in a significant amount of noise in the decomposed modal components. Furthermore, ICA faces challenges such as an unknown number of source signals and unstable separation performance. The Mutual Information-Based ICA Method (MIBSS) utilizes mutual information as an optimization criterion, improving the accuracy and stability of source signal separation, particularly demonstrating superior performance when processing EOG signals containing complex noise.

[0006] However, when using CEEMDAN or MIBSS alone to process electrooculogram signals, there are still problems such as noise interference and low system performance and quality. Summary of the Invention

[0007] The purpose of this invention is to address the technical problems of distortion and low accuracy in existing electrooculogram (EOG) signal processing technologies by creatively proposing an EOG signal processing method based on an improved MIBSS-CEEMDAN, which can effectively decompose EOG signals and remove noise, thereby improving the accuracy and signal integrity of EOG signal processing.

[0008] The innovative aspects of this invention include:

[0009] First, a fuzzy adaptive boundary information fitting method is proposed to solve the boundary defect phenomenon in EMD decomposition (Empirical Mode Decomposition).

[0010] Second, a denoising method that integrates two coefficients and two types was designed to solve the modal aliasing problem.

[0011] Third, an ICA intelligent optimization criterion based on mutual information is proposed to further decompose the IMF (Intrinsic Mode Function) components, separate different types of electrooculography signal components and noise components, and improve the accuracy and stability of source signal separation.

[0012] Beneficial effects

[0013] This method has the following advantages compared to existing technologies:

[0014] 1. This method proposes a fuzzy adaptive fitting boundary information method to generate flexible boundary events, which effectively solves the boundary loss problem in EMD decomposition, and more accurately fits and reconstructs the boundary information of the signal, thereby improving the accuracy and stability of the decomposition.

[0015] 2. This method designs a dual-coefficient, dual-type denoising method to solve the modal aliasing problem. By introducing dual-coefficient, dual-type fused noise, different modes can be effectively separated and distinguished, improving the accuracy and clarity of signal decomposition.

[0016] 3. This method proposes an ICA intelligent optimization criterion based on mutual information, which decomposes the signal matrix composed of IMF components and identifies and reduces artifact components unrelated to electrooculography (EOG) activity through mutual information analysis. This more accurately separates different types of EOG signal components and noise components, improving the accuracy and stability of source signal separation.

[0017] 4. This method has greater adaptability, can handle a wider range of electrooculogram (EOG) signal changes, adapts to the processing needs of EOG signals with different signal characteristics, solves the problem of large differences in EOG signals among different individuals and under different physiological states, and improves the accuracy and stability of source signal separation. Attached Figure Description

[0018] Figure 1 This is a schematic diagram of the process of the present invention;

[0019] Figure 2 This is a diagram illustrating the specific algorithm steps of the present invention;

[0020] Figure 3 This is a step diagram of the improved CEEMDAN decomposition method in this invention;

[0021] Figure 4 This is a detailed schematic diagram of the improved CEEMDAN decomposition in this invention;

[0022] Figure 5 The results of vertical / horizontal electrooculography (EOG) signal instance events detection;

[0023] Figure 6 This is a high / low mixed sampling result for a vertical electrooculogram signal example;

[0024] Figure 7 This is a high / low mixed sampling result for a horizontal electrooculography signal example;

[0025] Figure 8 The IMF component results obtained from the CEEMDAN decomposition of a vertical electrooculogram signal example are shown in the figure.

[0026] Figure 9 The image shows the IMF component results obtained from the CEEMDAN decomposition of a horizontal electrooculogram (EOG) signal. Detailed Implementation

[0027] The method of the present invention will be further described in detail below with reference to the accompanying drawings and embodiments.

[0028] The present invention is achieved using the following technical solution.

[0029] like Figure 1 , Figure 2 As shown, a method for processing electrooculograms based on an improved MIBSS-CEEMDAN includes the following steps:

[0030] Step 1: Fuzzy adaptive boundary information fitting.

[0031] The upper and lower boundaries of the preprocessed electrooculogram (EOG) signal are evaluated using a fuzzy decision tree algorithm to assess possible EOG events, thereby decomposing and generating reasonable signal boundaries and reducing boundary defects during the decomposition process.

[0032] Includes the following steps:

[0033] Step 1.1: Perform high / low mixed sampling on the original electrooculogram signal, where t is the continuous time, t is the low sampling rate, and t is the high sampling rate.

[0034] Step 1.1.1: Detect the start and end points of the electrooculogram (EOG) signal using the short-time energy change rate method.

[0035] First, frame and weight the EOG signal and calculate the energy of each frame. where \(x[i]\) is the sampled value in each frame of the signal, \(i = 1,2,\cdots,M\), \(M\) is the frame size, and \(w[i]\) is the Hamming window function.

[0036] Then, calculate the change rate between adjacent frames. where \(\epsilon\) is a very small positive number to avoid division by zero.

[0037] Define the short-time energy threshold \(STE=\mu\) E +\(\psi\sigma\) E , where \(\mu\) E and \(\sigma\) E are the mean and standard deviation of the frame energy \(E\) respectively, and \(\psi\) is used to adjust the increase degree of the threshold relative to the standard deviation.

[0038] For each frame index \(i\), if \(E\) i >\(STE\) and \(R\) i >\(\tau\), then mark it as the start of the event \(t\) s，i , until \(E\) i <\(STE\) or \(R\) i <1 / \(\tau\), mark it as the end of the event \(t\) e，i , where \(STE\) is the energy threshold and \(\tau\) is the change rate threshold.

[0039] Step 1.1.2: Perform position indexing according to the start and end points \(t\) s，i , \(t\) e，i of the event occurrence, and perform high / low mixed sampling on the EOG signal.

[0040] The sampled signal \(y\) f [t] is expressed as:

[0041]

[0042] where \(y\) L [n] is the low-sampling-rate \(f\) L sampled signal segment; \(y\) H [m] is the high-sampling-rate \(f\) H sampled signal segment;

[0043] Within the 2\(\Delta t\) range of the edge region between \(t\) s，i and \(t\) e，i of the event occurrence, linearly interpolate to obtain the sampled signal \(y\) l [t]:

[0044]

[0045] where \(t\) i is the start and end points \(t\) s，ior t e，i Δt is a very small sampling interval.

[0046] Step 1.2: Sample the signal y f [t], Design a fuzzy adaptive boundary information fitting method.

[0047] Step 1.2.1: Based on the event start and end points t obtained in Step 1.1.1 i Using the start and end positions of the signal as boundary points, calculate the boundary point t. b With the start and end points t of the event i Distance d, d = t b -t i Boundary point rate of change Where Δt is a very small sampling interval.

[0048] Step 1.2.2: Design fuzzy rules.

[0049] Define a distance fuzzy set A based on the magnitude of the distance d. far and A near Two fuzzy sets B of rate of change are defined based on the magnitude of the rate of change r at the boundary points. High and B Low Rules are formulated based on four possible scenarios:

[0050] If d∈A far r∈B Low If the signal changes smoothly at the boundary points, then the linear regression (LR) method is used.

[0051] If d∈A far r∈B High If the boundary point is in a situation where the signal changes drastically but the change is not directly related to the event itself, then the polynomial curve fitting (PCF) method is used.

[0052] If d∈A near r∈B Low If the boundary points are near the start and end points of the event but the signal changes relatively smoothly, then the support vector regression (SVR) method is used.

[0053] If d∈A near And r∈B High If the boundary point is located within the event occurrence area, then the non-linear least squares (NLS) method is used.

[0054] Right now:

[0055]

[0056] In this case, the distance membership function is set to a Gaussian function. e represents an exponential function; the membership function for the rate of change is the sigmoid function. r0 is the rate of change threshold; λ1 and λ2 are the thresholds for the two membership functions, respectively; D is the output variable decision.

[0057] Step 1.2.3: Design a fuzzy adaptive boundary fitting information method.

[0058] Specifically, let the training sample set be A = {(h1, g1), (h2, g2), ..., (h...}. l g l )}, where h i g i h is the i-th input / output vector. i =[y f (t i ), y f (t i+1 ), ..., y f (t N-l+i-1 )] T g i =[y f (t i ), y f (t i+1 ), ..., y f (t N-l+i )] T , 1≤i≤l, where l is the number of training samples, and T denotes matrix transpose. Test sample set B={(h N-s+1 g N-s+1 ), ..., (h N g N )}, 1≤s≤Nl is the number of test samples, where h N g N This is the Nth input / output vector.

[0059] The decision D is trained by selecting a fuzzy set that judges the distance d and the rate of change r. Based on Equation 3, the signal y f [t] performs right extension, outputting the right extension F = {y} f [t1],y f [t2],…,y f [t N ], y f [t N+1 ]}, where y f [t N+1 [ ] represents the first boundary signal of the rightward extension of the signal. Repeat the above steps to obtain the complete right-extended sequence y according to the required number of data points. f [tN+1 ]、y f [t N+2 ]、...、y f [t N+M Similarly, according to Equation 3, the signal is extended to the left to obtain the complete left-extended sequence y. f [t -M+1 ]、y f [t -M+2 ]、...、y f [t0], where M is the number of signal points extended to the left / right.

[0060] Finally, by sorting the left-extended sequence, the original signal sequence, and the right-extended sequence, a new signal sequence y′(t) is formed, where y′(t) = {y f [t -M+1 ],…,y f [t0],y f [t1],…,y f [t N ],…,y f [t N+M ]}.

[0061] Step 2: For Gaussian white noise with positive and negative pairs and 1 / f noise, a dual-coefficient dual-class denoising method is designed to improve the CEEMDAN decomposition of the electrooculogram signal.

[0062] like Figure 3 As shown, step 2 includes the following steps:

[0063] Step 2.1: Add paired Gaussian white noise and 1 / f noise N times to the extended signal y′(t) to be decomposed, and perform decomposition to ensure that the number of local extrema and zero-crossings is equal or differs by at most one throughout the entire time range, and the average of the upper and lower envelopes of the local extrema at any time point is zero.

[0064]

[0065] Where E(·) represents EMD decomposition, δ, For two noise amplitude coefficients, g i (t) represents a Gaussian white noise signal that follows a standard normal distribution, f i (t) represents 1 / f noise, i = 1, 2, ..., N, where i is the number of times noise is added, γ i，1 This is the remaining amount after signal decomposition.

[0066] The N modal components generated are averaged as a whole to obtain

[0067] Step 2.2: Through singular value decomposition, the residual signal γ is transformed.i，1 Remove the linear correlation of the signal, i.e., {γ i，1}=U∑V * , where U and V represent matrices spanned by eigenvectors, and ∑ represents the singular value matrix.

[0068] Select V * The columns corresponding to the non-zero singular values ​​are represented as a set of linearly independent vectors {v1, v2, ..., v} in the Hilbert space. m},pass The standard orthogonal system {e1, e2, ..., e} is calculated. m}, where m <N。

[0069] Step 2.3: Decomposition of the k-th time, for the remainder r k-1 (t) Add paired Gaussian white noise and 1 / f noise and decompose

[0070] The IMF components of the k-th decomposition are obtained by calculating the population mean.

[0071] Step 2.4: Calculate the residual margin r k (t)=r k-1 (t)-IMF k (t), using the method in step 2.2, the standard orthogonal system {e} is obtained. i，k Constructing orthogonal approximations Ensure that the L2 norm of the error between the residual signal and the fitted signal is less than ε; otherwise, reselect the linearly independent vector group, i.e.:

[0072]

[0073] Among them, E k This represents the error value between the residual signal and the fitted signal. <r k e i，k > represents the inner product in Hilbert space; ε represents the error threshold.

[0074] Step 2.5: Repeat the above steps. If the residual signal r after the Kth decomposition... K (t) Orthogonal approximation function If the algorithm is monotonic and has no local extrema, and cannot be further decomposed, then the iteration stops. The specific algorithm implementation steps are as follows: Figure 4 As shown.

[0075] Step 2.6: The decomposed signal is represented as follows:

[0076]

[0077] Step 3: Based on the ICA intelligent optimization criterion under mutual information, further decompose the IMF modal components to obtain different types of electrooculography signal components and noise components.

[0078] Includes the following steps:

[0079] Step 3.1: For each IMF modal component, the IMF... i (t), calculate the correlation coefficient with the original signal. The expression for calculating the correlation coefficient is:

[0080]

[0081] Where Ttp represents the total number of time points in the signal sequence y(t).

[0082] Based on the absolute value of the correlation coefficient |R i |<ρ2、ρ2≤|R i |<ρ1、|R i |≥ρ1, divided into three different demand levels: fine selection, secondary fine selection and coarse selection, where ρ1 and ρ2 are two different correlation coefficient thresholds.

[0083] Select the IMF components, combine them, and sort them by column to form new signal matrices X1, X2, and X3. Each column in the matrix represents a component, and the rows represent the time series. The signal matrix is ​​represented as follows:

[0084]

[0085] Step 3.2: Analyze the independent source components of the signal matrix.

[0086] First, the signal is processed by removing the mean and whitening it so that its mean is 0 and its variance is 1.

[0087] The correlation between observed signals is removed by linear transformation. After processing, a function is selected for approximation, and the comparison function is expressed as follows:

[0088]

[0089] Where a1 represents a constant variable that affects the steepness or proportion of the logarithm of the hyperbolic cosine function; u represents the independent variable.

[0090] Step 3.3: Use the direction in which the negative entropy increases the fastest, i.e. the gradient of the negative entropy, as the initial value of W, W←E{xG′(Wx)}, where G′ represents the derivative of G(Wx).

[0091] The output signal S is calculated using the current separation matrix W. i =WX i For each output signal component S i ={s1, s2, ..., s mEstimate the probability density function p i (s i ), m is each signal matrix X i Number of IMF components in China.

[0092] Step 3.4: Iterative optimization.

[0093] For signal matrices X1, X2, and X3, select the most suitable optimization criterion to improve the effectiveness of ICA decomposition and ensure that the decomposed components are as independent as possible at a specific level.

[0094] Step 3.4.1: For the selected signal matrix X1, select the log-likelihood maximization method to extract independent components with specific physiological significance from the complex electrooculogram data. Set the expression for the log-likelihood function:

[0095]

[0096] Where, p i (s i ) is the probability density function, det(W) represents the determinant of the separation matrix W, and λ is the regularization strength parameter.

[0097] Update W1 using an adaptive learning rate method, i.e.:

[0098]

[0099] Where α is the learning rate; ∈ is a minimal constant to ensure the denominator is not zero; and t is the number of iterations. and The first and second moment estimates after bias correction are respectively, and the update expressions are:

[0100]

[0101]

[0102]

[0103] in, This indicates a correction to the first-order moment estimate; β1 is the gradient of the log-likelihood function with respect to W1; β2 and β1 are the iteration parameters.

[0104] Step 3.4.2: For the second-selected signal matrix X2, select the method of minimizing mutual information to ensure that the extracted signal components have high statistical independence.

[0105] Based on the estimated probability density function, the mutual information between the components of the output signal is calculated. Considering the scale uncertainty problem, the mutual information I(S) is approximately calculated as follows:

[0106]

[0107] Where J(S) represents negative entropy, J(S) = [E{G(S)} - E{G(v)}] 2 , where S is the output signal, v represents a Gaussian random vector with zero mean and unit variance, and m represents the number of components in the output signal matrix.

[0108] Gradient of mutual information I(S) with respect to the separation matrix W Update the separation matrix W2, and obtain the iterative formula for W2 as follows:

[0109]

[0110] Where η is the learning rate, and h is a function that adjusts the learning rate according to the proportion of performance improvement between two consecutive iterations.

[0111] Step 3.4.3: For the coarse-selected signal matrix X3, update the separation matrix W3 using Newton's method to obtain the iterative formula:

[0112] W3←E{XG′(W3X)}-E{G″(W3X)}W3 (17)

[0113] Where G′ represents the first derivative of the comparison function G(·), G″ represents the second derivative of the comparison function G(·), and X represents the signal matrix.

[0114] Step 3.4.4: Normalize W:

[0115] Iterative calculation, if ΔW is less than the set convergence threshold If the maximum number of iterations is reached, the iteration will stop.

[0116] The final independent signal component is S i =W i X i W i For the unmixing matrix, X i It is a signal matrix composed of IMF components.

[0117] Step 4: Perform noise analysis on the obtained independent signal components to achieve signal denoising.

[0118] Observe the time series, spectral characteristics, and spatial distribution of each independent component. Set the components with obvious energy concentration in a specific frequency band, as well as those exhibiting spatial distribution characteristics consistent with Gaussian white noise and 1 / f noise, to zero, to obtain S1′, S2′, and S3′, respectively.

[0119] The independent signal component matrices are concatenated in sequence as S′ = [S1′ S2′ S3′]. The reconstructed new signal matrix is X′ = W -1 S′, where X′ = [IMF1′, IMF2′, …, IMF K ′].

[0120] The reconstruction combination of IMF components has the following calculation formula:

[0121]

[0122] where y′(t) is the denoised signal.

[0123] Example

[0124] This example describes the implementation process of applying the method of the present invention to the scenario of electrooculogram signal processing.

[0125] In this example, the data is from the publicly available dataset DEAP database. This database is based on physiological signals generated under the induction of music video materials, recording the physiological signals of 32 subjects watching 40 minutes of music videos (each music video is 1 minute) and the psychological scales of the subjects. The signal length of the collected electrooculogram signals is 63 s, the signals are downsampled to 128 Hz, and are preliminarily band-pass filtered to 4 - 45 Hz.

[0126] Figure 1 is the schematic diagram of the process of this method, Figure 2 is the specific algorithm step diagram of this method. It can be seen from the figure that the preprocessed electrooculogram signals are first obtained, and then the following steps are carried out:

[0127] Step 1: Extend the boundary of the collected electrooculogram signals by the method of fuzzy adaptive boundary fitting.

[0128] Step 1.1: Perform hybrid sampling on the electrooculogram signals.

[0129] Detect the start and end points of events of the electrooculogram signals according to the short-time energy change rate. The frame size M = 50, the step size of each frame is 10, the increase degree ψ of the threshold relative to the standard deviation is 0.8, and the energy threshold T = μ E + ψσ E . For each frame index i, if E i > T and R i > τ, it is marked as the start of the event t s，i , until E i < T or R i < 1 / τ, it is marked as the end of the event t e，i , where the change rate threshold τ = 1.2. The event detection results are as Figure 5 shown.

[0130] Use a low sampling rate f L =64Hz coarse sampling of the electrooculogram signal to obtain a preliminary downsampled signal y L [n]. Based on the start and end points of the electrooculogram (EOG) signal events, a higher sampling rate f is used. H =512Hz fine sampling to obtain y H [m]. Linear interpolation using Equation 2 is used to transition between high and low sampling rates. The final mixed-sample signal y is formed according to Equation 1. f [t].

[0131] In this specific example, the original number of sampling points was 8064, used to calculate short-time energy detection events. After downsampling, the number of sampling points was n = 4032. There were 9 events occurring in the horizontal electrooculogram (EOG) signal within 63 seconds, resulting in a final mixed sampling point count of k = 4386; there were 11 events occurring in the vertical EOG signal within 63 seconds, resulting in a final mixed sampling point count of k = 4252. The mixed sampling signal diagram is shown below. Figure 6 , 7 As shown.

[0132] Step 1.2: Perform fuzzy adaptive fitting boundary extension on the sampled signal.

[0133] Calculate the distance d between the boundary point and the start and end points, as well as the rate of change r of the boundary point, and use them as inputs to the fuzzy model. Then, use Equation 3 to obtain the output fuzzy decision D. Based on decision D, select the fitting and extension method, and extend left and right respectively to obtain the complete boundary extension sequence y′(t)={y f [t -M+1 ],…,y f [t0],y f [t1],…,y f [t N ],…,y f [t N+M ]}, where M is the number of signal points extended to the left / right.

[0134] In this specific example, the fuzzy set is set with λ1 = 0.683 and λ2 = 0.5; the membership function μ rat In this example, r0 is set to 0.6. 128 sampling points are extended outwards from both ends of the signal, i.e., M = 128, for a total of 256 sampling points.

[0135] Step 2: Improved CEEMDAN decomposition of electrooculogram signals is achieved using a dual-coefficient, dual-type fusion noise method.

[0136] Step 2.1: Using Equation 4, add the Gaussian white noise signal with paired positive and negative values ​​that satisfy the standard normal distribution and 1 / f noise to the extended signal to be decomposed, and decompose the signal after the i-th addition of noise, and calculate the overall mean IMF1. Integrate the linearly independent vector group to obtain the orthogonal system {e1, e2, ..., e... m}

[0137] In this specific example, the noise amplitude coefficient ε = 0.01. First, add it 100 times.

[0138] Step 2.2: Decomposition of the k-th time, for the remainder r k-1 (t) Adding paired Gaussian white noise and 1 / f noise yields a new signal. EMD decomposition is then performed on the new signal to obtain the IMF component. k，i (t), the IMF components are calculated from the overall mean, and the residuals are calculated using the standard orthogonal system {e i，k Constructing orthogonal approximations of signals Equation 5 is used to calculate the error between the residual signal and the fitted signal, ensuring that the error ε is less than 0.5.

[0139] Step 2.3: Perform the above processing and iteration on the residual signal to determine the orthogonal approximation signal function r. K If ′(t) is a monotonic function, then the iteration stops and the decomposition of the improved CEEMDAN algorithm ends.

[0140] In this specific example, CEEMDAN decomposition was performed on both the vertical and horizontal electrooculogram (EOG) signals, yielding 11 IMF modes respectively. The reconstructed signal was then calculated using Equation 6. The results of the signal decomposition are as follows: Figure 8 , 9 As shown.

[0141] Step 3: Based on the ICA intelligent optimization criterion under mutual information, the IMF component combination is further decomposed, and noise components are screened and removed.

[0142] Step 3.1: Calculate the correlation coefficient between the IMF modal components and the original signal, based on the magnitude of the absolute value of the correlation coefficient |R i |>0.8、0.5<|R i |<0.8、|R i |<0.5, select the modal component IMF respectively i (t), sorted by column to form new signal matrices X1, X2, X3, matrix representation 8, which are the selected matrix, the secondary selected matrix, and the coarse selected matrix, respectively.

[0143] Step 3.2: Perform ICA decomposition preprocessing on the signal, approximate the function using Equation 9, and use the gradient of the negative entropy as the initial value of W, i.e., W←E{xG′(Wx)}.

[0144] Step 3.3: Iterative Optimization. Use the current separation matrix W to calculate the output signal S = WX, where each output signal component S = {s1, s2, ..., s...} M Estimate the probability density function p i (s i ).

[0145] For the selected signal matrix X1, based on the estimated probability density function, the log-likelihood function is calculated using the method in Equation 10, W1 is updated using the adaptive learning rate method in Equation 11, and the parameters in Equation 11 are updated using Equations 12-14. Initialization: α = 0.001, m0 = 0, v0 = 0, β1 = 0.9, β2 = 0.999, ∈ = 10 -8 .

[0146] For the subselected signal matrix X2, the mutual information between the output signal components is calculated using Equation 15. The separation matrix W2 is updated using Equation 16. Initialize h = 1, η = 0.01.

[0147] For the coarsely selected signal matrix X3, the separation matrix is ​​updated using Newton's method, i.e., iterative formula 17.

[0148] Step 3.4: For W i Normalization Iterate through the calculations until the stopping condition is met.

[0149] Step 4: Perform noise analysis on the obtained independent signal components. Components with significant energy concentration within a specific frequency band, as well as components exhibiting spatial distribution characteristics consistent with Gaussian white noise and 1 / f noise, are set to zero. Signal reconstruction is then performed, and the reconstructed signal is obtained using Equation 18.

Claims

1. A method for processing electrooculograms based on the improved MIBSS-CEEMDAN, characterized in that, Includes the following steps: Step 1: Fuzzy adaptive boundary information fitting; The upper and lower boundaries of the preprocessed electrooculogram (EOG) signal are evaluated using a fuzzy decision tree algorithm to assess possible EOG events, thereby decomposing and generating reasonable signal boundaries and reducing boundary defects during the decomposition process. Step 2: For Gaussian white noise with positive and negative pairs and 1 / f noise, a denoising method with dual coefficients and dual categories is designed to improve the CEEMDAN decomposition of the electrooculogram signal. Step 2.1: Add paired Gaussian white noise and 1 / f noise N times to the extended signal y′(t) to be decomposed, and perform decomposition to ensure that the number of local extrema and zero-crossings is equal or differs by at most one throughout the entire time range, and the average of the upper and lower envelopes of the local extrema at any time point is zero. Where E(·) represents EMD decomposition, δ, For two noise amplitude coefficients, g i (t) represents a Gaussian white noise signal that follows a standard normal distribution, f i (t) represents 1 / f noise, i = 1, 2, ..., N, where i is the number of times noise is added, γ i，1 The remaining amount after signal decomposition; The N modal components generated are averaged as a whole to obtain Step 2.2: Through singular value decomposition, transform the residual signal γ i，1 Remove the linear correlation of the signal, i.e., {γ i，1 }=U∑V * , where U and V represent matrices spanned by eigenvectors, and Σ represents the singular value matrix; Select V * The columns corresponding to the non-zero singular values ​​are represented as a set of linearly independent vectors {v1, v2, ..., v} in the Hilbert space. m },pass The standard orthogonal system {e1, e2, ..., e} is calculated. m }, where m < N; Step 2.3: Decomposition of the k-th time, for the remainder r k-1 (t) Add paired Gaussian white noise and 1 / f noise and decompose The IMF components of the k-th decomposition are obtained by calculating the population mean. Step 2.4: Calculate the residual margin r k (t)=r k-1 (t)-IMF k (t), using the method in step 2.2, the standard orthogonal system {e} is obtained. i，k Constructing orthogonal approximations Ensure that the L2 norm of the error between the residual signal and the fitted signal is less than ε; otherwise, reselect the linearly independent vector group, i.e.: Among them, E k This represents the error value between the residual signal and the fitted signal. <r k e i，k > represents the inner product in Hilbert space; ε represents the error threshold. Step 2.5: Repeat the above steps; if the residual signal r after the Kth decomposition... K (t) Orthogonal approximation function If the value is monotonic and has no local extrema, and cannot be further decomposed, then the iteration stops. Step 2.6: The decomposed signal is represented as follows: Step 3: Based on the ICA intelligent optimization criterion under mutual information, further decompose the IMF modal components to obtain different types of electrooculography signal components and noise components; Step 3.1: For each IMF modal component, the IMF... i (t), calculate the correlation coefficient with the original signal. The expression for calculating the correlation coefficient is: Where Ttp represents the total number of time points in the signal sequence y(t); Based on the absolute value of the correlation coefficient |R i |<ρ2、ρ2≤|R i |<ρ1、|R i |≥ρ1, which is divided into three different demand levels: fine selection, secondary fine selection and coarse selection, where ρ1 and ρ2 are two different correlation coefficient thresholds; Select the IMF components, combine them, and sort them by column to form new signal matrices X1, X2, X3; each column in the matrix represents a component, and the rows represent the time series. The signal matrix is ​​represented as follows: Step 3.2: Analyze the independent source components of the signal matrix; First, the signal is processed by removing the mean and whitening it so that its mean is 0 and its variance is 1. The correlation between observed signals is removed by linear transformation. After processing, a function is selected for approximation, and the comparison function is expressed as follows: Where a1 represents a constant variable that affects the steepness or proportion of the logarithm of the hyperbolic cosine function; u represents the independent variable; Step 3.3: Take the direction of the fastest increase in negative entropy, i.e. the gradient of negative entropy, as the initial value of W, W←E{xG′(Wx)}, where G′ represents the derivative of G(Wx); The output signal S is calculated using the current separation matrix W. i =WX i For each output signal component S i ={s1, s2, ..., s m Estimate the probability density function p i (s i ), m is each signal matrix X i Number of IMF components in China; Step 3.4: Iterative optimization; For signal matrices X1, X2, and X3, select the most suitable optimization criterion to improve the effectiveness of ICA decomposition and ensure that the decomposed components are as independent as possible at a specific level. Step 3.4.1: For the selected signal matrix X1, select the log-likelihood maximization method to extract independent components with specific physiological significance from the complex electrooculogram data; set the expression for the log-likelihood function: Where, p i (s i ) is the probability density function, det(W) represents the determinant of the separation matrix W, and λ is the regularization strength parameter; Update W1 using an adaptive learning rate method, i.e.: Where α is the learning rate; ∈ is a minimal constant to ensure the denominator is not zero; and t is the number of iterations. and The first and second moment estimates after bias correction are respectively, and the update expressions are: in, This indicates a correction to the first-order moment estimate; β1 is the gradient of the log-likelihood function with respect to W1; β2 and β1 are the iteration parameters; Step 3.4.2: For the suboptimal signal matrix X2, select the method that minimizes mutual information; Based on the estimated probability density function, the mutual information between the components of the output signal is calculated; the mutual information I(S) is approximately calculated as follows: Where J(S) represents negative entropy, J(S) = [E{G(S)} - E{G(v)}] 2 Where S is the output signal, v represents a Gaussian random vector with zero mean and unit variance, and m represents the number of components in the output signal matrix; Gradient of mutual information I(S) with respect to the separation matrix W Update the separation matrix W2, and obtain the iterative formula for W2 as follows: Where η is the learning rate, and h is a function that adjusts the learning rate according to the performance improvement ratio between two consecutive iterations; Step 3.4.3: For the coarse-selected signal matrix X3, update the separation matrix W3 using Newton's method to obtain the iterative formula: W3←E{XG′(W3X)}-E{G″(W3X)}W3 Where G′ represents the first derivative of the comparison function G(·), G″ represents the second derivative of the comparison function G(·), and X represents the signal matrix; Step 3.4.4: Normalize W: Iterative calculation, if ΔW is less than the set convergence threshold If the maximum number of iterations is reached, then stop iterating; The final independent signal component is S i =W i X i W i For the unmixing matrix, X i The signal matrix is ​​composed of IMF components; Step 4: Perform noise analysis on the obtained independent signal components to achieve signal denoising.

2. The electrooculography processing method based on the improved MIBSS-CEEMDAN as described in claim 1, characterized in that, Step 1 includes the following steps: Step 1.1: Perform high / low mixed sampling on the original electrooculography signal y(t), where t is continuous time and f is the frequency of the signal. L For low sampling rate, f H For high sampling rate; Step 1.1.1: Use the short-time energy change rate method to detect the start and end points of events in the electrooculogram (EOG) signal; First, the electrooculogram (EOG) signal is segmented and weighted, and the energy of each frame is calculated. Where y i [n] represents the sampled value in each frame of the signal, n = 0, 1, ..., M-1, M is the frame size, and ω[n] is the Hamming window function; Then, calculate the rate of change for each adjacent frame. Where ζ is a very small positive number, avoid dividing by zero; Define the short-time energy threshold STE = μ E +ψσ E , where μ E σ E , where are the mean and standard deviation of the frame energy E, respectively, and ψ is used to adjust the degree of increase of the threshold relative to the standard deviation; For each frame index i, if E i >STE and R i If the value is greater than τ, then it is marked as the start of event t. s，i until E i <STE or R i <1 / τ, marked as the end of event t e，i , where STE is the energy threshold and τ is the rate of change threshold; Step 1.1.2: Based on the start and end points t of the event s，i t e，i Position indexing is performed, and high / low mixed sampling of the electrooculogram signal is performed; Sampled signal y f [t] is represented as: Among them, y L [n] represents the low sampling rate f L Sampling signal segment; y H [m] represents the high sampling rate f H Sampling signal segment; In t s，i and t e，i Within the 2Δt range of the edge region, the sampled signal y is obtained by linear interpolation. l [t]: Among them, t i Let t be the start and end points of the event. s，i or t e，i Δt is a very small sampling interval; Step 1.2: Sample the signal y f [t], Design a fuzzy adaptive boundary information fitting method; Step 1.2.1: Based on the event start and end points t obtained in Step 1.1.1 i Using the start and end positions of the signal as boundary points, calculate the boundary point t. b With the start and end points t of the event i Distance d, d = t b -t i Rate of change at boundary points Where Δt is a very small sampling interval; Step 1.2.2: Design fuzzy rules; Define a distance fuzzy set A based on the magnitude of the distance d. far and A near Two fuzzy sets B of rate of change are defined based on the magnitude of the rate of change r at the boundary points. High and B Low Rules are formulated based on four possible scenarios: If d∈A far r∈B Low If the signal changes gradually at the boundary points, then the linear regression method LR is used. If d∈A far r∈B High If the boundary point is in a situation where the signal changes drastically but the change is not directly related to the event itself, then the polynomial fitting method PCF is used. If d∈A near r∈B Low If the boundary point is near the start and end point of the event but the signal change is relatively smooth, then the support vector regression (SVR) method is used. If d∈A near And r∈B High If the boundary point is located within the event occurrence area, then nonlinear least squares (NLS) is used. In this case, the distance membership function is set to a Gaussian function. e represents an exponential function; the membership function for the rate of change is the sigmoid function. r0 is the rate of change threshold; λ1 and λ2 are the thresholds for the two membership functions, respectively; D is the output variable decision. Step 1.2.3: Design a fuzzy adaptive boundary fitting information method; Let the training sample set be A = {(h1, g1), (h2, g2), ..., (h... l g l )}, where h i g i For the i-th input / output vector, h i =[y f (t i ), y f (t i+1 ), ..., y f (t N-l+i-1 )] T g i =[y f (t i ), y f (t i+1 ), ..., y f (t N-l+i )] T , 1≤i≤l, where l is the number of training samples, and T represents the matrix transpose; test sample set B={(h N-s+1 g N-s+1 ), ..., (h N g N )}, 1≤s≤Nl is the number of test samples, where h N g N This is the Nth input / output vector; The decision D is trained by selecting a fuzzy set based on the distance d and the rate of change r; based on Equation 1, the signal y is... f [t] performs right extension, outputting the right extension F = {y} f [t1],y f [t2],…,y f [t N ], y f [t N+1 ]}, where y f [t N+1 [ ] represents the first boundary signal of the rightward extension of the signal; repeat the above steps to obtain the complete right-extended sequence y according to the required number of data points. f [t N+1 ]、y f [t N+2 ]、…、y f [t N+M Similarly, according to Equation 1, the signal is extended to the left to obtain the complete left-extended sequence y. f [t -M+1] y f [t -M+2 ]、…、y f [t0], where M is the number of signal points extended to the left / right; Finally, by sorting the left-extended sequence, the original signal sequence, and the right-extended sequence, a new signal sequence y′(t) is formed, where y′(t) = {y f [t -M+1 ],…,y f [t0],y f [t1],…,y f [t N ],…,y f [t N+M ]}.

3. The electrooculography processing method based on the improved MIBSS-CEEMDAN as described in claim 1, characterized in that, Step 4 includes the following steps: Observe the time series, spectral characteristics and spatial distribution of each independent component; set the components with obvious energy concentration in a specific frequency band and the components that exhibit characteristics consistent with Gaussian white noise and 1 / f noise in spatial distribution to zero, and obtain S1′, S2′ and S3′ respectively. The independent signal component matrices are concatenated sequentially as S′=[S1′ S2′ S3′]; the reconstructed new signal matrix is ​​X′=W⁻¹S′, where X′=[IMF1′, IMF2′, …, IMF K ′]; The IMF component reconstruction combination is given by the following formula: Where y′(t) is the denoised signal.

Images