Single-channel real-time artifact rejection method and system for electroencephalogram signals

Abstract

The application discloses a single-channel electroencephalogram signal real-time artifact suppression method and system, frames the collected signal first, extracts the power frequency stable component to establish a real-time calibration baseline, quantifies the baseline drift amplitude frame by frame according to the baseline, locates the drift steep change site, extracts the artifact segment with electromyographic characteristics in combination with the waveform and the shape consistency degree of the template, differentiates the trap wave depth according to the overlapping degree of the artifact and the electroencephalogram frequency band, and delimits the filter pollution boundary of the ringing tail area introduced by the trap wave itself, thereby defining the protection partition where the residual artifact and the brain wave coexist; interpolation repair is implemented on the damaged segment in the delimited safe segment in combination with the phase characteristics of the inherent rhythm of the electroencephalogram, the signal overall quality rating is given according to the length proportion of the suspected falling segment, local degradation marking is made on the overlapping paragraph which is difficult to separate in the frequency domain, and explicit gap annotation is made on the segment which cannot be reliably suppressed, and the accuracy of artifact identification and suppression is improved.

Description

Technical Field

[0001] This invention relates to the field of biomedical signal processing technology, and in particular to a method and system for real-time artifact suppression of single-channel electroencephalogram (EEG) signals. Background Technology

[0002] Electroencephalogram (EEG) signals reflect the electrical activity characteristics of the brain under different states and are widely used physiological signals in brain-computer interfaces, epilepsy monitoring, and anesthesia depth assessment. During acquisition and transmission, EEG signals are easily contaminated with artifacts from various sources, such as power supply interference, muscle activity (e.g., masticatory muscles), and unstable electrode contact. These artifacts often overlap with the true EEG rhythms in terms of amplitude and frequency. If they cannot be identified and suppressed in time, they will affect the reliability of subsequent signal interpretation and real-time decision-making. This artifact contamination problem is particularly prominent in applications requiring continuous, real-time output.

[0003] Existing artifact processing methods either rely on offline batch analysis of the entire recording, which is difficult to meet the processing latency requirements of real-time scenarios; or they use fixed-parameter filtering schemes, which may weaken the real EEG components while suppressing artifacts, making it difficult to balance the suppression effect and signal fidelity. In cases where artifacts cause signal loss, existing simple interpolation methods often ignore the inherent rhythmic characteristics of the EEG itself, resulting in significant differences between the repaired results and the real physiological waveforms. Furthermore, existing solutions generally lack quantitative evaluation methods for the overall reliability of the processed signal, making it difficult for downstream applications to determine whether the signal is reliable at a certain time. Summary of the Invention

[0004] This invention provides a method and system for real-time artifact suppression of single-channel EEG signals, aiming to solve the problems of existing solutions struggling to balance real-time performance, suppression accuracy, and rhythm fidelity, and lacking reliable result assessment. This method establishes a power frequency calibration baseline and identifies artifact characteristics, differentially suppresses interfering components in the signal, repairs damaged segments by combining EEG rhythm characteristics, and finally outputs an artifact-suppressed signal with quality assessment information. This improves the accuracy and usability of artifact suppression while maintaining real-time processing capabilities.

[0005] The first aspect of this invention proposes a method for real-time artifact suppression of single-channel electroencephalogram (EEG) signals, comprising: The single-channel real-time EEG acquisition sequence is processed into frames to obtain an EEG frame sequence. The power frequency stable segment is extracted from the EEG frame sequence to establish a real-time calibration baseline. Based on the real-time calibration baseline, frame-by-frame baseline drift amplitude quantization is performed to form a drift amplitude sequence. The drift abrupt change site is located according to the drift amplitude sequence to generate an EMG artifact candidate set. The template morphology highly consistent segments are labeled for the EMG artifact candidate set to obtain an artifact morphology feature set. Based on the artifact morphology feature set, the output notch filter parameters are reduced by depth differentiation in the high-overlap frequency band. The notch filter parameters are used to analyze the self-ringing tail span to obtain the filtered contamination section. Based on the filtered contamination section, the boundary between the coexistence of residual artifacts and brain waves is delineated to obtain the residual protection zone. The remaining protection partition is delineated according to the degree of separability to construct a repair candidate sequence by defining the boundary of the safe repair segment. Rhythmic phase anchoring completion is performed on the repair candidate sequence to form an interpolation repair sequence. The interpolation repair sequence is used to determine the suspicious silent segment based on the segment with the lowest repair rate. Based on the analysis of the suspected silent segments, the risk of detachment and coverage is analyzed to output signal quality rating. For the signal quality rating, the difficult-to-separate overlapping segments are downgraded and a graded EEG sequence is established. Based on the graded EEG sequence, the unsuppressed segments are explicitly marked with gaps to output a suppressed EEG sequence.

[0006] A second aspect of this invention provides a real-time artifact suppression system for single-channel electroencephalogram (EEG) signals, comprising: The calibration establishment unit is used to perform frame-segmentation processing on the single-channel real-time EEG acquisition sequence to obtain an EEG frame sequence, and to use the EEG frame sequence to extract the power frequency stable segment to establish a real-time calibration baseline. The feature collection unit is used to perform frame-by-frame baseline drift amplitude quantization based on the real-time calibration baseline to form a drift amplitude sequence, locate the drift abrupt change site based on the drift amplitude sequence to generate an EMG artifact candidate set, and mark the template morphology highly consistent segment for the EMG artifact candidate set to obtain an artifact morphology feature set. The filtering protection unit is used to implement differentiated reduction of output notch filter parameters in high-overlap frequency bands based on the artifact morphology feature set, analyze the self-ringing tail span through the notch filter parameters to obtain the filtered contamination section, and delineate the coexistence boundary between residual artifacts and brain waves based on the filtered contamination section to obtain the residual protection zone. The interpolation repair unit is used to delineate the boundaries of the safe repair segments according to the degree of separability of the residual protection partition to construct a repair candidate sequence, perform rhythmic phase anchoring completion on the repair candidate sequence to form an interpolation repair sequence, and use the interpolation repair sequence to determine the suspicious silent segment based on the segment with the lowest repair rate. The sequence output unit is used to analyze the risk of detachment and coverage based on the suspected silent segment, output signal quality rating, perform downgrading labeling of difficult-to-separate overlapping segments to establish a graded EEG sequence, and perform explicit gap labeling of unsuppressed segments to output a suppressed EEG sequence based on the graded EEG sequence.

[0007] The beneficial effects of this invention are reflected in the following points: 1. Traditional methods often rely on a single amplitude or slope threshold to determine whether artifacts have occurred, easily confusing gradual baseline drift with sudden muscle contraction; this invention first establishes a calibration baseline using a stable power frequency segment and then separates the power frequency, allowing the directional shift of baseline drift and the instantaneous jumps of electromyography to be separated on the same sequence, and then superimposed waveform morphology matching for secondary confirmation. Candidate segments must first trigger the jump criterion and then pass morphological verification before they can be included, no longer relying on single-dimensional numerical fluctuations, thus improving the accuracy of artifact recognition. 2. Existing notch filter designs suppress the target frequency band at a uniform depth, and consider the processing complete once the suppression intensity is sufficient, but they ignore the fact that the greater the depth, the more likely the filter's own impulse response tailing is to exceed expectations, thereby contaminating adjacent data; this invention, based on this, inversely calculates the upper limit of the notch filter depth according to the proportion of brain waves, and delineates the contaminated segment for each notch filter's measured impulse response, reducing the damage of the notch filter to the brainwave rhythm and preventing the filter's own distortion from being used as a valid signal. 3. Interpolation repair and defect labeling are not mutually exclusive processing methods. While completing the repair by planning the interpolation path according to the rhythm phase, this invention further judges which segments are unreliable based on the degree of the repair rate being lower than the median level and the silent characteristics of the original waveform. These segments are explicitly labeled as gaps instead of being quietly left in the processed data. Compared with simply using interpolation to cover up defects or discarding the entire segment, this improves the consistency between the repaired segments and the real EEG rhythm and enhances the traceability of output quality judgment. Attached Figure Description

[0008] Figure 1 This is a flowchart of a real-time artifact suppression method for single-channel EEG signals according to the present invention.

[0009] Figure 2 This is a schematic diagram of the measured waveforms of notch depth distribution and ringing trailing pollution of the present invention.

[0010] Figure 3 This is a structural block diagram of a single-channel EEG signal real-time artifact suppression system according to the present invention.

[0011] Among them: 1-scalp EEG signal waveform; 2-starting time of notch application; 3-ringing attenuation tail segment; 4-theoretical expected duration boundary; 5-measured trailing contamination segment; 6-frequency axis; 7-notch depth curve; 8-baseline depth line; 9-high-risk notch interval; 10-correspondence indicator line. Detailed Implementation

[0012] In the following description, specific details such as particular system architectures and techniques are set forth for illustrative purposes and not for limitation, in order to provide a thorough understanding of the embodiments of this application. However, those skilled in the art will understand that this application may also be implemented in other embodiments without these specific details. In other instances, detailed descriptions of well-known systems, apparatuses, circuits, and methods have been omitted so as not to obscure the description of this application with unnecessary detail.

[0013] References to "one embodiment" or "some embodiments" as described in this specification mean that one or more embodiments of this application include a specific feature, structure, or characteristic described in connection with that embodiment. Therefore, the phrases "in one embodiment," "in some embodiments," "in other embodiments," "in still other embodiments," etc., appearing in different parts of this specification do not necessarily refer to the same embodiment, but rather mean "one or more, but not all, embodiments," unless otherwise specifically emphasized. The terms "comprising," "including," "having," and variations thereof mean "including but not limited to," unless otherwise specifically emphasized.

[0014] The technical solutions of the embodiments of this application are described below.

[0015] like Figure 1 As shown, this embodiment of the invention provides a method for real-time artifact suppression of single-channel EEG signals, including the following steps S101-S105: Step S101: The single-channel real-time EEG acquisition sequence is processed into EEG frame sequence by framing, and the power frequency stable segment is extracted from the EEG frame sequence to establish a real-time calibration baseline.

[0016] Specifically, the single-channel real-time EEG acquisition sequence is framed to obtain the EEG frame sequence. The single-channel real-time EEG acquisition sequence is taken from the continuous output of a single lead electrode after pre-amplification and analog-to-digital conversion, with a sampling rate of 1000Hz and a quantization bit depth of 24 bits. No notch filtering or high-pass pre-processing is applied at the acquisition end, and the power frequency component and low-frequency drift caused by electrode polarization are preserved as is. The cutoff frequency of the pre-amplifier analog low-pass filter is set to 400Hz, which is kept within the 500Hz Nyquist frequency range to avoid aliasing of higher-frequency EMG and electromagnetic interference back to the 50Hz adjacent frequency band. For the framing of the single-channel real-time EEG acquisition sequence, the frame length is set to 200 milliseconds and the frame shift is set to 100 milliseconds: a single cycle of the 50Hz power supply frequency is 20 milliseconds, and 200 milliseconds is just enough to accommodate ten complete cycles. There are enough cycles within the frame for comparison, and the frame is short enough for the EEG rhythm to be approximately stable. The frame shift is half the frame length, so that any waveform segment falls within two adjacent frames. The processing result is updated every 100 milliseconds. Each frame of the EEG frame sequence is a 200-point voltage sampling segment. The sampling time interval is consistent with the acquisition end, which is 1 millisecond. A new frame is added every 100 milliseconds. The new frame shares the middle 100-millisecond sampling points with the previous frame, and the overlapping points are calculated independently in each frame. The input range of the single-channel real-time EEG acquisition sequence is ±375 millivolts. Under 24-bit quantization, the least significant bit corresponds to about 0.045 microvolts, which is three orders of magnitude lower than the power frequency coupling amplitude in the tens of microvolts range. A rectangular window is used for framing, not a Hanning window. Adding a window will suppress the samples near the frame edge. When the same power frequency cycle falls in the middle of the frame, the peaks and valleys are normal, but when it falls at the frame edge, the peaks and valleys are flattened, thus increasing the symmetry deviation rate. When a sampling interruption occurs in the single-channel real-time EEG acquisition sequence, the actual number of points in that frame is less than 200. The EEG frame sequence for that frame is delivered according to the actual number of points.

