Disclosure of Invention
Aiming at the defects of the prior art, the invention provides an automatic sleep stage method, an automatic sleep stage system and a storage medium for rodents, which can carry out rapid and accurate automatic sleep stage on long-time sleep monitoring data of rodents, can effectively reduce the time consumption of manual stage and improve the research efficiency of the related biological fields such as sleep and the like.
In order to achieve the above purpose, the technical scheme adopted by the invention is as follows:
An automated sleep staging method for a rodent comprising the steps of:
Preprocessing and feature extraction are carried out on the acquired sleep monitoring original data to obtain EEG feature values and EMG level feature values;
Determining an EMG level threshold and an EEG characteristic threshold by adopting an adaptive threshold calculation method, and performing sleep classification and classification rough classification according to EMG data, wherein the sleep classification and rough classification are divided into a awake period and a sleep period;
according to the rough classification result, threshold values are automatically set for all the characteristic values by combining the data distribution of the characteristic values, and the clear stage of sleep is realized by combining the threshold values.
In the above technical solution, further, the method includes collecting long-term sleep monitoring raw data of the rodent, wherein the data includes an electromyography EMG channel data and at least one electroencephalogram EEG channel data;
In preprocessing, when a plurality of EEG channels exist, one EEG channel with fewer artifacts is selected, EMG channel data and the selected EEG channel data are filtered and downsampled, the filtering is that the selected EEG channel data is subjected to low-pass filtering to remove high-frequency interference, the cut-off frequency is set to be 40Hz or lower, high-pass filtering is performed to remove extremely-low-frequency interference, the cut-off frequency is set to be 0.4Hz or lower, and the EMG channel data is subjected to 10-100Hz band-pass filtering.
Further, the feature extraction includes the following steps:
after FFT conversion of EEG channel data in a selected time window, delta band energy Pow_delta is calculated as an EEG characteristic value;
calculating the root mean square value or the average absolute value of the EMG channel data waveform in the selected time window to obtain an EMG level EMG_level which is an EMG level characteristic value;
and performing smoothing processing after extracting the characteristic values of all frames.
Further, the adaptive threshold calculation method includes the following steps:
(1) Determining an EMG level change baseline based on EMG level data, judging the level distinction of the EMG baseline, acquiring an EMG baseline level threshold value, detecting a high EMG section and calculating an EMG baseline dynamic threshold value;
(2) Performing phase two classification on sleep data by using the high EMG section and the EMG baseline dynamic threshold value, and roughly classifying into a waking period and a sleep period;
(3) And according to the rough classification result, combining the data distribution of the characteristic values to automatically determine an EMG level threshold and an EEG characteristic threshold.
Further, the specific steps in the step (1) are as follows:
s401, coarsening EMG level data by adopting a method of solving a median value by adopting a long-time moving window to obtain an EMG level change base line L;
s402, arranging the values of the base line L according to the size, taking the smaller and larger ranking values to respectively represent the low and high EMG levels, calculating the difference between the two values, and prompting if the difference is smaller than a preset discrimination decision threshold;
S403, acquiring a first-order differential signal D of a base line L and an absolute value D_abs thereof, and further acquiring a larger threshold value TH_D_hi and a smaller threshold value TH_D_lo according to the data distribution of the D_abs;
S404, detecting the starting and ending positions of the high EMG section based on a double-threshold endpoint detection method by combining the differential signal D obtained in the step S403 and two thresholds TH_D_hi and TH_D_lo;
s405, setting a vector LW to represent the awake level, wherein the length is the same as L, wherein the element value is 0 or 1,0 represents the sleep period, 1 represents the awake period, initializing all elements to 0, then initializing the element value corresponding to the index position of the high EMG section to 1 to obtain a coarse-granularity binary sleep stage sequence, initializing a coarse-granularity dynamic EMG baseline threshold vector BL to L, and correcting the BL value of each high EMG section to be the average value of L vector elements corresponding to the nearest 5 moments before the start of the section.
Further, the (2) specifically includes the following steps:
S406, sleeping time phase rough division
1) Restoring the coarse granularity vectors BL and LW acquired in the step S405 to fine granularity vectors corresponding to the frame length by using a linear interpolation mode, wherein the fine granularity vectors are BLrs and LWrs respectively, and the value higher than 0.5 in LWrs is set to be 1 and the value lower than 0.5 is set to be 0;
2) Calculating sleep frame index, and if the frame index i is satisfied
LWRs (i) =0 and EMG_level (i) -BLrs (i) < TH_D_lo
Wherein LWRs (i) is a fine granularity awake level value at the i-th moment, EMG_level (i) is a myoelectricity level value at the i-th moment, BLrs (i) is a fine granularity EMG baseline threshold at the i-th moment, and the frame is set as a sleep frame and frames except the sleep frame are set as awake frames;
3) And respectively acquiring EMG level characteristic values corresponding to sleep and awake frame indexes, and correspondingly defining the EMG level characteristic values as two vectors BLs and BLw.
Further, the specific steps in (3) are as follows:
s407, calculating an EMG level threshold, namely calculating a higher threshold TH_EMG_hi for distinguishing sleep from wakefulness according to the BLs data distribution characteristics, and a lower threshold TH_EMG_lo for distinguishing REM from NREM;
S408, adaptively obtaining EEG feature thresholds for distinguishing REM and NREM according to the statistical distribution features of the data:
1) Acquiring a value of the sleep frame index position calculated in the corresponding S406 of the EEG feature vector pow_delta, and setting the value as pow_delta_s;
2) Removing extreme points in Pow_delta_s, dividing the residual data into 10 equally divided intervals, counting the number of data points falling into each interval, assuming that the interval vector is x, the number of corresponding intervals is vector n, and calculating n-means M;
3) Normalizing the vector n by M to obtain a vector r;
4) Finding a first index position which satisfies the quantity higher than TH in the n vectors, and setting the first index position as p, wherein the product of the maximum value in the TH vector r subtracted by 0.5 and the product of the maximum value and M is obtained by downward rounding;
5) Setting a cyclic variable i, traversing p to 10, and calculating a high eigenvalue region, namely an eigenvalue F1 with larger probability in the NREM period, wherein the high eigenvalue region is the ratio of the accumulated product value of a vector x and a vector n to the sum value of the vector n from p to 10;
6) Obtaining a low eigenvalue area, namely an eigenvalue F2 with larger REM probability, wherein the eigenvalue F2 is the ratio of the sum of the accumulated product value of the vector x and the vector n to the sum of the vector n from 1 to p-1 of the variable i;
7) Taking the average value of F1 and F2 as a candidate threshold value TH1, taking the average value of the pow_delta_s vector plus one time of the standard deviation as a candidate threshold value TH2, and taking the smaller one of the TH1 and TH2 as an EEG characteristic threshold value TH_delta.
Further, detecting and correcting the region where the REM period is likely to be based on the extracted feature values of all frames and the determined EMG level threshold, including the following steps:
1) Setting EEG characteristic value component vectors corresponding to index positions with EMG characteristic values lower than a higher myoelectricity level threshold value TH_EMG_hi as Pow_delta_s1;
2) Performing peak detection on the reverse waveform of the Pow_delta_s1 waveform, and obtaining a peak position according to a set condition;
3) And searching the starting position and the ending position of the low delta energy region leftwards and rightwards for each detected peak point position respectively, judging whether the starting and ending position is required to be corrected according to the set conditions, and if so, correcting according to the correction rules until the conditions are met without correcting, so that the low delta energy region position index meeting the conditions is finally obtained.
Furthermore, sleep is specifically staged by combining an EMG level threshold value and an EEG characteristic threshold value, and the method comprises the following steps:
First, all frames are set to NREM period, those frames in which EMG level eigenvalue is higher than TH_EMG_hi are modified to awake period, then REM is set, and then, for each obtained start-stop index position of low delta energy region, the method further filters the continuous frames in the corresponding start-stop index position by removing frames in which wakefulness has been classified, removing frames in which delta energy is higher than TH_delta, removing frames in which myoelectricity level is higher than the average of both TH_EMG_lo and TH_EMG_hi, and classifying the rest frames to REM.
The invention also provides a rodent automated sleep staging system for use in implementing a rodent automated sleep staging method as defined in any one of the preceding claims.
The present invention also provides a computer readable storage medium having stored thereon a computer program, characterized in that the program when executed by a processor implements the rodent automated sleep staging method as defined in any one of the above.
The invention also provides an electronic device comprising:
one or more processors;
A memory for storing one or more programs;
the one or more programs, when executed by the one or more processors, cause the one or more processors to implement the rodent automated sleep staging method of any one of the preceding claims.
The invention has the following beneficial effects:
The invention provides a rodent automatic sleep stage method based on characteristic engineering and self-adaptive threshold setting. Compared with a method based on model learning, the method does not need to obtain a large number of marked training data sets in construction, has clear physical significance of the features and has strong interpretation of the stage judgment result. Meanwhile, the self-adaptive threshold setting mode provided by the invention can carry out sleep stage classification rough classification (namely waking and sleeping) according to EMG data acquired by each file, further automatically set reasonable thresholds for each characteristic value according to rough classification results, and be used for subsequent final sleep stage judgment. The method can effectively get rid of and avoid the trouble of setting the manual threshold value due to the need of considering individual and acquisition differences, and can realize stage rapidly and accurately.
Detailed Description
The technical scheme of the invention is further described in detail below with reference to the accompanying drawings and specific embodiments.
The invention discloses an automatic sleep stage method for rodents, which mainly comprises the following steps of data reading, data preprocessing, feature extraction, self-adaptive threshold setting and sleep stage judgment (shown in figure 1).
1. Data reading (S10)
Firstly, the acquired rodent long-term sleep monitoring original data are read into a computer memory, and the data format is not limited as long as the data protocol is disclosed. The raw data needs to contain at least one EEG channel and one EMG channel.
2. Data preprocessing (S20)
The preprocessing steps mainly include EEG channel selection, filtering and downsampling (as in fig. 2).
S201 channel selection
If multiple EEG channels exist in the original data, some simple rules can be set to select one channel with fewer artifacts for subsequent analysis. For example, a larger threshold value exceeding the conventional electroencephalogram amplitude range is set, and one of the channels whose total length of time is the shortest in all the sections exceeding the threshold value is selected.
S202 filtering
In view of the small upper limit of the cut-off frequency of the frequency band involved in the feature extraction of sleep electroencephalogram (EEG) stage, the cut-off frequency can be set at 40Hz or lower, the cut-off frequency can be set at 0.4Hz or lower, and the EMG channel can be subjected to bandpass filtering, such as 10-100Hz.
S203 downsampling
The sampling rate of the original collected data is generally higher and reaches 1kHz or even higher, the concerned characteristic frequency band in the sleep stage does not need to be characterized by the high sampling rate, and in order to reduce the calculation amount, save the memory occupation of the system and reduce the computer configuration requirement, the original EEG and EMG data can be downsampled to a certain degree, such as to 200Hz or even 100Hz.
3. Feature extraction (S30)
The feature extraction mainly comprises the steps of intercepting data, carrying out FFT conversion on EEG channel data in a selected time window, calculating delta frequency band, namely 0.5-4 Hz energy pow_delta, calculating EMG channel data in the selected time window to obtain EMG level EMG_level, and carrying out smoothing treatment on the extracted feature values of all frames (as shown in figure 3).
S301 data interception
The decision on sleep stage requires a small data area (analysis time window) to be selected first, the smaller the window length, the finer the sleep stage result, and for rodent sleep stage, more 10s or 4s are used, and generally non-overlapping windows are adopted.
S302:FFT
And converting the EEG channel waveform in the selected time window into a frequency domain through FFT, and calculating the energy corresponding to each frequency point.
S303 band energy calculation
And calculating delta frequency band (0.5-4 Hz) energy Pow_delta according to the frequency and energy corresponding to each frequency point after FFT, namely the energy sum of the frequency points in the frequency range.
S304 EMG level calculation
The EMG level emg_level may be obtained by calculating the root mean square value or average absolute value of the EMG waveform over a selected time window.
S305, eigenvalue smoothing
To reduce outlier interference and prevent too frequent switching of sleep stages, after extracting eigenvalues of all frames, the eigenvalues pow_delta and emg_level are moderately smoothed, which may be performed in a weighted moving average manner or in a convolution manner using a kernel function such as a gaussian kernel.
4. Adaptive threshold calculation method (S40)
The self-adaptive threshold setting provided by the invention firstly carries out two classifications on sleep data based on EMG level data, roughly classifies the sleep data into awake and sleep periods, and then determines thresholds of EMG level and brain electrical characteristics. The method mainly comprises the steps of calculating an EMG level change base line, calculating an EMG base line level distinguishing degree, calculating an EMG base line level threshold value, detecting an EMG level section, calculating an EMG base line dynamic threshold value, rough sleep phase division, calculating an EMG level threshold value and determining an EEG characteristic threshold value (shown in figure 4).
S401 calculating EMG level Change Baseline
Coarse granularity is carried out on the EMG level change trend by adopting a method of calculating the median value by using a moving window, the window length is set to 300s, and the window moving step length is set to 100s, so that an EMG level change baseline L is obtained.
S402, calculating EMG baseline level differentiation
In view of the fact that the adaptive threshold setting mode used in the invention is relatively dependent on the change of the EMG level, if the change of the EMG level is not distinguished obviously due to the data acquisition problem, the subsequent threshold setting is affected, and then the sleep stage result is affected. Therefore, the degree of discrimination of the EMG baseline level in the file is detected here by prompting if it is too small.
1) Arranging the values of the EMG level change base line L from small to large;
2) Values were obtained for the 10% and 90% ranks, representing low and high EMG levels, respectively;
3) Calculating the difference between the high and low EMG levels, namely the level differentiation of the EMG base line;
4) If the threshold value is lower than the threshold value, prompting is carried out, and the automatic stage result obtained under the condition is unreliable and a larger range of result correction is needed.
S403, calculating an EMG baseline level difference threshold
The two thresholds are calculated for detecting the high myoelectricity level section by adopting a double-threshold method, and the method is as follows
1) Acquiring a front-rear point differential signal D of L;
2) Obtaining an absolute value D_abs of the differential signal;
3) Calculating a smaller threshold value TH_D_lo which is the median value of D_abs;
4) Calculating a larger threshold value TH_D_hi, namely, firstly, carrying out ascending arrangement on D_abs, and then obtaining a 98% ranking value;
5) Correcting the TH_D_hi by firstly calculating 50% value of the TH_D_hi and average value of the TH_D_lo and the TH_D_hi, and taking larger value as corrected TH_D_hi compared with the average value of the TH_D_lo and the TH_D_hi;
6) The TH_D_lo is corrected by firstly calculating a3 times value of the TH_D_lo, and then comparing the value with the corrected TH_D_hi, and taking a smaller value as the corrected TH_D_lo, namely TH_D_lo is less than or equal to TH_D_hi.
S404, detecting a high EMG horizontal segment
The main calculation process is as follows:
1) Acquiring each moment of D > TH_D_hi as a candidate starting moment of detection of each high EMG section;
2) Whether each candidate time can be used as the starting time of the high EMG section detection is judged by the following method, assuming that the index position corresponding to the current candidate time is id, if the following condition is satisfied,
(L (id+1) -L (id) > TH_D_lo and L (id+2) -L (id+1) > TH_D_lo) or L (id+1) -L (id) > TH_D_hi
The starting time of a high EMG section is judged to be ST, the detection of the subsequent ending time is carried out, and otherwise, the judgment of the next candidate starting time is carried out;
3) The termination time of each high EMG section is calculated by starting from ST+1 index to L length, assuming the loop index is j, and determining whether the following condition is satisfied
L (j) -L (j-1) < = -th_d_lo and L (j+1) -L (j) > -th_d_lo and abs (L (j) -L (ST)) < th_d_hi
If the condition is met, finding a termination index SE of the high EMG section corresponding to the ST start, terminating the cycle of the termination index detection, returning to 2), and judging the next candidate start time;
s405 calculating EMG baseline dynamic threshold
The method mainly comprises the following steps:
1) Assuming that the vector LW represents the awake level, the length is equal to the length of L, the element values in the vector are 0 or 1,0 represents the sleep period, 1 represents the awake period, all elements are initialized to 0 first, and then the index positions of all the high myoelectricity level intervals detected in S404 are 1, so as to obtain a coarse-granularity binary sleep stage sequence.
2) Assume that another vector BL characterizes the coarse-grained dynamic EMG baseline threshold, initialized to the value of the L vector. The threshold value of the high EMG level interval obtained in S404 is then calculated by, for each high EMG level interval, assuming ST as its start index, SE as its end index, traversing the loop variable j from ST-1 to 5, if the following conditions are satisfied
LW (j) +LW (j-1) +LW (j-2) +LW (j-3) +LW (j-4) =0 and
Abs ((L (j) -L (j-1))) < TH_D_lo and
Abs ((L (j-1) -L (j-2))) < TH_D_lo and
Abs ((L (j-2) -L (j-3))) < TH_D_lo and
abs((L(j-3)-L(j-4)))<TH_D_lo
The method comprises the steps of taking the average value of the current index j and 5L vector elements corresponding to the four index positions in front of the current index j (taking the average value as BL values of all moments of the current high EMG zone (ST to SE moment)), and then entering BL value correction of the next high EMG zone.
S406, sleeping time phase rough division
1) The coarse-granularity vectors BL and LW acquired in S405 are restored to fine-granularity vectors corresponding to the frame length, assuming BLrs and LWrs, respectively, by using a linear interpolation method. A value of LWrs above 0.5 is set to 1 and a value below 0.5 is set to 0.
2) Calculating sleep frame index, assuming frame index i, if it is satisfied
LWRs (i) =0 and EMG_level (i) -BLrs (i) < TH_D_lo
Wherein LWRs (i) is a fine granularity awake level value at the i-th moment, EMG_level (i) is a myoelectricity level value at the i-th moment, BLrs (i) is a fine granularity EMG baseline threshold at the i-th moment, and the frame is set as a sleep frame and frames except the sleep frame are set as awake frames;
3) And respectively acquiring EMG level characteristic values corresponding to sleep and awake frame indexes, and correspondingly defining the EMG level characteristic values as vectors BLs and BLw.
S407 calculating EMG level threshold
To adaptively obtain two EMG level thresholds, a higher threshold for distinguishing sleep from wakefulness and a lower threshold for distinguishing REM from NREM, the following method is set:
1) Calculating an EMG level candidate threshold TH1 by taking the average value of the median value of BLs and the median value of BLw;
2) Calculating another EMG level candidate threshold TH2 by adding 3 times of standard deviation to the BLs mean;
3) Two thresholds of EMG levels are obtained for later sleep stage decisions. The higher threshold of EMG level TH_EMG_hi is taken as the smaller of TH1 and TH2 for distinguishing sleep from wakefulness, and the lower threshold of EMG level TH_EMG_lo is taken as the median of BLs for distinguishing REM from NREM.
S408 determining EEG characteristic threshold
Based on experience, the number of frames in NREM period is far greater than that in REM period, and delta energy in NREM period is higher than that in REM, the following method is set for adaptively obtaining threshold value for distinguishing EEG characteristics of REM and NREM
1) Acquiring a value of the sleep frame index position calculated in the corresponding S406 of the EEG feature vector, and setting the value as Pow_delta_s;
2) The data distribution of the Pow_delta_s is calculated, the extreme points in the Pow_delta_s are removed firstly, the rest data are divided into 10 equally divided intervals, the number of data points falling into each interval is counted, the interval vector is assumed to be x, and the number of corresponding intervals is assumed to be vector n;
3) Calculating the ratio of the number of each interval to the n average value, and setting the ratio as a vector r, wherein r represents the average value number of which the number of each interval is several times;
4) Setting a quantity threshold value TH, finding a first index position meeting the quantity higher than TH in an n vector, and setting p, wherein a calculation formula of TH is as follows
TH = mean(n)*(floor(max(r))-0.5)
5) Setting a loop variable i, traversing p to 10, obtaining the following two accumulated values S1 and C1, the formula is as follows
The eigenvalue F1 of the high eigenvalue area (with higher probability of NREM period) is calculated as follows
6) In a similar way, a feature value F2 of a low feature value region (with a large REM probability) is obtained as follows
7) The delta energy threshold, TH _ delta, is calculated, distinguishing NREM from REM, as follows,
Calculating a candidate threshold TH1, which is the average value of F1 and F2;
Calculating a candidate threshold value TH2, and adding one time of standard deviation to the average value of the pow_delta_s vector;
Th_delta is set to be the smaller of both TH1 and TH 2.
5. Sleep stage determination (S50)
1) Acquiring an index position of which the myoelectricity level characteristic value is lower than the higher myoelectricity level threshold value TH_EMG_hi calculated in the step S407, acquiring a vector consisting of delta energy characteristic values of corresponding positions, and setting the vector as Pow_delta_s1;
2) Performing peak detection on the reverse waveform of the Pow_delta_s1 waveform, setting a limit condition, wherein the detected peak value must be highlighted around the peak by at least P, the value of P is 0.3 times of the difference between the maximum value and the minimum value in the Pow_delta_s vector, the distance between the peaks needs to be at least 30 seconds, and the positions of the detected peak values can be the region where REM phase exists;
detection of low delta energy regions
A. for each detected peak position, let be loc, find the start and end positions of the low delta energy region to the left and right, respectively.
For the starting position search, a loop variable j is set, the index of loc-1 to 2 is traversed, if the following condition is satisfied,
Pow_delta (j) < Pow_delta (j+1) and Pow_delta (j-1) < Pow_delta (j)
The traversal is terminated, and the current index j is recorded as the starting position index ST;
for the end position seek, a loop variable j is set, the index of loc+1 to frame length-1 is traversed, if the following condition is satisfied,
Pow_delta (j) < Pow_delta (j-1) and Pow_delta (j+1) < Pow_delta (j)
The traversal is terminated, and the current index j-1 is recorded as a termination position index SE;
b. the start-stop position of the above-obtained section is corrected as follows
Calculating the absolute difference of the characteristic value pow_delta of the start and stop index position of the section, and if the absolute difference is higher than 0.8 times of the P value calculated in the step 2, correcting the start and stop index position by the following method
If the Pow_delta value corresponding to the start index is lower than the Pow_delta value of the end index, continuing to search the start position forwards, and traversing the indexes of ST-2 to 2 by assuming a loop variable j similarly to the above method, if the following conditions are satisfied
Pow_delta (j) < Pow_delta (j+1) and Pow_delta (j-1) < Pow_delta (j)
The traversal is terminated, and the current index j is recorded as the corrected starting position index ST;
if the Pow_delta value corresponding to the start index is higher than the Pow_delta value of the end index, continuing to search the start position backwards, and traversing the indexes from SE+2 to frame number-1 by assuming a loop variable j similar to the above method, if the following conditions are satisfied
Pow_delta (j) < Pow_delta (j-1) and Pow_delta (j+1) < Pow_delta (j)
The traversal is terminated, and the current index j-1 is recorded as a corrected termination position index SE;
And (3) circulating the process until the absolute difference value of the Pow_delta corresponding to the updated start and stop index position is smaller than P0.8, stopping correction, and entering detection and correction of the next REM area.
4) Sleep staging, the method is as follows
Firstly, all frames are set as NREM period, those frames in which EMG level characteristic value is higher than TH_EMG_hi are corrected as awake period, then REM is set, the method is as follows, for each low delta energy region start-stop index position obtained in step 3, the continuous frames in the corresponding start-stop index position are further screened, the frames in which wakefulness has been classified are removed, the frames in which delta energy is higher than TH_delta obtained by calculation in S408 are removed, the frames in which myoelectricity level is higher than the average value of both TH_EMG_lo and TH_EMG_hi are removed, and the rest frames are classified into REM. Thus, the sleep stage of the whole file is completed.
The sleep stage algorithm result provided by the invention is largely calculated based on the self-adaptive threshold, so that the sleep data cannot be too short and each sleep phase should be covered in order to make the stage result as accurate as possible.
Examples
The method of the invention is used for carrying out automatic sleep stage on the sleep monitoring data of a 4-hour mouse, the frame length is set to be 10s, and the results are shown in figure 5, namely an EEG waveform diagram, an EEG time-frequency diagram, an EMG level diagram and a sleep time-phase sequence diagram (W: awake, N: NREM, R: REM) from top to bottom. Therefore, the method can obtain more accurate sleep stage results. One example of frame data (fig. 6) automatically divided into NREM is that of the current frame, the top of the graph is EEG waveform of the current frame, the middle of the graph is displayed with the time-frequency chart of the EEG waveform of the current frame and 20 frames on the left and right of the current frame, the left-most side of the lower graph is displayed with the spectrogram of the current frame, the delta frequency band with higher energy can be seen, the middle graph below is displayed with the relative energy of 3 frequency bands (delta: 0.5-4 Hz, theta: 5-9 Hz, alpha: 9-12 Hz), the right-hand graph below is displayed with myoelectricity level of the current frame and the front and back 5 frames thereof, the current frame is marked with red, the two transverse dotted lines in the graph are respectively upper and lower EMG level thresholds calculated in an adaptive manner, namely, the rose line is TH_EMG_hi, and the green line is TH_EMG_lo. One example of frame data automatically divided into REM (FIG. 7) shows that the brain wave theta rhythm is obvious, and the threshold setting of myoelectric level can be seen reasonably in the myoelectric level diagram. One example of frame data (fig. 8) that is automatically divided into wakefulness is that the brain wave theta rhythm is obvious, and the central frequency of theta frequency band in wakefulness can be seen to be slightly higher than REM period in frequency spectrum, and the myoelectricity level diagram at the lower right shows that the myoelectricity level of the current frame is higher than TH_EMG_hi. The interval of fig. 5 is enlarged, and an automatic staging result (fig. 9) before and after the REM period is displayed, so that the sensitivity of the automatic identification result of the method can be seen to be higher.
It will be appreciated by those skilled in the art that embodiments of the present invention may be provided as a method, system, or computer program product. Accordingly, the present invention may take the form of an entirely hardware embodiment, an entirely software embodiment or an embodiment combining software and hardware aspects. Furthermore, the present invention may take the form of a computer program product embodied on one or more computer-usable storage media (including, but not limited to, disk storage, CD-ROM, optical storage, and the like) having computer-usable program code embodied therein.
The above embodiments are only some of the preferred embodiments of the present invention, but are not intended to limit the present invention. Various changes and modifications may be made by one of ordinary skill in the pertinent art without departing from the spirit and scope of the present invention. Therefore, all the technical schemes obtained by adopting the equivalent substitution or equivalent transformation are within the protection scope of the invention.