[0017] In some embodiments, the step of extracting power frequency stable segments from the EEG frame sequence to establish a real-time calibration baseline includes: verifying the intra-frame power frequency peak-valley symmetry deviation of the EEG frame sequence to obtain a power frequency cycle candidate set; performing cross-frame phase drift classification based on the power frequency cycle candidate set to form a power frequency stability classification; extracting the frame segment with the lowest phase jump variable using the power frequency stability classification to generate a power frequency stability benchmark set; and sparsely integrating the power frequency stability benchmark set after removing outlier frames to establish a real-time calibration baseline.

[0018] A candidate set of power frequency cycles is obtained by verifying the intra-frame power frequency peak-valley symmetry deviation of the EEG frame sequence. The following operations are performed independently within each frame of the EEG frame sequence: First, a 50Hz narrowband extraction is performed on the frame, with the passband center fixed at 50Hz and the bandwidth set to 5Hz, to suppress components that deviate far from 50Hz and prevent their zero-crossing points from being mistakenly identified as power frequency cycle boundaries. Then, the cycle boundaries are divided by the zero-crossing point of the rising edge of the bandpass waveform, with every 20 sampling points constituting a power frequency cycle. Within each cycle, the mean value of the cycle is subtracted, and the maximum positive value is taken as the positive peak amplitude, and the minimum negative value as the negative valley amplitude. The symmetry deviation rate is calculated using the formula δ=|V_p-V_v|÷((V_p+V_v)÷2), where V_p is the positive peak amplitude of the cycle, V_v is the absolute value of the negative valley amplitude of the cycle (both in microvolts), and δ is the symmetry deviation rate of the cycle, a dimensionless ratio that measures the relative degree of deviation between the positive and negative half-cycles. Power frequency coupled purely from the power supply line is symmetrical in both positive and negative directions. When electromyography, eye movement, or electrode shaking are superimposed, only one half-cycle is elevated or flattened. Therefore, δ does not need to estimate artifact energy cycle by cycle to identify contaminated cycles. Cycles with δ less than 0.08 pass the verification. This upper limit is calibrated based on the upper bound of the measured power frequency coupling asymmetry when the electrode contact impedance is in the normal range of 5 to 20 kΩ. The power frequency cycle candidate set is valid on a cycle-by-cycle basis. Within the same frame, only a portion of the cycles may pass the verification, or all ten cycles may pass. Short frames with insufficient points in the EEG frame sequence can yield fewer complete cycles, and the number of cycles obtained from the power frequency cycle candidate set is correspondingly reduced. The δ criterion is not relaxed due to short frames.

[0019] A power frequency stability classification is formed by classifying cross-frame phase drift based on a power frequency cycle candidate set. The power frequency stability classification categorizes frame segments according to their stability level: the highest-level frame segments show almost no shift in power frequency phase between adjacent frames, while the lowest-level frame segments show a continuous shift in the same direction in every frame. The classification does not consider the instantaneous phase reading of a single frame, but only whether the phase difference between adjacent frames within a frame segment, after being converted into a time value, continuously shifts in the same direction. Phase and amplitude statistics are taken from the cycles given by the power frequency cycle candidate set: for each cycle, the analytical phase and envelope amplitude of its 50Hz component at the midpoint of the cycle are taken, and the median is taken within the same frame as the representative power frequency phase and amplitude values ​​for that frame. The median is more tolerant of individual cycle skew than the mean. The difference in the representative phase values ​​between two adjacent frames is then converted into a time value, and the radian difference is converted to seconds at 50Hz to obtain how long the same power frequency waveform has shifted along the time axis between two frames. This quantity is easier to directly compare with the sampling interval than the raw phase difference. A frame segment consists of five consecutive frames, covering a waveform length of 600 milliseconds with half-frame overlap. The drift of a frame segment is calculated by dividing the absolute value of the algebraic sum of the reduced time values ​​of adjacent frames within the segment by the number of intervals. This calculation is only performed between adjacent frames that have passed the period verification in the power frequency cycle candidate set. If a frame has no period verification, this calculation between that frame and its adjacent frames is skipped, and the value is obtained from the remaining available frame intervals. The level difference for power frequency stability grading is 0.1 milliseconds: when the subject remains seated and the electrodes are in good contact, the reduced time value between frames usually fluctuates around 0.1 milliseconds. This value increases exponentially when the body moves slowly or the electrode edges begin to lift. 0.1 milliseconds is the boundary between the two states.

[0020] The lowest-level frame segment in the power frequency stability classification is used to generate a power frequency stability reference set. Even a frame segment that reaches the highest level in the power frequency stability classification may still contain a momentary phase jump. A subject biting their teeth or a lead wire being touched can cause a momentary shift in the power frequency phase, and the drift averaged over five frames barely reflects this jump. The phase jump variable is designed for this situation: it is the absolute value of the difference between the analytical phases of adjacent power frequency cycles within a frame segment, with the largest value being the phase jump variable for that frame segment. It measures a momentary jump rather than a slow directional shift. Extraction is performed on the two highest-level frame segments in the power frequency stability classification, calculating the phase jump variable frame by frame, and selecting the frame segment with the smallest value. The power frequency stability reference set is collected at a density of at least three frames per minute, determined by the period of power frequency amplitude step changes caused by power load switching. If the frame segment with the smallest jump variable within one minute still exceeds the 95th percentile of historical resting records, the frame segment length is reduced from five frames to three for recalculation. The recalculated and selected frames are marked as low confidence. In the power frequency stable reference set, frames marked with low confidence were removed from outlier frames at twice the median absolute deviation instead of three times, while the remaining frames were still removed at three times the median absolute deviation. Historical resting-state recordings were taken from a two-minute period after the subject had worn electrodes and before the formal data acquisition began, during which the subject sat still with eyes closed. This recording contained no conscious muscle activity and was used as baseline reference for subsequent thresholds. The upper limit for inclusion in the power frequency stable reference set was the 25th percentile of the phase jump variable in the historical resting-state recordings. Each included frame in the power frequency stable reference set simultaneously retained a 50Hz amplitude (µV) and a power frequency phase (radians), both describing the complete state of power frequency coupling at that moment.

[0021] After removing outlier frames from the stable power frequency reference set, sparse integration is used to establish a real-time calibration baseline. The median and median absolute deviation of the 50Hz amplitude for each frame in the stable power frequency reference set are calculated. Frames whose amplitude deviates from the median by more than three times the median absolute deviation are removed. These frames are often caused by sudden large movements of the subject or instantaneous electrode detachment from the skin, leading to abrupt changes in the power frequency coupling amplitude. Three times the median absolute deviation is approximately equivalent to two times the standard deviation under a normal distribution. The median absolute deviation is used here instead of the standard deviation because the frames to be removed would themselves raise the standard deviation. Sparse integration treats the center moment of the remaining frames as a node on the time axis. Linear interpolation is performed on the power frequency amplitude and phase between any two nodes. The values ​​from the last node onwards are retained. Integration is performed once per second. Before interpolation, unwinding is performed. When the phase difference between adjacent nodes exceeds π, the minimum difference is returned by adding or subtracting 2π to ensure that flips at the ±π boundary are not misinterpreted as true drift. In terms of node density, only three to six nodes per minute are typically needed to keep up with the step-like changes in power frequency coupling amplitude that last for more than ten seconds, eliminating the need for a reference for each frame. When the center times of two adjacent nodes are more than two seconds apart, the real-time calibration baseline maintains the value of the nearest node within this interval, without performing linear extrapolation across two seconds. The longer the interval, the more likely the power supply load may have changed direction, and the extrapolated slope is more likely to be opposite to the actual value. For each frame segment that passes the rejection process and is included in the power frequency stability reference set, the real-time calibration baseline only modifies the interpolation segment after that node, with a time resolution consistent with the acquisition end, both being 1 millisecond. Amplitude and phase are updated independently. If the amplitude of a node is refreshed but the phase has not yet been refreshed, the phase side still uses the previous interpolation result, and vice versa. Low-confidence markers in the power frequency stability reference set do not participate in the determination of three times the median absolute deviation during the outlier frame rejection phase.

[0022] Step S102: Based on the real-time calibration baseline, perform frame-by-frame baseline drift amplitude quantization to form a drift amplitude sequence. Based on the drift amplitude sequence, locate the drift abrupt change site to generate an EMG artifact candidate set. Mark the template morphology highly consistent segments in the EMG artifact candidate set to obtain the artifact morphology feature set.

[0023] In some embodiments, the step of performing frame-by-frame baseline drift amplitude quantization to form a drift amplitude sequence based on the real-time calibration baseline includes: performing frame-by-frame power frequency baseline stripping based on the real-time calibration baseline to form a power frequency-removed baseline frame set; identifying drift acceleration points with shortened reversal intervals based on the power frequency-removed baseline frame set to obtain drift direction records; identifying local anti-rise points within low-amplitude drift segments through the drift direction records to generate an anti-rise drift annotation set; and performing frame-by-frame amplitude quantization on the anti-rise drift annotation set to form a drift amplitude sequence.

[0024] Based on the real-time calibration baseline, a frame-by-frame power frequency baseline stripping process is performed to form a power frequency baseline-removed frame set. A sinusoidal waveform of the same frequency is reconstructed using the power frequency amplitude and phase given by the real-time calibration baseline at the corresponding time in each frame. This reconstructed waveform is subtracted from the original signal point by point. The reconstructed waveform and the original signal share the same clock reference at the acquisition end for time alignment. After subtraction, the power frequency coupling amplitude is significantly reduced, while fluctuations from other sources such as EEG rhythms, EMG, and eye movements are preserved intact in the power frequency baseline-removed frame set. The previously superimposed regular, fine jagged edges on the stripped waveform disappear, leaving only the slowly undulating baseline and its EEG rhythms. Stripping is performed independently within each frame, using the same framing method as the real-time calibration baseline: a 200-millisecond frame length and a 100-millisecond frame shift. Since the power frequency coupling itself fluctuates with the subject's posture and electrode contact state, the amplitude and phase used for reconstruction are taken from the corresponding instantaneous readings for each frame. Therefore, changes in power frequency amplitude caused by posture changes can be tracked in real-time by the stripping process. If a frame's time interval happens to cross an interpolation node, the reconstructed waveform is obtained point-by-point by taking the amplitude and phase of the interpolation curve at that point. The slight discontinuity caused by the change in slope at the node is preserved. No additional smoothing is needed here; preserving this slight step actually helps to accurately identify the direction of the frame crossing the node in subsequent drift direction recording. Each frame in the power frequency baseline removal frame set stores the complete waveform after stripping, and all subsequent direction determinations are performed on this stripped data. Before stripping, it is first checked whether there are valid nodes in the real-time calibration baseline within the frame's time interval; if a node is missing, the original waveform is directly filled into the power frequency baseline removal frame set. After the adjacent nodes are filled in, the frame is stripped again, and the new result overwrites the previously filled original waveform.

[0025] The drift direction is recorded by identifying drift acceleration points with shortened reversal intervals using the de-baseline frame set. The average waveform of each frame in the de-baseline frame set is taken as the baseline representative value for that frame, and the rising or falling direction of the representative values ​​of adjacent frames is recorded as one positive or negative unit. The determination of direction reversal only considers the sign of the rise or fall of the representative value itself, not the specific magnitude of the rise or fall. Even if the rise or fall amplitude varies greatly in some frames, as long as the sign direction is consistent, it is still considered that no reversal has occurred. The reversal interval is counted in frames. The number of frames between the same reversal and the previous reversal is the interval of this segment. During smooth drift, there is often only one reversal every few dozen frames. When the subject adjusts his or her posture or begins to sweat under the electrodes, the reversal interval will shorten by a factor of two within a few seconds, which is a precursor to the subsequent appearance of artifacts. The determination of shortened interval is based on the ratio of the average of the three most recent reversal intervals to the average of the three before that: if the ratio is less than 0.5, it is determined that the interval is shortened, and this point is recorded as a drift acceleration point. The boundary is derived from the normal fluctuation range near the median of the interval distribution of natural baseline drift. If the cumulative number of flips is less than six, the ratio determination is not performed temporarily. The calculation starts from the earliest six flips after the cumulative number of flips reaches six. A direction marker is written once at each drift acceleration point in the drift direction record. The marker only takes positive and negative values; amplitude information is not recorded here but is provided separately in the subsequent frame-by-frame amplitude quantization step. When the baseline always moves monotonically in the same direction and never flips, no marker is generated in the drift direction record during this interval, and the determination process continues frame-by-frame along the power frequency baseline frame set as usual. The markers in the drift direction record are arranged in the order of appearance of the drift acceleration points, with adjacent markers stored adjacent to each other. If no new spacing shortens in an intermediate frame segment, no marker is generated in this gap. Therefore, the actual number of frames between two adjacent markers may be much greater than one; a longer span indicates a more stable baseline between two accelerations.

[0026] Local reversal points within low-amplitude drift segments are identified by recording drift direction, generating a reversal drift marker set. The interval between two adjacent markers in the drift direction record is the segment to be examined. The direction that appears more frequently in each frame of the rising and falling symbols within this segment is taken as the dominant direction. If a single frame appears within the segment that is opposite to the dominant direction, and the two adjacent frames are still the dominant directions, this frame is identified as a local reversal point. This local reversal is different from the previously identified shortened reversal interval. The latter looks at the flipping rhythm itself, while the former looks at isolated reversals against the general trend. In an interval where the baseline is slowly decreasing in one direction, there may occasionally be one or two frames that briefly tilt upwards against this general trend, and then immediately fall back to the original direction. This is mostly due to momentary disturbances caused by slight electrode displacement or small movements of the subject. The definition of a low-amplitude drift segment is based on the difference between the baseline representative values ​​of the start and end frames of the segment taken from the de-power frequency baseline frame set. A low-amplitude segment is defined as one-third of the typical drift amplitude under resting conditions. This one-third limit is determined by the median level of the baseline change amplitude corresponding to the normal flip segment in the drift direction record. Significant changes exceeding this limit are not processed in this step; their amplitude information is retained in the de-power frequency baseline frame set and incorporated into the drift amplitude sequence in subsequent drift abrupt change point location steps. These isolated reverse shifts may not appear significant on their own, but they are often followed by a larger jump. Therefore, the reverse-rise drift annotation set records each local reverse-rise point along with its corresponding low-amplitude drift segment, leaving two frames before and after the point as an observation window. The width of this window is the upper limit of the number of frames the local reverse-rise itself typically lasts. Reverse fluctuations exceeding this upper limit are likely due to a shortened reversal interval already identified elsewhere in the drift direction record.

[0027] The anti-rise drift annotation set is subjected to frame-by-frame amplitude quantization to form a drift amplitude sequence. For each record in the anti-rise drift annotation set, the corresponding frame segment is first checked frame by frame to see if the baseline representative value is exactly the same as that of the adjacent frames. If the readings of adjacent frames are consistent and there is no change for three or more consecutive frames, it is judged as a zero-change flat frame and is removed from the corresponding interval of the anti-rise drift annotation set. Then, quantization is continued for the remaining frames. The frame-by-frame amplitude is calculated according to the formula D_k=B_k-B_(k-1), where B_k is the baseline representative value of the k-th frame, B_(k-1) is the baseline representative value of the previous frame, and D_k is the drift amplitude of the current frame relative to the previous frame. The unit is the same as that of the baseline representative value, which is microvolts. The value can be positive or negative, with positive and negative corresponding to rise and fall, respectively. Frame-by-frame amplitude quantization proceeds frame-by-frame along the entire set of frames from the power frequency baseline. In the anti-rise drift annotation set, only one anti-rise marker is added to the corresponding position of each observation window. Zero-change flat frames mostly correspond to signals repeatedly filled after brief acquisition interruptions; therefore, they are first removed before frame-by-frame quantization to avoid false zero values ​​lowering the statistical distribution of the entire drift amplitude sequence. Each position in the drift amplitude sequence corresponds to one frame, and the position stores the D_k calculated for that frame. The sequence is arranged in chronological order of frames. The D_k of the frame containing the local anti-rise point is usually opposite in sign to the preceding and following frames; this sign change is the trace left by the marker in the anti-rise drift annotation set on the drift amplitude sequence. For positions determined to be zero-change flat frames, their D_k is given by linear interpolation of the readings of the two preceding and following valid frames. Interpolation is only performed within a short window of no more than five frames after removal; long flat intervals exceeding five frames are considered acquisition anomalies and handled separately. The interpolated D_k is also written into the drift amplitude sequence, and an interpolation mark is added next to the value to distinguish it from the measured reading.

[0028] Based on the drift amplitude sequence, a candidate set of EMG artifacts is generated by locating abrupt change points. Taking each position on the drift amplitude sequence as the center, and including three positions before and after it (a total of seven points), the average absolute value of the six points excluding the center point is calculated. A point is identified as an abrupt change point when the absolute value of the center point exceeds four times this average. This four-fold multiple is determined by the ratio of the peak value to the mean of normal baseline fluctuations in historical resting records; normal fluctuations rarely exceed four times. EMG itself consists of short-term bursts with amplitudes far exceeding the baseline, and when superimposed on the baseline, this causes a sudden spike in the values ​​of the corresponding frames in the drift amplitude sequence, while the preceding and following frames remain at normal levels. Identifying abrupt change points involves selecting such isolated high values. Entirely sustained high-amplitude intervals are not considered in this step, as the drift in such intervals is already reflected in the amplitude trend of the drift amplitude sequence itself. For each abrupt change point, a corresponding entry is created in the candidate set of EMG artifacts, with the start and end times taken as two frames before and after that position, reserving complete waveform boundaries for the next step of template matching. Each record simultaneously records the ratio of the absolute value of the center point to the mean of the six points, called the trigger multiple. When the distance between two adjacent abrupt drift points is less than four frames, they are grouped into the same record. The start and end times are the union of their extended ranges, ensuring that a single EMG burst, even if it lasts for multiple frames, corresponds to only one record. If multiple consecutive frames within the coverage area of ​​a record exceed four times the mean, the start and end times cover all the frames exceeding the limit and are then extended by two frames. In this case, each record only contains the start and end times and the trigger multiple, without copying waveform data. The next step, template matching, directly returns to the time window corresponding to the power frequency baseline frame set. The role of the EMG artifact candidate set here is to define the time range to be matched.

[0029] In some embodiments, obtaining a set of artifact morphological features by labeling high-consistency segments of the template in the candidate set of electromyographic artifacts includes: performing candidate frame template matching quantization to generate a morphological deviation comparison set; identifying persistently low-value deviation frames in the morphological deviation comparison set to obtain high-consistency candidate segments; using the high-consistency candidate segments to filter out rhythmic morphologically consistent segments to form a pure artifact high-consistency set; and obtaining an artifact morphological feature set by categorizing the pure artifact high-consistency set according to the time span of occurrence.

[0030] A morphological deviation comparison set is generated by candidate frame template matching and quantization for the candidate set of EMG artifacts. The original waveform (taken from the baseline frame set) corresponding to each record in the candidate set of EMG artifacts is aligned with the template waveform by peak value. The two waveforms are then normalized to the same order of magnitude by their own peak values, and the difference is calculated point by point. The morphological deviation is calculated using the formula ε=(1÷N)×Σ(i=1 to N)|X_i-M_i|, where X_i is the normalized amplitude at point i of the record, M_i is the normalized amplitude at point i of the template waveform (both are dimensionless), N is the number of points resampled to the same length as the template, and ε is the morphological deviation of the record, which is the normalized dimensionless mean absolute difference. A smaller value indicates that the shape is closer to the template. The template waveform was derived from confirmed EMG recordings of the subject performing three to five teeth-clenching movements as instructed before the formal data acquisition began. Each recording was averaged point-by-point after peak alignment. The waveform exhibited a typical profile: a steep initial rise, a forward-shifting peak, and a slower decay than the rising phase. The length was fixed at 500 sampling points, consistent with the common span of an expanded single recording in the EMG artifact candidate set. Different contraction forces of the same movement resulted in different overall waveform heights. Normalization before calculating ε ensured that the deviation only reflected the degree of shape similarity, without incorporating differences in amplitude. When the recording length differed from the template, the template waveform was linearly resampled according to the actual recording length, changing only the time axis scale of the template waveform without altering its amplitude or shape. When the recording length differed from the template waveform length by more than 30%, the shape stretching caused by resampling became significant, and the ε value was marked as excessively stretched. The morphological deviation comparison set stored an overall ε value for each recording, with the frame order consistent with the frame order within the EMG artifact candidate set recordings. If a record spans a length far exceeding the typical length of a single electromyographic burst due to previous merging operations, the ε value is still calculated according to the formula and written into the morphological deviation comparison set. The determination of authenticity is then left to the next step based on the continuity of the ε value.

[0031] High-consistency candidate segments are obtained by identifying persistently low-value deviation frames in the morphological deviation comparison set. The template waveform slides frame by frame within each record of the morphological deviation comparison set, and the ε value corresponding to a window segment starting from that frame is calculated frame by frame. Adjacent records do not overlap. After labeling, a deviation trajectory is formed within each record arranged frame by frame. The gap frames between two adjacent records are not included in the continuity judgment. An ε value below 0.15 is considered a low value. This boundary is determined by the quartile of the deviation distribution of the template waveform and itself between different EMG burst samples. The deviation between normal EMG samples usually falls below this boundary. The larger shape deviation is mostly interference from shape resemblance rather than true EMG. Only frames with ε values ​​falling within the low range for three or more consecutive frames are considered high-consistency candidate segments. Singular low values ​​in one or two frames followed by a rise are not included. These sporadic low values ​​appear as isolated troughs on the trajectory given by the morphological deviation comparison set, often indicating coincidental local overlap of waveforms, and do not represent overall morphological consistency. The minimum span of three frames is determined by the shortest characteristic period of the template waveform itself. Windows shorter than three frames are insufficient to stably reproduce the template waveform outline on the morphological deviation comparison set, making inclusion in the high-consistency candidate segment uncertain. The start and end times of the high-consistency candidate segment are taken from the start and end times of the consecutive low-value frame segment itself, without additional extension. Extension has already been done in the previous step; here, only the range is defined according to the criterion itself. When a low-value frame segment crosses the boundary of two adjacent records, as long as the ε value remains low before and after the boundary, it is still considered the same high-consistency candidate segment. If an ε value rises above a low value once within a high-consistency candidate segment, it will break at the point of rise, even if it falls back to a low value in the next frame, and each segment will become an independent high-consistency candidate segment and proceed to the next step of judgment.

[0032] High-consistency candidate segments are used to filter out rhythm-morphologically consistent segments, forming a pure artifact high-consistency set. The original waveforms corresponding to the high-consistency candidate segments are first compared in frequency band energy. Those with a high-frequency energy ratio of more than 60% are retained and judged as real EMG; those with a ratio of less than 60% are judged as rhythm-morphologically consistent, that is, the fluctuations of this segment are composed of the EEG rhythm itself, and are filtered out from the high-consistency candidate segments. The high-frequency band is 20Hz to 400Hz, and the low-frequency band is 0.5Hz to 20Hz. The frequency grid containing 50Hz and its harmonics is first eliminated. The 20Hz boundary is determined by the frequency position where the EMG energy spectrum and the EEG paroxysmal rhythm energy spectrum first show a significant intersection in amplitude. There is a type of waveform that can also coincidentally resemble the EMG template in shape. The occasional K complex wave, spindle wave, or paroxysmal alpha rhythm onset when the subject is drowsy or switching attention is also steep. If all of these are included without discrimination, these real EEG rhythms will be mistakenly treated as artifacts. The 60% threshold is determined based on the high-frequency proportion distribution of confirmed EMG samples and confirmed paroxysmal EEG rhythm samples, ensuring that the two types of samples are distinct on either side of this threshold. The screening process is performed independently segment by segment. Within a high-consistency candidate segment, if the first half is EMG and the second half is rhythmically consistent, the true boundary is redrawn within this segment based on the high-frequency proportion. Only the rhythmic portion is removed, while the EMG portion remains in the pure artifact high-consistency set. Each segment in the pure artifact high-consistency set simultaneously records its start and end times and the high-frequency proportion value, with the time reference consistent with the high-consistency candidate segment. The high-frequency proportion equals the sum of the energy of all frequency grids above 20Hz divided by the total energy of that segment; it is a dimensionless ratio. The frequency grids are taken from the Fast Fourier Transform of the actual waveform of that segment. The closer the proportion is to 1, the more likely the segment is almost entirely composed of broadband high-frequency EMG components. A proportion just over 60% indicates that there is still some low-frequency energy residue in the EEG rhythm, which needs to be used to determine which frequency grids still have significant brainwave energy when constructing the frequency band overlap distribution map.

[0033] The artifact morphological feature set is obtained by categorizing the pure artifact high-consistency set according to the time span of occurrence. Based on the pure artifact high-consistency set, the start and end time difference of each segment is converted into frame numbers, which are used as the span value of that segment. The three-level boundaries are determined according to the ternary loci of the span distribution of existing artifact samples for the subject. The boundaries are updated monthly as samples accumulate, and no fixed frame number is used; the span itself carries clues to the cause of artifacts, and categorization preserves these clues. The artifact morphological feature set assigns each segment in the pure artifact high-consistency set to the corresponding level according to the categorization results. The start and end time, span frame number, and level of a segment together constitute the complete value of that segment in the artifact morphological feature set. Segments whose span falls exactly on the frame number of the two level boundaries are assigned to the level with the longer span. The boundaries may shift slightly up or down after being recalculated according to the sample quantiles, and the classification at the boundary will also be adjusted accordingly. The segments are categorized by frame span. Within each segment, the trigger multiples recorded in the candidate set of electromyographic artifacts are sorted from largest to smallest. A larger multiple indicates a more intense spurt and a higher requirement for subsequent notch depth. This sorting also determines the processing order of each segment when building the frequency band overlap distribution map. The frame span is converted to actual duration using a frame shift of 100 milliseconds. Short segments typically span three to four frames, corresponding to 400-500 milliseconds, close to the duration of a masseter muscle contraction; medium segments typically span five to eight frames, corresponding to 600-900 milliseconds, close to the duration of a frown or the first half of swallowing; long segments can reach more than ten frames, corresponding to around one and a half seconds, close to the duration of continuous speaking or a complete swallowing action. After the three segments are converted to duration, they generally correspond to the typical duration of common electromyographic artifact types. The artifact morphology feature sets are stored separately for these three segments, allowing subsequent frequency band overlap analysis to build independent frequency band overlap distribution maps based on artifact duration.

[0034] Step S103: Based on the artifact morphology feature set, implement the high overlap frequency band notch depth differential reduction output notch filter parameters, analyze the self-ringing tail span through the notch filter parameters to obtain the filtered contaminated section, and delineate the coexistence boundary between residual artifacts and brain waves according to the filtered contaminated section to obtain the residual protection zone.

[0035] In some embodiments, the step of implementing differentiated reduction of notch depth in high-overlap frequency bands based on the artifact morphological feature set to output notch filter parameters includes: identifying the overlap intervals between each artifact component and the brainwave frequency band based on the artifact morphological feature set to establish a frequency band overlap distribution map; using the frequency band overlap distribution map to suppress the upper limit of notch depth in a partition to form a restricted notch configuration; performing high-overlap frequency band notch center offset locking to generate a frequency-locked notch scheme for the restricted notch configuration; and verifying and analyzing the bandwidth boundary of the frequency-locked notch scheme by reducing depth band by band and outputting notch filter parameters.

[0036] Based on the artifact morphological feature set, an overlap distribution map is established by identifying the overlap intervals between each artifact component and the brainwave frequency band. For each segment of the artifact morphological feature set, the waveform corresponding to the baseline frame set is first subjected to a spectral analysis. The frequency band width is 1 Hz, covering 0.5 to 100 Hz, to obtain the energy of that segment in each frequency band grid. If the artifact energy and the reference energy of the subject's resting EEG baseline spectrum in the same grid both exceed one percent of their respective total energy, that grid is considered an overlap interval. Typical overlaps occur in the 8 to 30 Hz range, which covers major EEG components such as alpha and beta rhythms and falls within the frequency bands where facial muscle EEG energy is concentrated, such as that of the masseter and frontalis muscles. The resting EEG baseline spectrum is taken from historical resting records, and the energy distribution of each frequency band grid is recalculated weekly based on the latest closed-eye meditation segment. Grids with an energy percentage below one percent are considered noise floors. Grids containing 50 Hz and its harmonics were already removed in the previous step and are also not included in the judgment. The frequency band overlap distribution map is arranged by frequency band grid. Each grid is marked with whether it is an overlapping interval and two values: the artifact energy and the brainwave reference energy of the corresponding segment. Non-overlapping grids are only marked with whether they fall below the noise floor. The short, medium, and long segments in the artifact morphology feature set are analyzed independently in this step. The degree of overlap of the same frequency band grid may be different in segments of different levels. The frequency band overlap distribution map is stored separately for each segment, and the results of the three levels are not merged into a single map. The analysis window of a segment takes the number of sampling points corresponding to the actual span of the segment. If it is less than 1024 points, it is padded with zeros to 1024 points before transformation, so that the frequency band grid width is 1Hz. This avoids the frequency resolution loss caused by an excessively short window from obscuring frequency band grids that should be judged as overlapping. Zero padding only changes the frequency axis interpolation density and does not change the relative distribution of the original energy between the frequency band grids.

[0037] For example, the step of using the frequency band overlap distribution map to suppress the upper limit of the notch depth in a partition to form a restricted notch configuration includes: determining the real-time proportion of brain wave components in each overlapping interval based on the frequency band overlap distribution map to form a brain wave proportion distribution; identifying high-risk notch intervals with brain wave proportions exceeding a threshold through the brain wave proportion distribution to establish a high-risk notch list; performing depth difference analysis on each interval based on the high-risk notch list to obtain a partition depth constraint set; and integrating the partition depth constraint set with gradient widening to form a restricted notch configuration.

[0038] The brainwave proportion distribution is formed by determining the real-time proportion of brainwave components in each overlapping interval based on the frequency band overlap distribution map. The proportion is calculated using the formula R = E_b ÷ (E_b + E_a), where E_b is the brainwave reference energy of that frequency band grid after being converted according to the duration of that segment, E_a is the measured energy of the corresponding segment in the same frequency band grid minus E_b, both of which are energy units and are taken from the values ​​already stored in the frequency band overlap distribution map. R is the brainwave proportion of that grid, a dimensionless ratio between 0 and 1. The closer R is to 1, the more energy in that grid comes from the brainwave itself. If a deeper notch is applied here, it will suppress most of the brainwave rather than artifacts; the closer R is to 0, the opposite is true. The brainwave proportion distribution is calculated only for the overlapping interval grids marked on the frequency band overlap distribution map; non-overlapping grids are not included in this calculation step. The calculation proceeds segment by segment along the frequency band overlap distribution map. For the same frequency band cell, E_b has the same value across different segments, while E_a varies from segment to segment. Therefore, R is stored separately for each segment, and the R values ​​for the same frequency band cell may differ across segments. For individual frequency band cells, E_a is obtained by subtracting E_b from the measured energy. When the two values ​​are close, the estimate of E_a is easily affected by measurement fluctuations, deviating from the true value or even becoming negative. In this case, the calculated R result is no longer reliable. The cell is then recalculated using the average of E_a from the two adjacent frequency band cells, with the source indicated next to the recalculated result. E_b comes from the resting-state EEG baseline spectrum. The actual brainwave energy in the same frequency band cell may fluctuate under different test states. The resting-state EEG baseline spectrum is only a long-term averaged reference value; the R calculation does not attempt to track the true brainwave energy at a specific moment. The R values ​​in the brainwave proportion distribution retain the original calculation results for each cell, and the next threshold determination is performed based on these values.

[0039] A high-risk notch list is established by identifying high-risk notch intervals with brainwave proportions exceeding a threshold based on the brainwave proportion distribution. For each overlapping interval in the brainwave proportion distribution, the R value is compared to 0.4. Grids exceeding 0.4 are identified as high-risk notch intervals. The 0.4 threshold is determined based on the median level of the energy proportion of confirmed brainwave samples in each frequency band. Performing deep notch imaging beyond this level results in brainwave loss significantly outweighing the benefits of artifact suppression. The high-risk notch list is compiled for each high-risk notch interval, recording its segment affiliation, frequency band location, and corresponding R value. The R value is directly taken from the existing values ​​in the brainwave proportion distribution; this step is limited to numerical comparison and list registration. The 0.4 threshold is reviewed quarterly along with the resting-period EEG baseline spectrum. The review involves re-statistically analyzing the energy proportion distribution of confirmed brainwave samples in each frequency band. If the median level shifts significantly compared to the previous quarter, the threshold is adjusted accordingly. When multiple adjacent frequency bands within the same segment are classified as high-risk, the high-risk notch list merges adjacent cells into a single interval record, with the start and end frequencies taken from the boundaries of each adjacent cell. Merging only occurs when the R-value difference between adjacent cells does not exceed 0.1; if the difference exceeds 0.1, the interval is broken and a new interval is started. This prevents a single interval from being too wide, which could lead to subsequent depth constraints being treated with the same value, thus masking the inherent R-value differences within the interval. For segments where all frequency band cells have R-values ​​below 0.4, no new entries are added to the high-risk notch list, and the classification continues along the brainwave proportion distribution to the next segment. This situation often occurs in short, strong bursts where the artifact energy in each frequency band cell is much higher than the brainwave reference energy, and brainwave loss is almost negligible during suppression.

[0040] Based on the high-risk notch list, depth differentiation analysis is performed on each interval to obtain the partition depth constraint set. For example... Figure 2 As shown, each interval in the high-risk notch list is recorded on frequency axis 6 corresponding to high-risk notch interval 9. The maximum allowable notch depth is deduced by using the R value of this interval: the upper limit of depth A is given by multiplying the baseline depth A0 by (1 minus the R value), where A0 is... Figure 2The values ​​shown in reference depth line 8, where A and A0 are in decibels. The closer the R value is to 1, the closer the calculated upper limit of depth is to 0, meaning that notch filtering is almost impossible in this range, corresponding to the position of the notch depth curve 7 on frequency axis 6 that is closer to zero; the closer the R value is to the critical point of 0.4, the closer the upper limit of depth is to 60% of the reference depth, resulting in a difference between this and the reference depth used in non-high-risk ranges. The partition depth constraint set stores the upper limit of depth A item by item in the high-risk notch list. This reverse calculation relationship only holds true for ranges with R values ​​exceeding 0.4 and included in the high-risk notch list; ranges not included in the high-risk notch list do not undergo this reverse calculation. Non-high-risk overlapping and non-overlapping ranges uniformly use the reference depth A0 in the partition depth constraint set; if a segment has no high-risk ranges throughout, the partition depth constraint set gives the same reference depth for all frequency bands in that segment. If the calculated upper limit of depth for two adjacent intervals within the same segment differs by more than half of the baseline depth, the original calculation results are retained for each interval, and the difference is handled in the gradient integration step without any pre-smoothing. The baseline depth A0 adopts the value uniformly used in non-high-risk intervals of the subject's previous notch irradiation protocols, i.e., the value that appears most frequently in the most recently accepted protocols; A0 is set to 40 dB as the initial default value during the first acquisition, and is replaced by the measured statistical value after subsequent protocol acceptance. The upper limit of depth in the partition depth constraint set is shallower than the baseline, and the decimal value is retained for accurate interpolation during gradient integration.

[0041] A gradient-widening integration of the partition depth constraint set is used to form a restricted notch filter configuration. Progressing grid-by-grid along frequency axis 6, the depth difference between adjacent frequency band grids in the partition depth constraint set is limited to within 2 dB. When the difference exceeds 2 dB, the depths of each grid are rearranged from the deeper side outwards, decreasing by 2 dB per grid, until they connect with the shallower side. The rearranged grid depths are obtained through linear interpolation to prevent abrupt changes in the notch filter depth curve 7. Abrupt changes introduce additional sideband energy into the filter's frequency response, which can become a source of ringing. If the original step values ​​given by the partition depth constraint set are used directly before integration, the difference between adjacent grids will be retained in the notch filter parameter stage. Gradient integration aims to mitigate this potential problem from the design phase. The restricted notch configuration is stored in a grid-by-grid manner according to the final integrated depth value. Each grid is marked with the upper limit of the integrated depth and whether the grid has been rearranged. Rearranged grids are treated differently from grids directly taken from the partition depth constraint set when locking the notch center. Rearranged grids are generally not considered as candidate positions for the notch center. When the number of regular grids between two adjacent high-risk intervals is less than three, the two intervals and the regular grids in between are included in the gradient integration range and processed continuously according to the same transition curve. In this case, the restricted notch configuration will not have a clearly flat reference depth region in this continuous range. The whole segment has a gradual shape, which looks like a smooth transition arc rather than a step when viewed from the frequency axis. After integration, the final depth is checked grid-by-grid to see if it still does not exceed the upper limit given by the partition depth constraint set. If the overlap of the rearrangement ranges on both sides causes a grid to be assigned a value of the deeper side, the shallower value of the two rearrangement results is taken instead.

[0042] For restricted notch configurations, a frequency-locking notch center offset locking scheme is implemented to generate notch centers in high-overlap frequency bands. In a restricted notch configuration, frequency bands with an upper depth limit significantly lower than the baseline depth are identified as high-overlap frequency bands. These bands would normally fall on the notch center at the energy peak. Instead, the notch center is locked to the edge frequency of this overlap interval. The peak often coincides with the position where the brainwave proportion (R value) is highest, and artifacts and EEG rhythms almost completely overlap here. Locking to the edge avoids this position, and the tail of the artifact energy spreading to both sides falls precisely on the edge, thus being suppressed. The edge is taken from the side farther from the peak within the start and end frequencies of the overlap interval; when the distances on both sides are similar, the side with the lower upper depth limit is taken. The frequency-locking notch scheme assigns a locked center frequency to each high-overlap frequency band, while the notch center of non-high-overlap frequency bands remains at the original energy peak position. When the distance between the locked edge frequencies of two adjacent high-overlap frequency bands is less than 2Hz, the midpoint of their edge frequencies is taken as the unified locking center, and the lower of the two upper depth limits is taken. The locked center frequency, together with the upper limit of depth given by the restricted notch configuration, constitutes the complete value of each notch position in the frequency-locked notch scheme. The upper limit of depth is not readjusted at this step. The peak position is taken from the frequency band cell with the highest artifact energy E_a in the frequency band overlap distribution map. When taking the peak, the positions marked as rearranged cells by the restricted notch configuration are skipped. The depth value of the rearranged cells is linearly interpolated from the adjacent sides and does not correspond to any real energy reading. The offset notch can still suppress the tail of artifact energy spreading to both sides because the attenuation of the transition band of the notch has not yet returned to zero at the edge frequency, and it still has a certain suppressive ability on the frequency band cell immediately adjacent to the center. Frequency band cells in non-high overlap frequency bands where the artifact energy is significantly higher than the brain wave reference energy also need notching. Multiple consecutive cells are merged into one. These types of notches do not need to be offset from the center because the proportion of brain waves is low.

[0043] The notch filter parameters are output by reducing the depth band by frequency band according to the frequency-locked notch filter scheme and analyzing the bandwidth boundary. For each depth upper limit given by the frequency-locked notch filter scheme, it is checked whether it can be achieved using an integer-order notch filter. The depth increment brought by each increase in order is approximately equal within the design bandwidth. If the depth upper limit, after being converted to order, falls between two integers, the depth is reduced and rounded down, and a lower integer order is used to ensure that the actual depth does not exceed the upper limit. The difference between the actual depth after reduction and the original upper limit is the depth loss L, in decibels. L is written into the verification record for each frequency band, and subsequent ringing judgment and residual protection zone delineation are based on this input. The bandwidth boundary is calculated by inversely using both the order and the target depth: the higher the order and the deeper the target depth, the wider the transition band. The actual bandwidth is taken as the frequency span corresponding to the attenuation of the transition band to half the reference depth at that order. The notch filter parameters are given for each notch location in the frequency-locked notch scheme, including the center frequency, the actual order after depth reduction, and the bandwidth. These three parameters together constitute the complete design of the notch in the final filter. The center frequency is directly adopted from the value given in the frequency-locked notch scheme; the bandwidth, which was not present in the frequency-locked notch scheme before, is a newly resolved value in this step. When the bandwidth boundaries of two adjacent notch locations overlap, the overlapping part is based on the boundary of the deeper side, and the shallower side is correspondingly shrunken inward to avoid the same frequency being repeatedly included in excessive attenuation. Notches with L exceeding 1 dB are additionally marked as having excessive depth loss, and this marking will be used in subsequent ringing anomaly judgment and residual artifact ratio estimation. For adjacent frequency bands that have been determined to be notch locations by the restricted notch configuration, their actual bandwidth is given at once based on the overall coverage of multiple consecutive notches, and the corresponding notches share the same bandwidth value in the notch filter parameters.

[0044] In some embodiments, obtaining the filtered contaminated section by resolving the self-ringing trailing span through the notch filter parameters includes: calculating the time for each order impulse response to decay to the background using the notch filter parameters to obtain an impulse response duration set; identifying notch filters whose decaying duration exceeds the theoretical duration based on the impulse response duration set to form a ringing over-limit group; measuring the over-limit ringing extension span for the ringing over-limit group order by order to generate a trailing extension map; and marking the trailing extension map with ringing contaminated segments outside the boundary to obtain the filtered contaminated section.

[0045] The impulse response duration set is obtained by calculating the decay time to the background of each order of impulse response using notch filter parameters. For each notch filter parameter, filter coefficients are constructed based on its center frequency, actual order, and bandwidth. A unit impulse is input, and the decay process of the output amplitude over time is recorded point by point until the output amplitude falls below the background noise level and remains below it for twenty consecutive sampling points. The duration corresponding to this moment is recorded. Brief low-value fluctuations occurring before the output amplitude falls back to the background noise level are not considered part of the decay completion. The background noise level is determined by the measured root mean square value of the clean segment uncontaminated by power frequency and electromyography in the subject's historical resting-state recordings. This value is updated every minute with the latest resting-state segment to avoid using outdated background readings for extended periods. The impulse response duration set is stored sequentially according to the notch filter parameters, with the decay time T_m being the notch number in milliseconds. In addition to calculating the actual order for each notch, the attenuation time for each of the two adjacent orders above and below it is also calculated once, resulting in the attenuation time for the three orders at the same center frequency. This facilitates comparison within the group in the next step. The order and duration must be calculated individually to obtain precise values; a fixed proportional relationship cannot be assumed. If two adjacent notches share the same bandwidth value in the depth-of-dip verification, the calculation is still performed separately for their independent center frequency and order. If the bandwidth of a notch is significantly wider than the conventional value due to the merging of adjacent notches, the calculation window is extended to twice the conventional window to ensure that the attenuation tail is not truncated by an excessively short observation window after the bandwidth widens. Each T_m in the impulse response time set is also marked with whether the notch uses a conventional window or an extended window, so that the corresponding reference benchmark can be selected in the next step of over-limit judgment.

[0046] Notch filters with attenuation times exceeding theoretical durations are identified based on the impulse response duration set, forming a ringing-out group. Notch filters of different orders at the same center frequency are grouped together in the impulse response duration set. Within each group, the attenuation durations are arranged from lowest to highest order; under normal circumstances, the duration should monotonically increase with increasing order. The measured durations of lower orders within the group are extrapolated to higher orders using the conventional engineering relationship between order and bandwidth, yielding the expected duration of the higher-order notch filter calculated from the lower-order filters of the same class. If the measured duration of a higher-order filter exceeds this expected duration by more than double, it is determined to be a ringing-out group. The extrapolated baseline value is then adjusted downwards using the depth loss recorded in the notch filter parameter verification record. A larger depth loss indicates that the actual notch is shallower than the design depth, storing less energy, and should attenuate faster. The baseline value is adjusted downwards by 5% for every 1 dB of depth loss L, with the adjustment not exceeding half the baseline value. The adjusted value is denoted as T_ref, which is the theoretically expected duration calculated from similar low-order filters used in this comparison step. Notch filters with extended windows are marked in the impulse response duration set. The measured duration has already captured the tail using a double window, so no further scaling is needed before comparison. If such notch filters also fall into the ringing exceedance group, it indicates that the exceedance is not an illusion caused by an insufficient observation window, making the judgment more reliable. The ringing exceedance group summarizes the parameters of each exceedance group into an entry, writing the group's center frequency, the range of orders involved, and the exceedance factor. The exceedance factor is given as T_m divided by T_ref minus one; a larger factor indicates that the ringing deviates further from the expected value. If notch filter parameters generated at the same center frequency in different time periods all exceed the limits, they should be recorded separately. Each record should also indicate the original entry position of the corresponding notch filter parameter.

[0047] For the ringing exceedance group, the excessive ringing extension span is measured order by order to generate a trailing extension map. Within each record of the ringing exceedance group, the number of orders is progressively increased. The extension span of each order is calculated using the formula ΔT = T_m - T_ref, where T_m is the measured attenuation time of that order in the impulse response time set, T_ref is the theoretical expected time after extrapolation from lower orders of the same order and correction by depth loss, and ΔT is the time span for the ringing to extend to the theoretical value for that order. All three are in milliseconds, and ΔT is always positive. Different orders within the same group calculate their own ΔT, and the average of the entire group is not used to represent the whole group, because the causes of ringing at different orders within the group may not be the same, and using the average would mix together deviations from different causes. The trailing extension plot arranges the ΔT values ​​of each order along two dimensions: center frequency and order. The trend of ΔT changing with order at the same center frequency—whether it continuously increases or suddenly jumps at a certain order—is also reflected on the plot. When determining the contamination range, this trend is used as a reference, not just a single point value. For some orders, because the filter poles are closer to the unit circle, the accuracy loss of coefficient values ​​increases, and the attenuation rate is significantly slower than adjacent orders in the same group. The jump amplitude can be more than twice the average extension span of the group. The order with the largest ΔT value often corresponds to the location with the largest depth loss L in the verification record. When delineating the filter contamination section, the ΔT of this order is taken as the contamination range at that location. Other orders are only used as references for comparison and are not demarcated separately. Although the ΔT values ​​of other orders are not involved in demarcation, if they show an abnormal trend of jumping with order on the trailing extension plot, they can still serve as clues for ringing cause analysis. The T_ref should be of the same origin and have the same value as the theory used when identifying ringing exceeding limits, and can be directly adopted.

[0048] The filtered contamination section is obtained by labeling the ringing contamination section outside the boundary of the trailing extension map. For example... Figure 2As shown, each ΔT in the tail extension diagram is first converted to the actual time axis. The starting time 2 of the notch application in the scalp EEG signal waveform 1, plus the theoretical duration of that order, is the theoretical duration boundary 4 of the ringing attenuation tail segment 3 that should converge. The interval extending ΔT after boundary 4 is the measured tail contamination segment 5. The waveform in this interval is the residual vibration of the filter itself before it has completely decayed. Its shape is very similar to a regular rhythm with gradually weakening amplitude, which is easily mistaken for real EEG. The correspondence between the high-risk notch interval 9 on the frequency axis 6 and the ringing attenuation tail segment 3 generated on the time axis is marked by the correspondence indicator line 10. When multiple orders in the same center frequency and the same time period exceed the limit, the interval marked by the order with the largest ΔT is taken as the contamination range at this point. Each order is not marked separately to avoid overlapping counting, because the filter finally applied to the signal has only one actual order, and the other orders are intermediate comparison objects in the calibration process. The filtered contaminated sections are listed one by one according to their contamination locations, with start and end times and corresponding center frequencies recorded. The start and end time references are consistent with the actual signal time axis where the trailing extension diagram is located. For each notch, the start time 2 is taken as the starting point of the artifact segment to which the notch belongs, and it is from the same source as the clock reference at the acquisition end. When the time spans of two adjacent contaminated sections overlap, they are included in the same record, and the start and end times are taken as the union of the two, while their respective center frequencies are retained. When delineating residual protection zones, each center frequency under this record is checked back and verified, and the one with the largest depth loss is taken as the reference. When ΔT is extremely small, such that the calculated contaminated section is less than a sampling point, the filtered contaminated section is still recorded according to the width of the minimum sampling point.

[0049] Based on the delineation of the residual artifact and brainwave coexistence boundary within the filtered contamination section, a residual protection zone is obtained. For each record within the filtered contamination section, the upper limit of depth A is retrieved by referring back to the frequency-locked notch filter scheme at the center frequency, and the depth loss L is retrieved by referring back to the notch filter parameter verification record. The larger the depth loss, the less thoroughly the artifacts are suppressed within that interval, and the higher the proportion of residual artifacts. The proportion of residual artifacts is calculated using the formula W = 10^(-(A-L)÷20), where A is the upper limit of depth at that location, L is the depth loss at that location (both in decibels), and W is the proportion of residual artifacts, a dimensionless ratio between 0 and 1. The W given by the formula is the peak proportion of residual artifacts at the instant of notch filter application. Within the filtered contamination section, the amplitude of residual artifacts attenuates at each sampling point along the ringing tail. Intervals with W exceeding 30% are classified as residual artifact-dominated, where residual artifact fluctuations in the waveform are still visually identifiable. Intervals with less than 30% are classified as coexisting brain waves and residual artifacts, but with brain waves dominating; in this case, the residual artifacts are submerged within the normal amplitude variations of the EEG rhythm. The residual protection partition is registered line by line for each defined interval, recording the start and end times, dominant category, and corresponding W value. When the boundary between the two categories falls within a record in the filtered contamination section, the boundary position is taken as the time point corresponding to the ringing amplitude envelope within that record dropping to 30% of the initial amplitude, without rounding to the start and end times of the interval. Intervals classified as brain wave-dominated and residual artifact-dominated within the residual protection partition retain their own start and end times and W values; even if the two types of intervals are temporally adjacent, they are not merged into the same record. The role of brain wave-dominated intervals in this residual protection partition is to indicate that this interval is not purely residual artifact; normal brain wave components dominate here, and it should not be treated as a complete artifact for repair.

[0050] Step S104: Delineate the boundaries of safe repair segments according to the degree of separability of the residual protection partition to construct a repair candidate sequence. Perform rhythmic phase anchoring completion on the repair candidate sequence to form an interpolation repair sequence. Use the interpolation repair sequence to determine the suspicious silent segment based on the segment with the lowest repair rate.

[0051] Specifically, the boundaries of safe repair segments are defined according to the degree of separability of the residual protection zone to construct a candidate repair sequence. The degree of separability is determined at each boundary of the residual protection zone. Returning to the attenuation envelope of the corresponding notch impulse response, the rate of decrease of the ringing amplitude near the boundary time point is examined: if the ringing attenuation rate is fast, and the amplitude drops significantly below 30% within two or three sampling points, the boundary is considered easily separable. This is often seen in short ringing caused by high-order narrowband notch filters. If the attenuation is slow, and the amplitude fluctuates around 30% for more than ten sampling points before barely crossing the boundary, the boundary is considered difficult to separate. This is often seen in long-tailed ringing caused by low-order broadband notch filters. The boundaries of safe repair segments are directly taken from the original start and end times given by the residual protection zone at easily separable boundaries; at difficult-to-separate boundaries, the boundary is contracted inwards by half the span of the fluctuation, so that the repair segment boundary avoids the region where the reading itself is volatile. Each interval identified as dominated by residual artifacts in the residual protection partition is extracted and segmented for repair candidate sequences. The boundary of the shrunken, safe repair segment and the separability marker at that location are then written into the sequence. After determining the boundary of the safe repair segment, the energy proportion of the dominant rhythm frequency band is verified on the waveform inside the boundary. Only segments confirmed to be dominated by artifact residues and with insignificant normal EEG rhythm characteristics are ultimately included in the repair candidate sequence. When the distance between the boundaries of two adjacent residual artifact-dominated intervals after shrinkage is less than three sampling points, they are included in the same repair candidate sequence record, with the interval portion merged and processed uniformly as a whole. The separability marker is the one that is more difficult to separate, ensuring a more cautious strategy in subsequent processing. The upper limit of the shrinkage is set to one-third of the original span of the interval; even when the tug-of-war span exceeds this ratio, only one-third is included. When the shrunken span is extremely short, less than twenty sampling points, and the anchor points at both ends cannot accommodate at least one complete rhythm cycle, the filtered waveform of this interval is output as is, equivalent to not interpolating this extremely short residual.

[0052] In some embodiments, performing rhythmic phase anchoring and completion to form an interpolated repair sequence for the repair candidate sequence includes: identifying complete rhythmic cycle boundaries before and after each segment of the repair candidate sequence to obtain a rhythmic anchoring point set; planning a phase continuity constraint path across defect segments based on the rhythmic anchoring point set to generate a phase transition scheme; performing frame-by-frame defect filling using the phase transition scheme to obtain a phase-filled frame set; and embedding the phase-filled frame set after performing over-amplitude filling frame downgrading annotation to form an interpolated repair sequence.

[0053] A set of rhythm anchor points is obtained by identifying the boundaries of complete rhythmic cycles before and after each segment of the candidate repair sequence. The dominant rhythm frequency band is the frequency corresponding to the highest energy band in the subject's resting EEG baseline spectrum. For each record of the candidate repair sequence, an observation window is taken before and after the boundary of the safe repair segment. The window length is three cycle lengths corresponding to the peak frequency of the frequency band. The zero-crossing point of the rhythm frequency band component is found within the window. The interval between two adjacent zero-crossing points in the same direction constitutes a complete cycle. The end point of the complete cycle closest to the boundary before the boundary and the start point of the complete cycle closest to the boundary after the boundary are each taken as an anchor point. Thus, a candidate repair sequence record has one anchor point before and after the boundary, for a total of two anchor points. The interval between the two anchor points is the entire safe repair segment marked by this record. The anchor point itself does not fall inside the segment, but only adheres to the two sides of the segment. The rhythm anchor point set stores a pair of anchor points one by one according to the recording order of the candidate repair sequences. Each anchor point records its time position and phase reading. The number of anchor points in the rhythm anchor point set maintains a strict two-to-one correspondence with the number of records in the candidate repair sequence, and this correspondence remains unchanged regardless of the observation window's extension or substitution. If the number of zero-crossing points of the dominant rhythm frequency band component within the observation window is less than three, the window is extended to five cycles before a new search is conducted. If the number is still less than three, the anchor point is replaced by the sampling point within the observation window whose amplitude is closest to zero and whose slope sign is consistent with the adjacent zero-crossing point. The replacement point is marked with an approximation label to distinguish it from the directly found zero-crossing point. The phase reading of the anchor point is taken as the analytical phase of the dominant rhythm frequency band component at that point, and uniformly converted to the starting phase of the corresponding cycle, serving as the reference zero point for subsequent phase transition planning. The analytical phase is obtained through Hilbert transform, which can provide an instantaneous phase estimate within a window of less than one cycle, suitable for situations where the anchor point falls precisely at an asymmetrical position in the rhythm waveform.

[0054] For example, the step of planning a phase continuity constraint path across a defect segment and generating a phase transition scheme based on the rhythm anchor point set includes: verifying the period deviation of adjacent anchor points in the rhythm anchor point set to obtain a period duration deviation set; identifying abnormal anchor point segments with abrupt changes in period duration based on the period duration deviation set to obtain an abnormal rhythm boundary set; enumerating and calculating the stability ratio of each candidate path using the abnormal rhythm boundary set to establish a stable path classification record; and performing a minimum phase difference fluctuation path selection on the stable path classification record to generate a phase transition scheme.

[0055] The cycle duration deviation set is obtained by verifying the cycle deviation between adjacent anchor points in the rhythm anchor point set. For each pair of anchor points, three complete cycles are reviewed outwards along the rhythm anchor point set, and the duration of each of these three cycles is recorded. The cycle duration deviation is calculated using the formula Δp_i = p_i - p_med, where p_i is the measured duration of the i-th cycle reviewed, p_med is the median cycle duration of the first three minutes of resting time after the start of the current recording for that subject, and Δp_i is the deviation of that cycle from the median level. All three values ​​are in milliseconds; a positive value indicates that the cycle is longer than normal, and a negative value indicates that it is shorter. i ranges from 1 to 3, corresponding to the first, second, and third review cycles from near the boundary to far from the boundary. p_med is then recalculated every ten minutes based on the latest resting time, without using statistical values ​​from other subjects. The cycle duration deviation set consists of three Δp_i values ​​for each pair of anchor points in the rhythm anchor point set, arranged in review order from near the boundary to far from the boundary. If the three Δp_i values ​​consistently move in the same direction and their absolute values ​​gradually increase, with the periodic deviation increasing closer to the boundary, it indicates that the rhythm has begun to distort as it approaches the defective segment. If the three Δp_i values ​​alternate between positive and negative without a single directional trend, it indicates that this rhythm is within the normal range of periodic fluctuations. Both are accurately recorded in the periodic duration deviation set, and their distribution differences are used as the initial basis for stability grading in subsequent path enumeration. When the anchor point in the rhythm anchor point set is one of those selected close to the approximate substitution point, the three cycles reviewed are counted outward from the approximate substitution point. The first cycle itself carries an approximation error, so an additional approximation mark is added during Δp_1 calculation. This mark is used to separately distinguish the confidence level in the next step of identifying abrupt changes.

[0056] Abnormal rhythm boundary sets are obtained by identifying anomalous anchor point segments with abrupt changes in cycle duration based on the cycle duration deviation set. Each pair of anchor points has three Δp_i in the cycle duration deviation set, progressing sequentially: if the absolute value of the difference between two adjacent pairs exceeds 20 milliseconds, it is determined to be a cycle duration abrupt change; the 20-millisecond limit is also statistically calculated based on the same historical resting record, taken from the same batch of cycles as p_med, and recalculated together with p_med. The differences between the first and second Δp_i and the differences between the second and third Δp_i in the cycle duration deviation set are compared separately. If either of them exceeds 20 milliseconds, a judgment is triggered. It is not required that both exceed the limit simultaneously, because abrupt changes often only occur in the one or two cycles immediately adjacent to the boundary, while the distant cycles usually maintain normal fluctuations. The abnormal rhythm boundary set independently creates a bar for each anchor point segment determined to be abrupt change, recording the record number of the repair candidate sequence to which the segment belongs, the pair of adjacent cycles in which the abrupt change occurs, and the specific difference and sign of the abrupt change. A positive sign corresponds to a sudden increase in period, indicating a brief deceleration of the rhythm; a negative sign corresponds to a sudden decrease in period, indicating a brief acceleration of the rhythm. The next step is to determine whether the abrupt change and the path expansion / contraction direction are in the same direction. If they are in the same direction, it means that the period jump of the fill waveform at that point will be more drastic, and the expansion / contraction allocation needs to be more careful. The two sides before and after a pair of anchor points are judged independently, and the abrupt change on the front side does not affect the judgment on the back side. If abrupt changes occur on both sides, the aberrant rhythm boundary sets are set on the two sides before and after the same point, and the difference and positive / negative sign are recorded separately. Δp_1 with an approximate label participates in the difference calculation in this step, but the abrupt change judgment calculated by it is additionally marked with a low confidence label to distinguish it from the abrupt change calculated by directly measuring the period difference.

[0057] A stable path grading system is established by enumerating and calculating the stability percentage of each candidate path using an abnormal rhythm boundary set. A candidate path is a scheme that fills the defective segment between the pre-boundary anchor point and the post-boundary anchor point given by the rhythm anchor point set with a complete EEG rhythm. This is analogous to filling a fixed-length gap with two to five bricks of varying lengths. One path is enumerated for each of the two to five cycles; the more cycles, the smaller the required expansion / contraction amplitude per cycle. The percentage of each candidate path is obtained by dividing the number of cycles through which the expansion / contraction direction remains consistent with the previous cycle by the total number of cycles minus one. This is a dimensionless ratio; a higher percentage indicates a smoother path. If a segment marked as a mutation in the abnormal rhythm boundary set falls within the period range traversed by a candidate path, the proportion of this path is additionally deducted. The deduction amount is calculated by dividing the absolute value of the mutation difference recorded in the abnormal rhythm boundary set by p_med and converting it into a dimensionless ratio. When the mutation direction is the same as the scaling direction of the path itself in this segment, the deduction amount is multiplied by 1.5, because such unidirectional mutations will overlap with the scaling of the path itself, making the actual jump more dramatic than the difference number appears. No additional bonus is added when the directions are opposite. If there are multiple mutations within the range of a path, the deduction amounts for each mutation are calculated separately and then added together. The largest deduction is not taken to represent the whole; the net proportion after the sum is the final value registered in the stable path classification record. The stable path classification record arranges each candidate path recorded under each repair candidate sequence in descending order of proportion. Even if the candidate path with the lowest proportion is less than 30%, it is still recorded. The classification here only serves a sorting function.

[0058] A phase transition scheme is generated by selecting the path with the least phase difference fluctuation from the stable path classification records. The phase increment of each step of a path refers to the cumulative phase traversed within that step cycle. The two candidate paths with the highest proportion in the stable path classification records are selected. The difference between the phase increment of each step of the path and the phase increment of its adjacent step is calculated cycle by cycle. The standard deviation of the difference is the phase difference fluctuation of the path, and the unit is the same as the phase increment, both in radians. If the phase difference fluctuation of the path with the highest proportion in the stable path classification records is significantly greater than that of the second most common path (more than double), the second most common path is selected. This is because the proportion only reflects whether the direction is consistent, not whether the expansion and contraction amplitude of each step is uniform. Although the proportion of paths with large fluctuations is high, the expansion and contraction of each step are uneven, and a visible waveform jump is easily left at the junction. The phase transition scheme provides a final selected path for each repair candidate sequence record, and writes the target phase increment for each cycle traversed by the path. The increment is evenly distributed along the path according to the phase difference between the start and end points of the cycle to each sampling point within the cycle. The next frame-by-frame filling is based on this value. When the phase difference fluctuations of the two most prevalent paths differ by less than double, the path with the higher prevalence is taken. In this case, the phase transition scheme provides the complete value for the path with the highest prevalence. The number of cycles for the path with the highest prevalence differs from that of the path with the second highest prevalence. The standard deviation of the phase difference fluctuations is calculated separately for each path based on its actual number of cycles, without forcibly aligning them to the same number of cycles before comparison. When both paths have fewer than three cycles, the maximum value of their single-step phase increments is directly compared, and the path with the smaller maximum value is taken. The fewer the number of cycles, the larger the average phase difference that needs to be crossed per step. This comparison method is more intuitive than standard deviation when the number of cycles is small.

[0059] A phase-filled frame set was obtained by performing frame-by-frame defect filling of the rhythmic phase reference using a phase transition scheme. The target phase increment given by the phase transition scheme is advanced point by point. The phase of each point to be filled is the starting phase of the path plus the cumulative increment. The amplitude is converted according to the value of the typical waveform template of the dominant rhythm frequency band of the subject at the corresponding phase. The template waveform is taken from the average waveform of multiple complete cycle waveforms of this frequency band in historical resting records after phase alignment. It is uniformly resampled into 360 interpolation points at 0 to 2π radians, with each point corresponding to one degree of phase resolution. The value is obtained by linear interpolation according to the target phase. The trough phase reference refers to the phase corresponding to the lowest point of the template waveform amplitude. The phase scale of the template waveform is calibrated with the trough as the zero point. During filling, the anchor point phase is converted according to the same scale and used as the starting point before accumulating the increment. Therefore, all repair candidate sequence records share the same template waveform reference, but the starting phase is still different from the anchor point. The starting phase of the path originates from the rhythm anchor point set. The trough calibration step only standardizes the phase scale of the template waveform; the number of cycles traversed by the path and the target phase increment for each cycle remain unchanged. Each candidate repair sequence record does not share the same cumulative increment and must be filled independently. The phase-filling frame set stores a filling amplitude for each sampling point within the safe repair segment boundary of each record in the candidate repair sequence, arranged in the order of sampling time. The filling amplitude itself is also superimposed with a random perturbation of the subject's current local background noise level, the perturbation amplitude being taken from the root mean square noise value of a segment of real signal immediately outside the safe repair segment boundary. Each filling point in the phase-filling frame set is also labeled with its corresponding path cycle number, indicating which transition cycle in the phase transition scheme it belongs to. This is used to locate the low-quality transition segment during subsequent overamplitude determination.

[0060] After downgrading and annotating the phase-filled frame set, an interpolation repair sequence is formed. For each fill point in the phase-filled frame set, if the difference between the amplitude of the fill point and the corresponding amplitude of the same phase on the template waveform exceeds 40% of the peak-to-peak value of the template waveform, it is determined to be an over-amplitude fill frame. This 40% threshold is determined by the upper bound of the normal fluctuation range of the true rhythm waveform relative to its own template waveform mean in historical resting records. If three or more consecutive sampling points under the same path period number are judged as over-amplitude, the entire transition segment corresponding to this period number is marked as a low-quality transition, allowing for tracing back to which specific period in the phase transition scheme the target phase increment was improperly calculated. Downgrading annotation does not change the amplitude of the fill point itself; it only adds a quality marker next to the point. Over-amplitude points are marked as low-quality, and non-over-amplitude points are marked as normal. If the same sampling point triggers annotation due to both amplitude exceeding the limit and consecutive adjacent over-limit frames, only one low-quality marker is recorded, without repeated superposition. The interpolation repair sequence embeds all fill points in the phase-filled frame set, along with their respective quality markers, into the time axis of the original signal according to the corresponding safe repair segment boundaries. The original signal outside the boundaries is preserved as is, without any modification. Each record simultaneously stores the original waveform reading of the segment before embedding, one copy for the fill result and one copy for subsequent silent feature verification, which is performed directly on this original data. The number of overamplitude sampling points and the total number of fill points recorded in each repair candidate sequence are written into the interpolation repair sequence along with the embedding action, and the repair rate is calculated from these two values. When covering multiple low-quality transition cycles, the number of overamplitude sampling points in each cycle is first summed and then divided uniformly by the total number of fill points.

[0061] Suspicious silent segments are identified using the segment with the lowest repair rate in the interpolation repair sequence. The repair rate is calculated using the formula η = 1 - n_low ÷ n_total, where n_low is the number of fill sampling points marked as low quality in a record of the interpolation repair sequence, n_total is the total number of fill points in that record (both in units of 1), and η is the repair rate of that record, a dimensionless ratio between 0 and 1; the lower the η, the less reliable the repair. All records in the interpolation repair sequence are arranged from low to high according to η. The record with the lowest η is taken as the benchmark. The η of the remaining records is compared with this benchmark. Records with a difference of no more than 0.05 are included in the judgment range and are considered to belong to the category of abnormally low repair rates, along with the benchmark. Within the judgment range, the original waveform readings before filling are retained alongside the interpolated repair sequence for each record. The amplitude is checked point by point to see if it consistently matches the background noise level of the acquisition channel. If the similarity exceeds 80% of the sampling points, the record is marked as a suspicious silent segment. For example, the gradual drying of conductive paste or the lifting of electrode edges causing a sudden increase in impedance likely points to a genuine signal loss caused by poor electrode contact or lead detachment. The start and end times and corresponding durations of the safe repair segment are written into each record in the suspicious silent segment. Records with a similarity of less than 80%, although having a low repair rate, still exhibit waveform fluctuations and are not classified as suspicious silent segments. These often occur during the period when the subject transitions from wakefulness to drowsiness, or when the dominant rhythm itself fluctuates in strength, making it difficult to find a smooth transition in path planning, which differs from the nature of signal loss at the hardware level. The degree of proximity is given by dividing the number of sampling points falling within a range of plus or minus twice the background noise level by the total number of sampling points in that segment. This is a dimensionless ratio. After a sudden increase in impedance, the amplitude of the residual signal is usually compressed to near the background noise level and loses its periodic fluctuations. This is quite different from the amplitude distribution of normal EEG rhythms. This difference in distribution is the key to distinguishing between true dropout and rhythm attenuation.

[0062] Step S105: Based on the analysis of suspected silent segments, output signal quality rating for the risk of detachment and coverage, implement downgrading labeling of difficult-to-separate overlapping segments to establish a graded EEG sequence, and perform explicit gap labeling of unsuppressed segments to output a suppressed EEG sequence according to the graded EEG sequence.

[0063] Specifically, the quality rating of the output signal is based on the analysis of suspected silent segments to assess the risk of electrode detachment and poor contact. The detachment coverage rate is calculated using the formula ρ=Σ(i=1 to n)t_i÷T_rec, where t_i is the duration of the i-th record in the suspected silent segment, n is the number of records, T_rec is the total duration of this record, and both t_i and T_rec are in seconds. ρ is the detachment coverage rate, a dimensionless ratio between 0 and 1, which measures the proportion of the entire record segment judged as suspected electrode detachment or poor contact coverage. A ρ below 5% is rated as excellent, between 5% and 20% is rated as acceptable, and above 20% is rated as insufficient quality and requires manual review of the electrode contact status collected. The rating is not based solely on the ρ value. If the longest duration of a single record within a suspected silent segment exceeds 10 seconds, even if ρ itself is below 5%, the rating will be directly downgraded to acceptable. This is because a single suspected detachment lasting more than 10 seconds is more likely to correspond to a genuine lead detachment rather than a brief contact jitter than detachments of the same total duration but scattered across multiple segments. For detachment records scattered across multiple segments, as long as none of the individual segments exceed 10 seconds, the rating will still be based on ρ when the total duration approaches the 5% limit. The signal quality rating provides a result for the entire record, while also recording the specific value of ρ and the specific record that triggered the downgrade, making the rating basis traceable. The calculation window for ρ is the cumulative duration of the current record from its start to the current moment. The rating is recalculated every minute as the record progresses. In the early stages, when the record duration is still short and the proportion of a single suspected silent segment is naturally higher, ρ is calculated using the larger of the recorded duration and the ten-minute duration as the denominator within the first ten minutes of the rating, while the boundary remains the complete three-level system. After the record duration reaches ten minutes, the calculation reverts to using the actual cumulative duration.

[0064] A graded EEG sequence was established to classify difficult-to-dissociate overlapping segments based on signal quality rating. The overall signal quality rating assigned to each segment was the initial grade, and subsequent downgrades were built up segment by segment above this initial grade. Downgrading was based on a frequency overlap distribution map and a list of high-risk notches. When the number of frequency bands covered by a high-risk notch interval in a segment exceeded 40% of the number of overlapping intervals in that segment, it indicated that the artifact was extensively entangled with the brainwaves in the frequency domain. This is commonly seen in subjects with sustained small movements rather than sudden jumps. In this case, no amount of differential downsampling can completely suppress the notch, and such segments were marked as difficult-to-dissociate overlapping segments. The graded EEG sequence used the waveform of the interpolated repair sequence, the repair rate of each recording, quality markers, and annotations of suspicious silent segments as its basis. The overall signal quality rating served as the baseline, and the same rating was used throughout. Downgrading was only applied to the time intervals corresponding to difficult-to-dissociate overlapping segments. Excellent segments were downgraded to acceptable, while acceptable and insufficient segments maintained their original grades without further downgrading. Downgrading only affects the local labeling of the segment itself, without changing the overall rating value given by the signal quality rating. Both are independently labeled and retained in the graded EEG sequence. When the time ranges of difficult-to-disseminate overlapping segments and suspicious silent segments overlap within a record, the overlapping part is labeled based on the suspicious silent segment, indicating a more severe signal loss than incomplete suppression; the non-overlapping parts retain their original labels. If there are no difficult-to-disseminate overlapping segments in the entire record, the graded EEG sequence maintains the original rating given by the signal quality rating throughout. The restoration rate of each record is written into the graded EEG sequence along with the waveform, and this is used as input for the next step of determining unsuppressed segments. The 40% proportion is determined by the distribution of high-risk segments confirmed to be completely suppressed by the notch filter scheme and segments confirmed to be incompletely suppressed; the two types of samples can be distinguished at around 40%.

[0065] The graded EEG sequence is processed by explicitly notching unsuppressed segments and outputting a suppressed EEG sequence. Unsuppressed segments are determined by the repair rate of the graded EEG sequence. If a record's repair rate is less than half the median repair rate of all records in the sequence, and it also contains a label indicating a suspicious silent segment, this segment is classified as unsuppressed and is no longer considered normal, clean data. Records with low repair rates but no labeling of suspicious silent segments are not considered unsuppressed segments; these segments are more likely due to irregular rhythms causing repair difficulties and remain in the graded EEG sequence with their original downgraded labeling. Notch labeling extracts the time interval corresponding to the unsuppressed segment from the graded EEG sequence, without filling in any interpolated or repaired values, and rewrites it as an explicit notch marker. The marker includes the start and end times and the specific criteria for classifying it as unsuppressed. The suppressed EEG sequence retains the data from non-unsuppressed segments in the graded EEG sequence as is, replaces the unsuppressed segments with notch markers, and outputs them sequentially. When encountering gap markers in suppressed EEG sequences, they should be treated directly as gap intervals. Treating them as zero values ​​will lead to misinterpretations of deep suppression downstream; treating them as interpolated waveforms will result in misinterpretations of normal rhythm. Both misinterpretations are more dangerous than explicit gaps. Gap markers only contain start and end times and the criteria for judgment, without any amplitude readings. If the signal quality rating indicates insufficient quality, or if the duration of degraded markers in the graded EEG sequence exceeds 30% of the entire recording, an additional overall prompt will be added to the output of the suppressed EEG sequence, indicating the rating result and the corresponding ρ value, suggesting that this recording is generally unsuitable for downstream analysis with high requirements for signal continuity.

[0066] To execute the methods described in the above embodiments and achieve the corresponding functions and technical effects, this application also provides the following system. See [link to system]. Figure 3 , Figure 3 This paper illustrates a structural block diagram of a single-channel EEG signal real-time artifact suppression system 300 provided in an embodiment of this application, comprising: The calibration establishment unit 301 is used to perform frame-segmentation processing on the single-channel real-time EEG acquisition sequence to obtain an EEG frame sequence, and to use the EEG frame sequence to extract the power frequency stable segment to establish a real-time calibration baseline. The feature collection unit 302 is used to perform frame-by-frame baseline drift amplitude quantization based on the real-time calibration baseline to form a drift amplitude sequence, locate the drift abrupt change site based on the drift amplitude sequence to generate an EMG artifact candidate set, and mark the template morphology highly consistent segment for the EMG artifact candidate set to obtain an artifact morphology feature set. The filtering protection unit 303 is used to implement differentiated reduction of the output notch filter parameters based on the artifact morphology feature set in the high overlap frequency band, analyze the self-ringing tail span through the notch filter parameters to obtain the filtered contaminated section, and delineate the coexistence boundary between residual artifacts and brain waves based on the filtered contaminated section to obtain the residual protection zone. The interpolation repair unit 304 is used to delineate the boundaries of the safe repair segments according to the degree of separability of the residual protection partition to construct a repair candidate sequence, perform rhythmic phase anchoring completion on the repair candidate sequence to form an interpolation repair sequence, and use the interpolation repair sequence to determine the suspicious silent segment based on the segment with the lowest repair rate. The sequence output unit 305 is used to analyze the risk of detachment and coverage of the suspected silent segment, output signal quality rating, perform downgrading labeling of difficult-to-separate overlapping segments to establish a graded EEG sequence, and perform explicit gap labeling of unsuppressed segments to output a suppressed EEG sequence based on the graded EEG sequence.

[0067] The single-channel EEG signal real-time artifact suppression system 300 described above can implement a single-channel EEG signal real-time artifact suppression method according to the above method embodiments. The options in the above method embodiments are also applicable to this embodiment, and will not be detailed here. The remaining contents of this application embodiment can be referred to the contents of the above method embodiments, and will not be repeated in this embodiment.

[0068] The purpose of the above embodiments is to reproduce and derive the technical solution of the present invention by way of example, and to fully describe the technical solution, purpose and effect of the present invention. The purpose is to enable the public to have a more thorough and comprehensive understanding of the disclosure of the present invention, and not to limit the scope of protection of the present invention.

Claims

1. A method for real-time artifact suppression of single-channel EEG signals, characterized in that, include: The single-channel real-time EEG acquisition sequence is processed into frames to obtain an EEG frame sequence. The power frequency stable segment is extracted from the EEG frame sequence to establish a real-time calibration baseline. Based on the real-time calibration baseline, frame-by-frame baseline drift amplitude quantization is performed to form a drift amplitude sequence. The drift abrupt change site is located according to the drift amplitude sequence to generate an EMG artifact candidate set. The template morphology highly consistent segments are labeled for the EMG artifact candidate set to obtain an artifact morphology feature set. Based on the artifact morphology feature set, the output notch filter parameters are reduced by depth differentiation in the high-overlap frequency band. The notch filter parameters are used to analyze the self-ringing tail span to obtain the filtered contamination section. Based on the filtered contamination section, the boundary between the coexistence of residual artifacts and brain waves is delineated to obtain the residual protection zone. The remaining protection partition is delineated according to the degree of separability to construct a repair candidate sequence by defining the boundary of the safe repair segment. Rhythmic phase anchoring completion is performed on the repair candidate sequence to form an interpolation repair sequence. The interpolation repair sequence is used to determine the suspicious silent segment based on the segment with the lowest repair rate. Based on the analysis of the suspected silent segments, the risk of detachment and coverage is analyzed to output signal quality rating. For the signal quality rating, the difficult-to-separate overlapping segments are downgraded and a graded EEG sequence is established. Based on the graded EEG sequence, the unsuppressed segments are explicitly marked with gaps to output a suppressed EEG sequence.

2. The method according to claim 1, characterized in that, The step of using the EEG frame sequence to extract the power frequency stable segment and establish a real-time calibration baseline includes: The intra-frame power frequency peak-valley symmetry deviation of the EEG frame sequence is verified to obtain a power frequency cycle candidate set; Based on the power frequency cycle candidate set, cross-frame phase drift classification is performed to form a power frequency stability classification; The lowest frame segment of phase jump variable is extracted using the power frequency stability classification to generate a power frequency stability reference set; After removing outlier frames from the power frequency stable reference set, a real-time calibration baseline is established by sparse integration.

3. The method according to claim 1, characterized in that, The step of performing frame-by-frame baseline drift amplitude quantization to form a drift amplitude sequence based on the real-time calibration baseline includes: Based on the real-time calibration baseline, a frame-by-frame power frequency baseline stripping is performed to form a power frequency baseline-free frame set; Based on the power frequency baseline frame set, the drift direction is recorded by identifying the drift acceleration point with shortened reversal spacing. The drift direction record is used to identify local anti-rise points within the low-amplitude drift segment and generate an anti-rise drift annotation set. The inverse drift annotation set is subjected to frame-by-frame amplitude quantization to form a drift amplitude sequence.

4. The method according to claim 1, characterized in that, The process of annotating the template morphology of the candidate electromyographic artifacts to obtain the artifact morphology feature set includes: For the aforementioned electromyography artifact candidate set, candidate frame template matching and quantization are performed to generate a morphological deviation comparison set; The morphological deviation ratio pair is used to identify frame segments with persistently low deviation values ​​to obtain high-consistency candidate segments; The highly consistent candidate segments are used to filter out rhythmic morphology self-consistent segments to form a pure pseudo-trace highly consistent set; Based on the high consistency set of pure artifacts, the artifact morphology feature set is obtained by categorizing them according to the time span of their occurrence.

5. The method according to claim 1, characterized in that, The step of implementing differentiated reduction of output notch filter parameters based on the artifact morphology feature set in the high overlap frequency band includes: Based on the artifact morphology feature set, an overlap distribution map of frequency bands is established by identifying the overlap intervals between each artifact component and brainwave frequency bands. The frequency band overlap distribution map is used to suppress the upper limit of the notch depth in the partition to form a restricted notch configuration. For the restricted notch configuration, a high-overlap frequency band notch center offset locking generation frequency-locked notch scheme is implemented; The frequency-locked notch filter scheme is verified and the bandwidth boundary is analyzed to output notch filter parameters by frequency band-by-frequency depth reduction.

6. The method according to claim 1, characterized in that, The step of obtaining the filtered contamination section by analyzing the self-ringing tail span through the notch filter parameters includes: The impulse response duration set is obtained by calculating the time it takes for each order of impulse response to decay to the baseline using the notch filter parameters. Based on the aforementioned impulse response duration set, notch waves with attenuation durations exceeding the theoretical duration are identified to form ringing over-limit groups; For the ringing exceeding the limit group, the excessive ringing extension span is measured step by step to generate a trailing extension diagram; The trailing extension map is labeled with ringing pollution segments outside the boundary to obtain filtered pollution segments.

7. The method according to claim 1, characterized in that, The step of performing rhythmic phase anchoring completion to form an interpolated repair sequence for the repair candidate sequence includes: The rhythm anchor point set is obtained by identifying the complete rhythm cycle boundaries before and after each segment of the repair candidate sequence; Based on the set of rhythmic anchor points, a phase transition scheme is generated by planning a phase continuity constraint path across the defect segment. The phase transition scheme is used to perform frame-by-frame defect filling of the rhythmic phase reference to obtain a phase-filled frame set; The phase-filled frame set is then labeled with over-amplitude filler frames and embedded to form an interpolation repair sequence.

8. The method according to claim 5, characterized in that, The method of using the frequency band overlap distribution map to suppress the upper limit of the notch depth in a partitioned area to form a restricted notch configuration includes: Based on the frequency band overlap distribution map, the real-time proportion of brain wave components in each overlap interval is determined to form the brain wave proportion distribution. The high-risk notch intervals with brainwave proportions exceeding the threshold are determined by the brainwave proportion distribution, and a high-risk notch list is established. Based on the high-risk notch list, perform depth difference analysis for each interval to obtain the partition depth constraint set; The partition depth constraint set is integrated by gradient widening to form a restricted notch configuration.

9. The method according to claim 7, characterized in that, The step of generating a phase transition scheme by planning a phase continuity constraint path across the defect segment based on the rhythm anchor point set includes: The periodic duration deviation set is obtained by verifying the periodic deviation between adjacent anchor points in the rhythmic anchor point set. Based on the aforementioned period duration deviation set, abnormal anchor point segments with abrupt changes in period duration are identified to obtain the abnormal rhythm boundary set; The stability ratio of each candidate path is enumerated and calculated using the abnormal rhythm boundary set to establish a stable path hierarchical record; A phase transition scheme is generated by performing a path selection process with the minimum phase difference fluctuation on the stable path classification record.

10. A single-channel EEG signal real-time artifact suppression system, characterized in that, include: The calibration establishment unit is used to perform frame-segmentation processing on the single-channel real-time EEG acquisition sequence to obtain an EEG frame sequence, and to use the EEG frame sequence to extract the power frequency stable segment to establish a real-time calibration baseline. The feature collection unit is used to perform frame-by-frame baseline drift amplitude quantization based on the real-time calibration baseline to form a drift amplitude sequence, locate the drift abrupt change site based on the drift amplitude sequence to generate an EMG artifact candidate set, and mark the template morphology highly consistent segment for the EMG artifact candidate set to obtain an artifact morphology feature set. The filtering protection unit is used to implement differentiated reduction of output notch filter parameters in high-overlap frequency bands based on the artifact morphology feature set, analyze the self-ringing tail span through the notch filter parameters to obtain the filtered contamination section, and delineate the coexistence boundary between residual artifacts and brain waves based on the filtered contamination section to obtain the residual protection zone. The interpolation repair unit is used to delineate the boundaries of the safe repair segments according to the degree of separability of the residual protection partition to construct a repair candidate sequence, perform rhythmic phase anchoring completion on the repair candidate sequence to form an interpolation repair sequence, and use the interpolation repair sequence to determine the suspicious silent segment based on the segment with the lowest repair rate. The sequence output unit is used to analyze the risk of detachment and coverage based on the suspected silent segment, output signal quality rating, perform downgrading labeling of difficult-to-separate overlapping segments to establish a graded EEG sequence, and perform explicit gap labeling of unsuppressed segments to output a suppressed EEG sequence based on the graded EEG sequence.

Images