WO2017135116A1 - 生体信号処理方法および生体信号処理装置 - Google Patents

生体信号処理方法および生体信号処理装置 Download PDF

Info

Publication number
WO2017135116A1
WO2017135116A1 PCT/JP2017/002483 JP2017002483W WO2017135116A1 WO 2017135116 A1 WO2017135116 A1 WO 2017135116A1 JP 2017002483 W JP2017002483 W JP 2017002483W WO 2017135116 A1 WO2017135116 A1 WO 2017135116A1
Authority
WO
WIPO (PCT)
Prior art keywords
biological signal
data
time
value
averaging
Prior art date
Application number
PCT/JP2017/002483
Other languages
English (en)
French (fr)
Inventor
松浦　伸昭
小笠原　隆行
啓 桑原
Original Assignee
日本電信電話株式会社
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by 日本電信電話株式会社 filed Critical 日本電信電話株式会社
Priority to CN201780009657.6A priority Critical patent/CN108601546B/zh
Priority to JP2017565497A priority patent/JP6687645B2/ja
Priority to US16/075,116 priority patent/US10918302B2/en
Publication of WO2017135116A1 publication Critical patent/WO2017135116A1/ja

Links

Images

Classifications

    • A—HUMAN NECESSITIES
    • A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
    • A61B—DIAGNOSIS; SURGERY; IDENTIFICATION
    • A61B5/00—Measuring for diagnostic purposes; Identification of persons
    • A61B5/24—Detecting, measuring or recording bioelectric or biomagnetic signals of the body or parts thereof
    • A61B5/316—Modalities, i.e. specific diagnostic methods
    • A61B5/318—Heart-related electrical modalities, e.g. electrocardiography [ECG]
    • A61B5/346—Analysis of electrocardiograms
    • A61B5/349—Detecting specific parameters of the electrocardiograph cycle
    • A61B5/352—Detecting R peaks, e.g. for synchronising diagnostic apparatus; Estimating R-R interval
    • A—HUMAN NECESSITIES
    • A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
    • A61B—DIAGNOSIS; SURGERY; IDENTIFICATION
    • A61B5/00—Measuring for diagnostic purposes; Identification of persons
    • A61B5/24—Detecting, measuring or recording bioelectric or biomagnetic signals of the body or parts thereof
    • A61B5/316—Modalities, i.e. specific diagnostic methods
    • A61B5/318—Heart-related electrical modalities, e.g. electrocardiography [ECG]
    • A61B5/346—Analysis of electrocardiograms
    • A61B5/349—Detecting specific parameters of the electrocardiograph cycle
    • A—HUMAN NECESSITIES
    • A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
    • A61B—DIAGNOSIS; SURGERY; IDENTIFICATION
    • A61B5/00—Measuring for diagnostic purposes; Identification of persons
    • A61B5/45—For evaluating or diagnosing the musculoskeletal system or teeth
    • A61B5/4538—Evaluating a particular part of the muscoloskeletal system or a particular medical condition
    • A61B5/459—Evaluating the wrist
    • A—HUMAN NECESSITIES
    • A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
    • A61B—DIAGNOSIS; SURGERY; IDENTIFICATION
    • A61B5/00—Measuring for diagnostic purposes; Identification of persons
    • A61B5/72—Signal processing specially adapted for physiological signals or for diagnostic purposes
    • A61B5/7203—Signal processing specially adapted for physiological signals or for diagnostic purposes for noise prevention, reduction or removal
    • A—HUMAN NECESSITIES
    • A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
    • A61B—DIAGNOSIS; SURGERY; IDENTIFICATION
    • A61B5/00—Measuring for diagnostic purposes; Identification of persons
    • A61B5/72—Signal processing specially adapted for physiological signals or for diagnostic purposes
    • A61B5/7235—Details of waveform analysis
    • A—HUMAN NECESSITIES
    • A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
    • A61B—DIAGNOSIS; SURGERY; IDENTIFICATION
    • A61B5/00—Measuring for diagnostic purposes; Identification of persons
    • A61B5/72—Signal processing specially adapted for physiological signals or for diagnostic purposes
    • A61B5/7271—Specific aspects of physiological measurement analysis

Definitions

  • the present invention relates to a biological signal processing method and a biological signal processing apparatus for appropriately removing a noise component mixed in a biological signal obtained from an electrocardiogram waveform to improve the analysis accuracy of the biological signal.
  • the pulsating rhythm of the heart is known to fluctuate under the influence of the autonomic nerve, that is, the sympathetic nerve and the vagus nerve.
  • the autonomic nerve that is, the sympathetic nerve and the vagus nerve.
  • heart rate fluctuation respiratory sinus arrhythmia
  • the respiratory rate at this time can be estimated by analyzing the time interval of the heartbeat (R wave) time extracted from the electrocardiogram waveform, that is, the RR interval.
  • the influence of respiration also appears on the amplitude of the electrocardiogram waveform.
  • the influence on the electrocardiogram waveform amplitude is considered to be due to the fluctuation of the impedance viewed from the electrocardiogram measurement system due to the expansion and contraction of the lung and thorax along with the respiration.
  • FIG. 17 (a) is a diagram showing an example of an electrocardiogram waveform when a subject practices a breathing method called 4-7-8 breathing method (4 seconds inhaled-7 seconds in breath hold-8 seconds in exhalation).
  • FIG. 17A it can be easily seen that the height of the electrocardiogram waveform is fluctuating.
  • a diagram in which data of an RR interval, which is a time interval of an R wave and an immediately preceding R wave, is extracted from this electrocardiogram waveform and plotted in time series is shown in (b) of FIG.
  • the RR interval is longer at exhalation and shorter at inspiration (heart rate is slower at exhalation and faster at inspiration).
  • the RR interval fluctuates in synchronization with the electrocardiogram waveform.
  • (c) of FIG. 17 is a diagram in which data of the amplitude from the peak value of the R wave to the peak value of the S wave, that is, RS amplitude is extracted for each heartbeat and plotted in time series.
  • the RS amplitude can be acquired, for example, by searching for an electrocardiogram waveform in a certain section before and after the extracted R wave, and calculating the difference between the maximum value and the minimum value of the section.
  • (c) of FIG. 17 it can be seen that the time width of the above-mentioned inspiratory-respiratory-expiratory appears in the fluctuation pattern of the RS amplitude.
  • FIG. 18 shows a spectrum obtained by analyzing the data of the RR interval shown in (b) of FIG. 17 by the Maximum Entropy Method (MEM), and (b) of FIG. It is a figure which shows the spectrum which analyzed the data of RS amplitude shown to (c) of FIG. 17 by MEM.
  • MEM Maximum Entropy Method
  • noise when measuring an electrocardiogram waveform, noise may be added to the waveform.
  • noise due to body movement or the like is likely to be introduced.
  • Such noise may also be mixed into the RS amplitude.
  • an error in R-wave extraction is induced, and as a result, it becomes a factor causing improper data to enter into the RR interval or the like.
  • Patent Document 1 in the method of performing respiration estimation based on the peak value of the T wave of the electrocardiogram waveform, the case where the T wave is not generated or the peak value of the T wave is very small is detected by a predetermined threshold, An arrangement for calibrating smaller peak values is disclosed.
  • Patent Document 1 when large noise is superimposed on an electrocardiogram waveform, it is impossible to correct biological signals such as RS amplitude and RR interval.
  • Patent No. 5632570 gazette
  • the present invention has been made in view of the above-described points, and noise components mixed in biological signals such as RS amplitude and RR interval can be appropriately removed, and the accuracy of analysis of biological signals can be reduced.
  • An object of the present invention is to provide a biological signal processing method and a biological signal processing device that can be improved.
  • the biological signal processing method of the present invention comprises a first step of extracting a biological signal from an electrocardiogram waveform of a living body, and a second step of calculating averaged data using time-series data of the biological signal extracted in the first step. And determining whether the data of the biological signal extracted in the first step is appropriate for each data based on the averaged data calculated using the data of the biological signal at a time earlier than this data And a fourth step of either deleting or interpolating the data of the biological signal determined to be inappropriate in the third step.
  • the biological signal processing apparatus of the present invention uses a biological signal extraction unit configured to extract a biological signal from an electrocardiogram waveform of a living body, and time-series data of the biological signal extracted by the biological signal extraction unit.
  • Averaging processing unit configured to calculate averaged data, and whether the data of the biological signal extracted by the biological signal extraction unit is appropriate, using data of biological signal at a time earlier than this data
  • an abnormal value determination unit configured to make a determination for each data based on the averaged data calculated as described above, or deletion or interpolation of data of a biological signal determined to be inappropriate by the abnormal value determination unit
  • an abnormal value processing unit configured to perform the following.
  • whether or not the data of the biological signal extracted in the first step is appropriate is determined for each data based on the averaged data calculated using the data of the biological signal at a time earlier than this data.
  • the noise component mixed in the biological signal can be appropriately removed by either deleting or interpolating the data of the biological signal determined to be inappropriate, and the accuracy of analysis of the biological signal can be improved. it can.
  • the value is inherently variable, and in order to determine if inappropriate data may be mixed, it is not a fixed value but based on the averaged value of the previous values. It is preferable to treat.
  • FIG. 1 is a diagram showing an example of time-series data of an electrocardiogram waveform and an RR interval.
  • FIG. 2 is a diagram showing an example of time-series data of RS amplitude.
  • FIG. 3 is a block diagram showing the configuration of the biological signal processing apparatus according to the first embodiment of the present invention.
  • FIG. 4 is a flowchart for explaining the operation of the biological signal processing apparatus according to the first embodiment of the present invention.
  • FIG. 5 is a diagram showing an example of RS amplitude time-series data, RS amplitude averaged data, and RS amplitude data after interpolation.
  • FIG. 1 is a diagram showing an example of time-series data of an electrocardiogram waveform and an RR interval.
  • FIG. 2 is a diagram showing an example of time-series data of RS amplitude.
  • FIG. 3 is a block diagram showing the configuration of the biological signal processing apparatus according to the first embodiment of the present invention.
  • FIG. 4 is a flowchar
  • FIG. 6 is a diagram showing an example of time-series data of RS amplitude and data obtained by resampling time-series data of RS amplitude.
  • FIG. 7 is a diagram showing an example of a spectrum of RS amplitude.
  • FIG. 8 is a block diagram showing the configuration of a biological signal processing apparatus according to a second embodiment of the present invention.
  • FIG. 9 is a flowchart for explaining the operation of the biological signal processing apparatus according to the second embodiment of the present invention.
  • FIG. 10 is a diagram showing an example of an electrocardiogram waveform, time series data of an RR interval, averaged data of an RR interval, and data of an RR interval after interpolation.
  • FIG. 10 is a diagram showing an example of an electrocardiogram waveform, time series data of an RR interval, averaged data of an RR interval, and data of an RR interval after interpolation.
  • FIG. 11 is a view showing an example of averaged data of RR intervals according to the second embodiment of the present invention and data of RR intervals after interpolation.
  • FIG. 12 is a flow chart for explaining the operation of the biological signal processing apparatus according to the third embodiment of the present invention.
  • FIG. 13 is a diagram illustrating an example of first and second order differential values of time-series data of RS amplitude, and determination results of the change amount decrease determination unit.
  • FIG. 14 is a diagram showing an example of a spectrum obtained by frequency analysis of RS amplitude by fast Fourier transform in the fourth embodiment of the present invention.
  • FIG. 15 is a diagram showing an example of a spectrum in which RS amplitude is frequency analyzed by the maximum entropy method in the fourth embodiment of the present invention.
  • FIG. 12 is a flow chart for explaining the operation of the biological signal processing apparatus according to the third embodiment of the present invention.
  • FIG. 13 is a diagram illustrating an example of first and second order differential values of time-series data
  • FIG. 16 is a block diagram showing a configuration example of a computer for realizing the biological signal processing apparatus according to the first to fifth examples of the present invention.
  • FIG. 17 is a diagram showing an example of an electrocardiogram waveform, time series data of an RR interval, and time series data of an RS amplitude.
  • FIG. 18 shows a spectrum of an RR interval and a spectrum of an RS amplitude.
  • FIG. 1 (a) shows an example of an electrocardiogram waveform and time series data of RR intervals extracted from the electrocardiogram waveform
  • FIG. 1 (b) is an enlarged view of a part of FIG. 1 (a).
  • 100 indicates an electrocardiogram waveform (unit [ ⁇ V])
  • 101 indicates an RR interval (unit [ms]).
  • the horizontal axis of (a) of FIG. 1 and (b) of FIG. 1 is time.
  • FIG. 1 is an enlarged view of (a) of FIG. 1, in the heartbeat seen in the portion 102, the R wave overlaps the noise peak, and the S wave is on the noise slope I understand.
  • FIG. 2 is a diagram in which RS amplitude data is extracted for each heartbeat from the electrocardiogram waveform of (a) of FIG. 1 and plotted in time series.
  • the RS amplitude is calculated from the difference between the maximum value and the minimum value of the electrocardiogram waveform in an interval of ⁇ 25 ms of the heartbeat time.
  • the point of the portion 103 is affected by noise and has an inappropriate value.
  • FIG. 3 is a block diagram showing the configuration of the biological signal processing apparatus according to the present embodiment.
  • the biological signal processing apparatus averages the time series data of the biological signal extracted by the biological signal extraction unit 2 which extracts the biological signal from the electrocardiogram waveform measured by the electrocardiograph 1 and the biological signal extraction unit 2
  • the averaging processing unit 3 to be performed compares the data of the biological signal extracted by the biological signal extracting unit 2 with the averaged data calculated using data of the time before this data, and the biological signal extracting unit
  • the abnormal value determination unit 4 that determines whether the data of the biological signal extracted in 2 is appropriate and the data of the biological signal determined to be inappropriate by the abnormal value determination unit 4 are deleted or replaced with appropriate data
  • Value interpolation unit 5 that interpolates and the differentiation unit 6 that calculates the first and second derivative values of the biological signal after being processed by the abnormal value processing unit 5, first derivative and second derivative values
  • the frequency change analysis unit 7 performs frequency analysis on the time-series data of the biological signal acquired by the resampling unit
  • the resampling unit 8 samples the time-series data of the biological signal processed by the change amount decrease determination unit 7 and the abnormal value processing unit 5.
  • the electrocardiograph 1 measures an electrocardiogram waveform of a subject (living body) not shown.
  • a specific method of measuring an electrocardiogram waveform is a known technique, and thus the detailed description is omitted.
  • the biological signal extraction unit 2 extracts a biological signal (RS amplitude in the present embodiment) from the electrocardiogram waveform measured by the electrocardiograph 1 (step S1 in FIG. 4).
  • the biological signal extraction unit 2 of this embodiment corresponds to the maximum value (maximum value of the R wave) of the electrocardiogram waveform in a predetermined section (here, a section of ⁇ 25 ms of the heartbeat time) centered on the heartbeat time.
  • the minimum value of the electrocardiogram waveform (corresponding to the minimum value of the S wave) is extracted for each heartbeat as the RS amplitude.
  • the technology disclosed in Japanese Patent Laid-Open No. 2015-156936 can be used.
  • sampling data of two points of an electrocardiogram waveform straddling a predetermined threshold is detected between a representative point of R wave and a representative point of S wave existing after this point,
  • the time at which the straight line connecting the two points of sampling data intersects with the threshold is calculated as the heartbeat time.
  • this heartbeat time is taken as the time of the data of RS amplitude.
  • the abnormal value determination unit 4 performs averaging using data of the biological signal (RS amplitude) extracted by the biological signal extraction unit 2 and data of the biological signal (RS amplitude) up to the immediately preceding time by the averaging processing unit 3 The data is compared with each other to determine whether the data of the biological signal extracted by the biological signal extraction unit 2 is appropriate for each data (step S2 in FIG. 4).
  • the abnormal value determination unit 4 centers on the averaged data X ′ (i ⁇ 1) of the biological signal up to the immediately preceding time, as the value X (i) of the data of the biological signal at a certain time If it is within the range of the predetermined normal value, it is determined that the data X (i) of the biological signal is appropriate, and if the data X (i) is out of the range of the normal value, the data X (i) is not Determine that it is appropriate. In the present embodiment, the range of ⁇ 30% of the averaged data X ′ (i ⁇ 1) is taken as the range of the normal value.
  • the abnormal value determination unit 4 determines whether the data X (i) of the biological signal is appropriate by using the averaged data X ′ (i-1) of the biological signal at the past time, the biological signal extraction is performed. The determination process is not performed on the data of the first biological signal extracted by the unit 2. The abnormal value determination unit 4 performs the determination process on the second and subsequent biological signal data extracted by the biological signal extraction unit 2.
  • the averaging processing unit 3 performs averaging processing on time-series data of the biological signal (RS amplitude) extracted by the biological signal extraction unit 2 (step S3 in FIG. 4). Assuming that the value of the i-th biological signal before the averaging process is X (i), and the value obtained by averaging the biological signals up to the i-th is X ′ (i), the averaging processor 3 gives The averaging process of the biological signal is performed.
  • X '(i) r * X (i) + (1-r) * X' (i-1) (1)
  • the averaging processing unit 3 does not take in the data of the biological signal (RS amplitude) determined to be inappropriate by the abnormal value determination unit 4 in the averaging processing, in order to prevent the erroneous value from mixing into the averaging processing. Let's do it. For example, when the data X (i) of the biological signal is determined to be inappropriate, the averaged data X '(i-1) of the biological signal up to the immediately preceding time is averaged as it is without using the data X (i) What is necessary is to set it as the conversion data X '(i). This makes it possible to further stabilize the transition of the value of the averaged data.
  • the abnormal value processing unit 5 interpolates by replacing the data of the biological signal (RS amplitude) determined as inappropriate by the abnormal value determination unit 4 with appropriate data (step S4 in FIG. 4).
  • a method of interpolation at this time there is linear interpolation in which data of a biological signal determined to be inappropriate is interpolated using appropriate data before and after that.
  • an interpolation method such as spline interpolation may be used instead of linear interpolation.
  • the biological signal extraction unit 2, the abnormal value determination unit 4, the averaging processing unit 3, and the abnormal value processing unit 5 perform steps S1 to S4 until an instruction to end the measurement is given from the subject (YES in FIG. 4).
  • the process is performed every predetermined cycle (for example, every sampling of the electrocardiograph 1).
  • the display unit 10 displays time-series data of the biological signal (RS amplitude) processed by the abnormal value processing unit 5 (step S6 in FIG. 4).
  • FIG. 5 is a view showing an example of time series data of RS amplitude, data obtained by averaging processing section 3 averaging time series data of RS amplitude, and data of RS amplitude interpolated by abnormal value processing section 5.
  • Circles 50 in FIG. 5 indicate time series data of the same RS amplitude as in FIG. 2, broken lines 51 indicate data obtained by averaging the time series data, and asterisks 52 indicate data of RS amplitudes after interpolation. . According to FIG. 5, it can be seen that data of inappropriate RS amplitude affected by noise can be properly interpolated.
  • the averaging processing unit 3 targets the averaging processing with the data of the biological signal extracted by the biological signal extraction unit 2 and not the data interpolated by the abnormal value processing unit 5 but the interpolated data is , Will not be incorporated into the subsequent averaging process.
  • the reason for not using the data after interpolation for the averaging process is that the data after interpolation is estimated based on the averaged data, and the data after interpolation is used to derive the data itself It is not appropriate to reincorporate it into the value that had been.
  • FIG. 6 is a diagram showing an example of time series data of RS amplitude and data obtained by resampling time series data of RS amplitude. Circles 60 in FIG. 6 indicate time series data of RS amplitudes, and asterisks 61 indicate data resampled by the resampling unit 8 after interpolation by the abnormal value processing unit 5.
  • the frequency analysis unit 9 analyzes the frequency of the time series data of the biological signal (RS amplitude) acquired by the resampling unit 8 by the fast Fourier transform or the maximum entropy method (MEM) to obtain the spectrum of the biological signal (step S8 in FIG. 4). ).
  • the display unit 10 displays the spectrum of the frequency analysis result by the frequency analysis unit 9 (step S9 in FIG. 4).
  • FIG. 7 is a diagram showing an example of a spectrum analyzed by MEM.
  • a solid line 70 in FIG. 7 indicates a spectrum obtained by frequency analysis of the time-series data (61 in FIG. 6) of the RS amplitude acquired by the resampling unit 8 using MEM.
  • a broken line 71 indicates a spectrum obtained by performing frequency analysis on time series data (60 in FIG. 6) of the RS amplitude before the interpolation processing by the MEM.
  • the spectrum is different between the spectrum obtained from the RS amplitude data before interpolation processing and the spectrum obtained from the RS amplitude data after interpolation of inappropriate data, and the RS amplitude data before interpolation processing is used. It can be seen that the spectrum obtained includes the component or distribution (72, 73, 74 in FIG. 7) of the spectrum due to the inappropriate data. It is considered that the spectrum obtained from the data of RS amplitude after interpolation of inappropriate data shows a result that correctly reflects the condition of the living body.
  • noise components mixed in a biological signal such as RS amplitude can be appropriately removed, and the accuracy of analysis of the biological signal can be improved.
  • interpolation is performed by replacing the data of the biological signal determined as inappropriate by the abnormal value determination unit 4 with data that is likely to be inappropriate, but the present invention is not limited to this. Data of a biological signal determined to be inappropriate by the value determination unit 4 may be deleted (data loss).
  • FIG. 8 is a block diagram showing the configuration of a biological signal processing apparatus according to this embodiment.
  • the biological signal processing apparatus according to this embodiment includes a biological signal extraction unit 2, an averaging processing unit 3 a, an abnormal value determination unit 4, an abnormal value processing unit 5, a differentiation unit 6, and a change amount reduction determination unit 7. , A resampling unit 8, a frequency analysis unit 9, and a display unit 10.
  • a biological signal extraction unit 2 an averaging processing unit 3 a
  • an abnormal value determination unit 4 an abnormal value processing unit 5, a differentiation unit 6, and a change amount reduction determination unit 7.
  • a resampling unit 8 a frequency analysis unit 9
  • a display unit 10 In the present embodiment, an example will be described in which inappropriate data is removed from time-series data of an RR interval.
  • the biological signal extraction unit 2 extracts a biological signal (RR interval in this embodiment) from the electrocardiogram waveform measured by the electrocardiograph 1 (step S10 in FIG. 9).
  • a biological signal for example, the document "ECG Implementation on the TMS320C5515 DSP Medical Development Kit (MDK) with the ADS1298 ECG-FE", Texas Instruments Incorporated, ⁇ http://www.ti.com/lit The technology disclosed in /an/sprabj1/sprabj1.pdf>, 2011 "can be used.
  • the RR interval is obtained based on the change in the time difference value of the electrocardiogram waveform.
  • the heartbeat time may be determined using the technique disclosed in Japanese Patent Application Laid-Open No. 2015-156936, and the interval of the heartbeat time may be taken as the RR interval.
  • the abnormal value determination unit 4 uses the data of the biological signal (RR interval) extracted by the biological signal extraction unit 2 and the data of the biological signal (RR interval) up to the time immediately before the averaging processing unit 3a. The calculated averaged data is compared, and it is determined for each data whether the data of the biological signal extracted by the biological signal extraction unit 2 is appropriate (FIG. 9, step S11).
  • the value X (i) of the data of the biological signal at a certain time is a predetermined multiple of the averaged data X ′ (i-1) of the biological signal up to the immediately preceding time
  • the value exceeds 1.35 times it is determined that the data X (i) is inappropriate, and the data X (i) is not more than a predetermined number of times the averaged data X '(i-1).
  • the range equal to or less than a predetermined multiple of the averaged data X ′ (i ⁇ 1) is the range of the normal value.
  • the abnormal value determination unit 4 performs the determination process on the second and subsequent biological signal data extracted by the biological signal extraction unit 2.
  • the averaging processing unit 3a performs averaging processing on the time-series data of the biological signal (RR interval) extracted by the biological signal extraction unit 2 (steps S12 and S13 in FIG. 9). The averaging process of this embodiment will be described later.
  • the abnormal value processing unit 5 interpolates by replacing data of the biological signal (RR interval) determined by the abnormal value determination unit 4 as inappropriate with appropriate data (step S14 in FIG. 9).
  • the abnormal value processing unit 5 determines that the abnormal value determination unit 4 generates inappropriate data in the biological signal (RR interval), the time at which the inappropriate data was generated and the time immediately before the time And estimate the number of data to be inserted (corresponding to the heart rate).
  • the abnormal value processing unit 5 sets a time interval (t 2 ⁇ t 1 ) between a time t 2 when inappropriate data occurs in the biological signal (RR interval) and a time t 1 of data immediately before that. ), and by dividing the value of the average data of R-R interval until time t 1, determining a number N of data to be inserted between the time t 2 and time t 1.
  • the abnormal value processing unit 5 equally divides the time interval (t 2 ⁇ t 1 ) between the time t 2 and the time t 1 by the determined number N of data into the time t 2 and the time t 1 . Calculated as the probable value of the RR interval to be inserted between.
  • the abnormal value processing unit 5 between times t 2 to inappropriate data occurs in R-R interval between the time t 1 immediately before, insert a plausible value of the number N by R-R intervals
  • the RR interval can be interpolated.
  • FIG. 10A shows an example of an electrocardiogram waveform and time series data of an RR interval extracted from the electrocardiogram waveform.
  • 80 indicates an electrocardiogram waveform (unit [ ⁇ V])
  • 81 indicates an RR interval (unit [ms]).
  • the RR interval largely fluctuates due to the respiratory motion of the living body at around 16:31:40.
  • FIG. 10B is a diagram showing an example of data obtained by averaging the time-series data of the RR interval shown in FIG. 10A and data interpolated by the abnormal value processing unit 5.
  • the broken line 82 in (b) of FIG. 10 indicates data obtained by averaging the time series data of the RR interval according to the method described in the first embodiment, and the * mark 83 indicates the data after interpolation.
  • averaging is performed using the values of the RR intervals themselves. Also, in order to prevent erroneous values from being included in the averaging, do not take into the averaging in the case of values separated by a fixed percentage (more than ⁇ 30% in this example) compared to the previous averaged values. And, the value of the RR interval showing the fluctuation derived from the original living condition is regarded as an inappropriate value and is not taken into the averaging, so that it becomes a form of sticking to the past value. It turns out that it leads to an incorrect result, such as continuing unnecessary interpolation.
  • the averaging processing unit 3a performs averaging on the values based on the reciprocal of the RR intervals. And calculating the averaged data of the RR interval from the reciprocal of the value obtained by this averaging process, for each piece of data.
  • the processing target up to the i-th A value C ′ (i) obtained by averaging the data can be obtained by the following equation similar to the equation (1).
  • C '(i) r * C (i) + (1-r) * C' (i-1) (2)
  • the processing target data C (i) is inappropriate and is not taken into the averaging process.
  • the value C ′ (i ⁇ 1) obtained by averaging the processing target data up to the immediately preceding time may be taken as C ′ (i).
  • the range of ⁇ 30% of the averaged data C ′ (i ⁇ 1) is taken as the normal value range.
  • steps S14, S15, S16, S17, S18, and S19 in FIG. 9 are the same as those in steps S4, S5, S6, S7, S8, and S9 in FIG.
  • FIG. 11 is a diagram showing an example of data obtained by averaging the time-series data of the RR interval shown in FIG. 10A by the method of the present embodiment, and data interpolated by the abnormal value processing unit 5. Circles and * marks 91 in FIG. 11 indicate time series data of the RR interval and data of the RR interval after interpolation, and a broken line 92 indicates averaged data of the RR interval.
  • the variation range of the RR interval of the steepest variation is about 650 ms ⁇ 900 ms (42% variation), while the variation range of the heart rate proportional to the reciprocal of the RR interval is about 92 bpm ⁇ 67 bpm (The fluctuation amount is 27%). Therefore, according to the method of the present embodiment, even in the place where the RR interval fluctuates, those fluctuations are taken into the averaging, and the data of the interpolated RR interval is R- before the interpolation process. There is no deviation from the R interval data. That is, by using the heart rate, the scale of fluctuation can be suppressed, so that the averaging process is stabilized and the process on the data sequence of the biological signal is performed correctly.
  • the fluctuation range of the averaged data is performed by performing averaging processing on the value based on the reciprocal of the data value rather than the data value itself in a certain numerical range.
  • the averaging process can be stabilized. According to this embodiment, in the time-series data of the biological signal, it is possible to remove inappropriate data caused by noise etc. and restore it in a plausible manner, thereby linking it to more accurate analysis of the condition of the living body. it can.
  • the differentiating unit 6 calculates first and second derivative values of the biological signal (RS amplitude) processed by the abnormal value processing unit 5 (step S20 in FIG. 12).
  • the first-order differential value f ′ (t k ) can be expressed as equation (3).
  • f '(t k ) ⁇ f (t k + 1 ) -f (t k ) ⁇ / (t k + 1 -t k ) (3)
  • the second-order differential value f ′ ′ (t k ) can be expressed as in equation (4).
  • f "(t k) ⁇ f (t k + 1) -2f (t k) + f (t k-1) ⁇ / (t k + 1 -t k) 2 ... (4)
  • the change amount decrease determination unit 7 determines that the state of the first-order differential value calculated by the differential unit 6 and the second-order differential value within a predetermined range centered on 0 continues for a predetermined time or more.
  • the display unit 10 displays the determination result of the change amount reduction determination unit 7 (step S22 in FIG. 12). This makes it possible to warn that the fluctuation of the biological signal is low.
  • a warning message may be displayed on the display unit 10 as a notification means, a notification means comprising a light emitter such as an LED may be separately provided and blinked, or a notification means such as a speaker may be provided. It may be separately provided to issue a warning by voice.
  • FIG. 13 is a diagram showing an example of the first-order differential value and the second-order differential value of the time-series data of RS amplitude shown in FIG. 5 and the determination result of the change amount reduction determination unit 7.
  • the ⁇ mark 110 in FIG. 13 indicates the first-order derivative value
  • the ⁇ mark 111 indicates the second-order derivative value
  • the solid line 112 indicates the determination result of the change amount reduction determination unit 7.
  • the biological signal fluctuation is low and the determination result is “1”.
  • at least one of the first derivative and the second derivative is out of the range of ⁇ 80 ⁇ V, or the first derivative and the second derivative are both in the range of ⁇ 80 ⁇ V If the duration of time is less than 3 seconds, it is determined that the fluctuation of the biological signal is normal and the determination result is "0". Thus, in this embodiment, it can be determined whether the breathing movement is changing or stopping with time, and the breathing movement of the living body can be monitored.
  • the frequency analysis unit 9 of the biological signal processing apparatus frequency analyzes time series data of the biological signal (RS amplitude) acquired by the resampling unit 8 to obtain the spectrum of the biological signal ( The display unit 10 displays this spectrum (step S9 in FIG. 4).
  • FIG. 14 is a diagram showing an example of a spectrum obtained by frequency analysis of RS amplitude by fast Fourier transform.
  • the horizontal axis is the frequency
  • the vertical axis is the spectral energy density normalized so that the maximum peak value is 1.
  • FIG. 15 is a diagram showing an example of a spectrum analyzed by the maximum entropy method using data of the same RS amplitude as in the case of FIG.
  • the maximum entropy method is different from fast Fourier transform, and frequency resolution is independent of measurement time, so high frequency resolution can be obtained without increasing measurement time.
  • the period of normal breathing is 3 to 4 seconds, and if the measurement time is lengthened to increase the frequency resolution with fast Fourier transform, data for multiple cycles will be included, and the obtained frequency distribution will be those statistics. In the end, information on a single respiratory activity is buried. However, since maximum resolution can be used to obtain high resolution in a short measurement time, it is possible to obtain frequency information for a single non-statistical breath.
  • the calculation of the frequency by the maximum entropy method is carried out by the following method (refer to the document, "N. Shigeo," Processing of waveform data for scientific measurement ", CQ publication, p. 173-174, 1986").
  • the maximum entropy method there are the Burg method and the Yule-Walker method.
  • the Burg method is taken as an example.
  • the spectral energy density S ( ⁇ ) to be determined is given by the following equation.
  • ⁇ t is a sampling rate, which is 0.2 seconds in the example of FIG.
  • the model order m takes an arbitrary integer, but here, the maximum model order has a value of 16 or more (see “Inoue Hiroshi,” Circular disease and autonomic nervous function, Medical School, p. 85-86, 2010 ". ) To 20, and the order m is 1 to 20.
  • N is the number of RS amplitude data, which is 32 here.
  • initial values of the coefficients b mi and b ' mi in the equations (7) to (9) are expressed by the following equations.
  • x i represents the i-th data out of N data.
  • the regression relationship equation of equations (13) and (14) is used to find the coefficient a mi of the autoregressive model and the variance P m of the prediction error.
  • Equation (14) P 0 used in equation (14) can be obtained by equation (15), assuming that the average value of data of RS amplitude is x ave .
  • the statistic Q m is calculated from Q 1 to Q 20 using equation (16), and m giving the smallest Q m among them is the model order m used in equation (6).
  • the frequency distribution can be obtained by Equation (6).
  • the determination process is performed using the above example. It is not limited.
  • the range of averaged data X ′ (i ⁇ 1) ⁇ ⁇ ( ⁇ is a specified value) may be set as the range of normal values.
  • the abnormal value determination unit 4 of the first to fourth embodiments uses the data X (i) of the biosignal to be determined and the biosignal data of a time earlier than this data (up to the immediately preceding time).
  • the variance ⁇ 2 obtained from the averaged data calculated using the biological signal data) may be calculated.
  • the abnormal value determination unit 4 determines that the variance ⁇ 2 is out of the range of the predetermined normal value based on the variance ⁇ p 2 obtained from the averaged data calculated using the data of the biological signal at the past time.
  • the data X (i) of the biological signal may be determined to be inappropriate. For example, if the range of 2 ⁇ p 2 or less is a normal range, and the variance ⁇ 2 exceeds 2 ⁇ p 2 , it may be determined that the data X (i) of the biological signal is inappropriate.
  • the biological signal processing apparatus described in the first to fifth embodiments can be realized by a computer having a central processing unit (CPU), a storage device, and an interface, and a program for controlling these hardware resources.
  • the computer includes a CPU 40, a storage device 41, and an interface device (hereinafter abbreviated as I / F) 42.
  • I / F interface device
  • the electrocardiograph 1 and the like are connected to the I / F 42.
  • a program for realizing the biological signal processing method of the present invention is provided in the state of being recorded on a recording medium such as a flexible disk, a CD-ROM, a DVD-ROM, a memory card, etc.
  • Stored in The CPU 40 executes the processing described in the first to fifth embodiments in accordance with the program stored in the storage device 41.
  • the present invention can be applied to a technique for analyzing a biological signal obtained from an electrocardiogram waveform.
  • SYMBOLS 1 electrocardiograph, 2 ... biomedical signal extraction part, 3, 3 ... averaging processing part, 4 ... abnormal value determination part, 5 ... abnormal value processing part, 6 ... differentiation part, 7 ... change amount reduction determination part, 8 ... Resampling unit, 9 ... Frequency analysis unit, 10 ... Display unit, 30 ... Inverse number average processing unit, 31 ... Averaged data calculation unit.

Abstract

生体信号処理装置は、心電計（１）が測定した心電図波形から生体信号を抽出する生体信号抽出部（２）と、生体信号抽出部（２）で抽出された生体信号の時系列データを用いて平均化データを算出する平均化処理部（３）と、生体信号抽出部（２）で抽出された生体信号のデータが適切かどうかを、このデータよりも過去の時刻の生体信号のデータを用いて算出された平均化データに基づいてデータ毎に判定する異常値判定部（４）と、異常値判定部（４）で不適切と判定された生体信号のデータの削除および補間のいずれかを行う異常値処理部（５）とを備えている。

Description

生体信号処理方法および生体信号処理装置
　本発明は、心電図波形から得られる生体信号に混入したノイズ成分を適切に除去して、生体信号の解析精度を向上させる生体信号処理方法および生体信号処理装置に関するものである。
　心臓の拍動リズムは、自律神経すなわち交感神経・迷走神経の影響を受けて変動することが知られている。例えば、安静かつ寛いだ状態では、迷走神経が亢進し、呼吸に伴う心拍変動（呼吸性洞性不整脈）が顕著にみられるようになる。このときの呼吸数は、心電図波形から抽出される心拍（Ｒ波）時刻の時間間隔、すなわちＲ－Ｒ間隔を分析することで、推定することができる。
　一方、心電図波形の振幅にも呼吸の影響が現れる。この心電図波形振幅への影響は、呼吸に伴って肺や胸郭が伸縮することで、心電図測定系から見たインピーダンスが変動することによると考えられる。
　図１７の（ａ）は、被験者が４－７－８呼吸法（吸気４秒－息止め７秒－呼気８秒）と呼ばれる呼吸方法を実践した際の心電図波形の例を示す図である。図１７の（ａ）によると、心電図波形の高さが変動していることが容易に見て取れる。この心電図波形から、Ｒ波と１つ前のＲ波の時間間隔であるＲ－Ｒ間隔のデータを抽出して時系列にプロットした図を図１７の（ｂ）に示す。Ｒ－Ｒ間隔は、呼気時に長くなり、吸気時に短くなる（心拍は、呼気時に遅くなり、吸気時に速くなる）。図１７の（ｂ）によると、Ｒ－Ｒ間隔が、心電図波形と同期して変動していることが分かる。
　また、図１７の（ｃ）は、Ｒ波のピーク値からＳ波のピーク値までの振幅、すなわちＲＳ振幅のデータを心拍毎に抽出して時系列にプロットした図である。ＲＳ振幅は、例えば、抽出したＲ波に対して、前後一定の区間の心電図波形を探索し、その区間の最大値と最小値の差をとる、といった方法で取得できる。図１７の（ｃ）では、上述の吸気－息止め－呼気の時間幅がＲＳ振幅の変動パターンに現れていることが分かる。
　図１８の（ａ）は、図１７の（ｂ）に示したＲ－Ｒ間隔のデータを最大エントロピー法（ＭＥＭ：Maximum Entropy Method）で解析したスペクトルを示す図、図１８の（ｂ）は、図１７の（ｃ）に示したＲＳ振幅のデータをＭＥＭで解析したスペクトルを示す図である。ＭＥＭを用いることにより、６０秒という短いデータからでも滑らかなスペクトル曲線が得られる。図１８の（ａ）、図１８の（ｂ）の例では、ＲＳ振幅から得られたスペクトルは、Ｒ－Ｒ間隔から得られたスペクトルよりも明確に、１９秒周期≒０．０５２Ｈｚのピークを示している。このように、Ｒ－Ｒ間隔と併せて、ＲＳ振幅からも呼吸についての情報を得ることができる。
　ところで、心電図波形を計測する際、波形にノイズが加わることがある。特に、携帯型の装置やウエアラブルな装置を用いて日常生活の中での心電図波形を取得する場合には、体動などによるノイズが入りやすい。そうしたノイズは、ＲＳ振幅にも混入する可能性がある。また、Ｒ波抽出の誤りを誘発し、その結果、Ｒ－Ｒ間隔などに不適切なデータを入り込ませる要因になる。
　特許文献１には、心電図波形のＴ波のピーク値を基に呼吸推定を行う方法において、Ｔ波が発生しない場合もしくはＴ波のピーク値が非常に小さい場合を所定閾値により検出し、所定閾値より小さいピークの値を校正する構成が開示されている。しかしながら、特許文献１に開示された技術では、心電図波形に大きなノイズが重畳した場合、ＲＳ振幅やＲ－Ｒ間隔などの生体信号を補正することは不可能であった。
　ＲＳ振幅やＲ－Ｒ間隔などの生体信号を基に、被験者の呼吸等に関する分析を行う場合、それら生体信号の、呼吸等に同期した変動成分に着目することになる。そこで、ＭＥＭや高速フーリエ変換（ＦＦＴ：Fast Fourier Transform）といったスペクトル解析などの手法を生体信号に適用し、呼吸等に関する情報を抽出するといった方法が採られる。ところが、生体信号のデータ列に、ノイズ等に由来する不適切な値が含まれていると、解析の結果が本来のものとは乖離してしまうという問題点があった。
特許第５６３２５７０号公報
　本発明は、上記のような点に鑑みてなされたものであり、ＲＳ振幅やＲ－Ｒ間隔などの生体信号に混入したノイズ成分を適切に除去することができ、生体信号の解析の精度を向上させることができる生体信号処理方法および生体信号処理装置を提供することを目的とする。
　本発明の生体信号処理方法は、生体の心電図波形から生体信号を抽出する第１のステップと、この第１のステップで抽出した生体信号の時系列データを用いて平均化データを算出する第２のステップと、前記第１のステップで抽出した生体信号のデータが適切かどうかを、このデータよりも過去の時刻の生体信号のデータを用いて算出した前記平均化データに基づいてデータ毎に判定する第３のステップと、この第３のステップで不適切と判定した生体信号のデータの削除および補間のいずれかを行う第４のステップとを含むことを特徴とするものである。
　また、本発明の生体信号処理装置は、生体の心電図波形から生体信号を抽出するように構成された生体信号抽出部と、この生体信号抽出部で抽出された生体信号の時系列データを用いて平均化データを算出するように構成された平均化処理部と、前記生体信号抽出部で抽出された生体信号のデータが適切かどうかを、このデータよりも過去の時刻の生体信号のデータを用いて算出された前記平均化データに基づいてデータ毎に判定するように構成された異常値判定部と、この異常値判定部で不適切と判定された生体信号のデータの削除および補間のいずれかを行うように構成された異常値処理部とを備えることを特徴とするものである。
　本発明によれば、第１のステップで抽出した生体信号のデータが適切かどうかを、このデータよりも過去の時刻の生体信号のデータを用いて算出した平均化データに基づいてデータ毎に判定し、不適切と判定した生体信号のデータの削除および補間のいずれかを行うことにより、生体信号に混入したノイズ成分を適切に除去することができ、生体信号の解析の精度を向上させることができる。生体信号においては本来その値は変動するものであり、不適切なデータが混入しうる場合にそれを判別するためには、固定的な値ではなく、それまでの値の平均化した値に基づいて処理することが好ましい。
図１は、心電図波形およびＲ－Ｒ間隔の時系列データの例を示す図である。 図２は、ＲＳ振幅の時系列データの例を示す図である。 図３は、本発明の第１の実施例に係る生体信号処理装置の構成を示すブロック図である。 図４は、本発明の第１の実施例に係る生体信号処理装置の動作を説明するフローチャートである。 図５は、ＲＳ振幅の時系列データ、ＲＳ振幅の平均化データ、および補間後のＲＳ振幅のデータの例を示す図である。 図６は、ＲＳ振幅の時系列データ、およびＲＳ振幅の時系列データを再サンプリングしたデータの例を示す図である。 図７は、ＲＳ振幅のスペクトルの例を示す図である。 図８は、本発明の第２の実施例に係る生体信号処理装置の構成を示すブロック図である。 図９は、本発明の第２の実施例に係る生体信号処理装置の動作を説明するフローチャートである。 図１０は、心電図波形、Ｒ－Ｒ間隔の時系列データ、Ｒ－Ｒ間隔の平均化データ、および補間後のＲ－Ｒ間隔のデータの例を示す図である。 図１１は、本発明の第２の実施例によるＲ－Ｒ間隔の平均化データ、および補間後のＲ－Ｒ間隔のデータの例を示す図である。 図１２は、本発明の第３の実施例に係る生体信号処理装置の動作を説明するフローチャートである。 図１３は、ＲＳ振幅の時系列データの１階微分値および２階微分値と、変化量低下判定部の判定結果の例を示す図である。 図１４は、本発明の第４の実施例においてＲＳ振幅を高速フーリエ変換により周波数解析したスペクトルの例を示す図である。 図１５は、本発明の第４の実施例においてＲＳ振幅を最大エントロピー法により周波数解析したスペクトルの例を示す図である。 図１６は、本発明の第１～第５の実施例に係る生体信号処理装置を実現するコンピュータの構成例を示すブロック図である。 図１７は、心電図波形、Ｒ－Ｒ間隔の時系列データおよびＲＳ振幅の時系列データの例を示す図である。 図１８は、Ｒ－Ｒ間隔のスペクトルおよびＲＳ振幅のスペクトルを示す図である。
［第１の実施例］
　図１ないし図７を参照して本発明の第１の実施例について説明する。図１の（ａ）は心電図波形、および心電図波形から抽出したＲ－Ｒ間隔の時系列データの例を示す図であり、図１の（ｂ）は図１の（ａ）の一部を拡大した図である。図１の（ａ）、図１の（ｂ）において、１００は心電図波形（単位［μＶ］）を示し、１０１はＲ－Ｒ間隔（単位［ｍｓ］）を示している。図１の（ａ）、図１の（ｂ）の横軸は時間である。
　心電図波形にはノイズが重畳しているものの、Ｒ波に相当する棘波が確認でき、Ｒ－Ｒ間隔は正しく抽出されている。図１の（ａ）を拡大した図１の（ｂ）によると、１０２の部分に見られる心拍では、Ｒ波がノイズのピークと重なっており、Ｓ波はノイズのスロープにかかっていることが分かる。
　図２は、図１の（ａ）の心電図波形からＲＳ振幅のデータを心拍毎に抽出して時系列にプロットした図である。ここでは、ＲＳ振幅は、心拍時刻の±２５ｍｓの区間での心電図波形の最大値と最小値の差から算出している。１０３の部分の点はノイズの影響を受けたものであり、不適切な値になっている。
　図３は本実施例に係る生体信号処理装置の構成を示すブロック図である。生体信号処理装置は、心電計１が測定した心電図波形から生体信号を抽出する生体信号抽出部２と、生体信号抽出部２で抽出された生体信号の時系列データに対して平均化処理を行う平均化処理部３と、生体信号抽出部２で抽出された生体信号のデータとこのデータよりも前の時刻のデータを用いて算出された平均化データとを比較して、生体信号抽出部２で抽出された生体信号のデータが適切かどうかを判定する異常値判定部４と、異常値判定部４で不適切と判定された生体信号のデータを削除するか、もしくは適切なデータに置き換えて補間する異常値処理部５と、異常値処理部５によって処理された後の生体信号の１階微分値と２階微分値を算出する微分部６と、１階微分値と２階微分値を基に生体信号の変動が低いかどうかを判定する変化量低下判定部７と、異常値処理部５によって処理された生体信号の時系列データをサンプリングするリサンプリング部８と、リサンプリング部８が取得した生体信号の時系列データを周波数解析して生体信号のスペクトルを求める周波数解析部９と、異常値処理部５によって処理された後の生体信号、変化量低下判定部７の判定結果および周波数解析部９の周波数解析結果を表示する表示部１０とを備えている。
　次に、図４を参照して本実施例の生体信号処理装置の動作を説明する。本実施例では、図２に示したＲＳ振幅の時系列データから不適切なデータを取り除く例について説明する。なお、微分部６と変化量低下判定部７の動作については他の実施例で説明する。
　心電計１は、図示しない被験者（生体）の心電図波形を測定する。心電図波形の具体的な測定方法は周知の技術であるので、詳細な説明は省略する。
　生体信号抽出部２は、心電計１が測定した心電図波形から生体信号（本実施例ではＲＳ振幅）を抽出する（図４ステップＳ１）。本実施例の生体信号抽出部２は、上記のとおり、心拍時刻を中心とする所定の区間（ここでは心拍時刻の±２５ｍｓの区間）での心電図波形の最大値（Ｒ波の最大値に相当）と心電図波形の最小値（Ｓ波の最小値に相当）との差をＲＳ振幅として心拍毎に抽出する。
　なお、心拍時刻を求める方法としては、例えば特開２０１５－１５６９３６号公報に開示された技術を用いることができる。この文献に開示された技術では、Ｒ波の代表点とこの点の後に存在するＳ波の代表点との間で、所定の閾値を跨ぐ心電図波形の２点のサンプリングデータを検出して、この２点のサンプリングデータを結ぶ直線が前記閾値と交わる時刻を心拍時刻として算出する。また、この心拍時刻を、ＲＳ振幅のデータの時刻とする。
　異常値判定部４は、生体信号抽出部２が抽出した生体信号（ＲＳ振幅）のデータと平均化処理部３が直前の時刻までの生体信号（ＲＳ振幅）のデータを用いて算出した平均化データとを比較し、生体信号抽出部２が抽出した生体信号のデータが適切かどうかをデータ毎に判定する（図４ステップＳ２）。
　具体的には、異常値判定部４は、ある時刻での生体信号のデータの値Ｘ（ｉ）が、直前の時刻までの生体信号の平均化データＸ´（ｉ－１）を中心とする所定の正常値の範囲内の場合、生体信号のデータＸ（ｉ）は適切であると判定し、データＸ（ｉ）が前記正常値の範囲から外れている場合、データＸ（ｉ）は不適切であると判定する。本実施例では、平均化データＸ´（ｉ－１）の±３０％の範囲を正常値の範囲とする。なお、異常値判定部４は、過去の時刻の生体信号の平均化データＸ´（ｉ－１）を用いて、生体信号のデータＸ（ｉ）が適切かどうかを判定するので、生体信号抽出部２が抽出した１個目の生体信号のデータについては判定処理を実施しない。異常値判定部４が判定処理を実施するのは、生体信号抽出部２が抽出した２個目以降の生体信号のデータについてである。
　次に、平均化処理部３は、生体信号抽出部２が抽出した生体信号（ＲＳ振幅）の時系列データに対して平均化処理を行う（図４ステップＳ３）。平均化処理部３は、平均化処理前のｉ番目の生体信号の値をＸ（ｉ）、ｉ番目までの生体信号を平均化した値をＸ´（ｉ）としたとき、次式のようにして生体信号の平均化処理を行う。
　Ｘ´（ｉ）＝ｒ×Ｘ（ｉ）＋（１－ｒ）×Ｘ´（ｉ－１） ・・・（１）
　式（１）におけるｒは予め定められた係数である。係数ｒの値を小さくするほど、生体信号のデータ列の値の細かな変動が抑えられるものの、一方で生体信号の大まかな変化に追随し難くなるので、このような点を鑑みて、例えばｒ＝０．２などとすると、生体信号の瞬時の変動が抑えられ、適度に平均化された生体信号のデータ列が得られる。
　ただし、平均化処理部３は、誤った値が平均化処理に混入することを防ぐために、異常値判定部４が不適切と判定した生体信号（ＲＳ振幅）のデータを平均化処理に取り込まないようにする。例えば生体信号のデータＸ（ｉ）が不適切と判定された場合、データＸ（ｉ）を使用せずに、直前の時刻までの生体信号の平均化データＸ´（ｉ－１）をそのまま平均化データＸ´（ｉ）とすればよい。これにより、平均化データの値の推移をより安定させることができる。
　異常値処理部５は、異常値判定部４が不適切と判定した生体信号（ＲＳ振幅）のデータを適切なデータに置き換えて補間する（図４ステップＳ４）。このときの補間の方法としては、不適切と判定された生体信号のデータを、その前後の適切なデータを用いて補間する線形補間がある。また、線形補間に代わってスプライン補間などの補間方法を用いてもよい。
　生体信号抽出部２と異常値判定部４と平均化処理部３と異常値処理部５とは、例えば被験者から測定終了の指示があるまで（図４ステップＳ５においてＹＥＳ）、ステップＳ１～Ｓ４の処理を所定の周期毎（例えば心電計１のサンプリング毎）に行う。
　表示部１０は、異常値処理部５によって処理された生体信号（ＲＳ振幅）の時系列データを表示する（図４ステップＳ６）。
　図５はＲＳ振幅の時系列データ、このＲＳ振幅の時系列データを平均化処理部３が平均化したデータ、および異常値処理部５が補間したＲＳ振幅のデータの例を示す図である。図５の○印５０は図２と同じＲＳ振幅の時系列データを示し、破線５１はこの時系列データを平均化したデータを示し、＊印５２は補間後のＲＳ振幅のデータを示している。図５によれば、ノイズの影響を受けた不適切なＲＳ振幅のデータを適切に補間できていることが分かる。
　なお、平均化処理部３が平均化処理の対象とするのは、生体信号抽出部２が抽出した生体信号のデータであって、異常値処理部５が補間したデータではなく、補間したデータが、その後の平均化処理に組み込まれることはない。補間後のデータを平均化処理に使用しない理由は、補間後のデータは、平均化されたデータに基づいて推定されたものであり、その補間後のデータを、このデータ自身を導くのに用いた値に再び組み込むことは適切ではないからである。
　次に、リサンプリング部８は、異常値処理部５によって処理された生体信号（ＲＳ振幅）の時系列データを、心電計１のサンプリング周波数よりも遅いサンプリング周波数（例えば１秒間隔）でサンプリングする（図４ステップＳ７）。
　図６はＲＳ振幅の時系列データ、およびＲＳ振幅の時系列データを再サンプリングしたデータの例を示す図である。図６の○印６０はＲＳ振幅の時系列データを示し、＊印６１は異常値処理部５による補間後にリサンプリング部８が再サンプリングしたデータを示している。
　周波数解析部９は、リサンプリング部８が取得した生体信号（ＲＳ振幅）の時系列データを高速フーリエ変換または最大エントロピー法（ＭＥＭ）により周波数解析し、生体信号のスペクトルを求める（図４ステップＳ８）。
　表示部１０は、周波数解析部９による周波数解析結果のスペクトルを表示する（図４ステップＳ９）。
　図７はＭＥＭにより周波数解析したスペクトルの例を示す図である。図７の実線７０はリサンプリング部８が取得したＲＳ振幅の時系列データ（図６の６１）をＭＥＭにより周波数解析したスペクトルを示している。ここでは、比較のため、補間処理前のＲＳ振幅の時系列データ（図６の６０）をＭＥＭにより周波数解析したスペクトルを破線７１で示している。
　補間処理前のＲＳ振幅のデータから得られたスペクトルと、不適切なデータを補間した後のＲＳ振幅のデータから得られたスペクトルとでは様相が異なっており、補間処理前のＲＳ振幅のデータから得られたスペクトルには、不適切なデータに起因するスペクトルの成分ないしは分布（図７の７２，７３，７４）が含まれていることが見て取れる。不適切なデータを補間した後のＲＳ振幅のデータから得られたスペクトルの方が、生体の状態を正しく反映した結果を示していると考えられる。
　以上のように、本実施例では、ＲＳ振幅などの生体信号に混入したノイズ成分を適切に除去することができ、生体信号の解析の精度を向上させることができる。
　なお、本実施例では、異常値判定部４が不適切と判定した生体信号のデータをもっともらしいデータに置き換えて補間しているが、これに限るものではなく、異常値処理部５は、異常値判定部４が不適切と判定した生体信号のデータを削除（データ欠損）してもよい。
［第２の実施例］
　次に、本発明の第２の実施例について説明する。図８は本実施例に係る生体信号処理装置の構成を示すブロック図である。本実施例の生体信号処理装置は、生体信号抽出部２と、平均化処理部３ａと、異常値判定部４と、異常値処理部５と、微分部６と、変化量低下判定部７と、リサンプリング部８と、周波数解析部９と、表示部１０とを備えている。本実施例では、Ｒ－Ｒ間隔の時系列データから不適切なデータを取り除く例について説明する。
　次に、図９を参照して本実施例の生体信号処理装置の動作を説明する。生体信号抽出部２は、心電計１が測定した心電図波形から生体信号（本実施例ではＲ－Ｒ間隔）を抽出する（図９ステップＳ１０）。Ｒ－Ｒ間隔を求める方法としては、例えば文献「“ECG Implementation on the TMS320C5515 DSP Medical Development Kit (MDK) with the ADS1298 ECG-FE”，Texas Instruments Incorporated，＜http://www.ti.com/lit/an/sprabj1/sprabj1.pdf＞，2011」に開示された技術を用いることができる。この文献に開示された技術では、心電図波形を時間差分した値の変化を基にＲ－Ｒ間隔を求めている。あるいは特開２０１５－１５６９３６号公報に開示された技術等を用いて心拍時刻を求め、心拍時刻の間隔をＲ－Ｒ間隔としてもよい。
　異常値判定部４は、生体信号抽出部２が抽出した生体信号（Ｒ－Ｒ間隔）のデータと平均化処理部３ａが直前の時刻までの生体信号（Ｒ－Ｒ間隔）のデータを用いて算出した平均化データとを比較し、生体信号抽出部２が抽出した生体信号のデータが適切かどうかをデータ毎に判定する（図９ステップＳ１１）。
　本実施例の異常値判定部４は、ある時刻での生体信号のデータの値Ｘ（ｉ）が、直前の時刻までの生体信号の平均化データＸ´（ｉ－１）の所定数倍（本実施例では１．３５倍）の値を超える場合、データＸ（ｉ）は不適切であると判定し、データＸ（ｉ）が平均化データＸ´（ｉ－１）の所定数倍以下の値である場合、データＸ（ｉ）は適切であると判定する。すなわち、本実施例では、平均化データＸ´（ｉ－１）の所定数倍以下の範囲が正常値の範囲となる。第１の実施例で説明したとおり、異常値判定部４が判定処理を実施するのは、生体信号抽出部２が抽出した２個目以降の生体信号のデータについてである。
　次に、平均化処理部３ａは、生体信号抽出部２が抽出した生体信号（Ｒ－Ｒ間隔）の時系列データに対して平均化処理を行う（図９ステップＳ１２，Ｓ１３）。本実施例の平均化処理については後述する。
　異常値処理部５は、異常値判定部４が不適切と判定した生体信号（Ｒ－Ｒ間隔）のデータを適切なデータに置き換えて補間する（図９ステップＳ１４）。まず、異常値処理部５は、異常値判定部４が生体信号（Ｒ－Ｒ間隔）に不適切なデータが生じていると判定した場合、不適切なデータが生じた時刻とその直前の時刻との間に挿入すべきデータの個数（心拍数に相当）を推定する。具体的には、異常値処理部５は、生体信号（Ｒ－Ｒ間隔）に不適切なデータが生じた時刻ｔ2とその直前のデータの時刻ｔ1との時間間隔（ｔ2―ｔ1）を、時刻ｔ1までのＲ－Ｒ間隔の平均化データの値で除することによって、時刻ｔ2と時刻ｔ1との間に挿入するデータの個数Ｎを決める。
　そして、異常値処理部５は、決定したデータの個数Ｎで時刻ｔ2と時刻ｔ1との時間間隔（ｔ2―ｔ1）を等分した値を、時刻ｔ2と時刻ｔ1との間に挿入するＲ－Ｒ間隔のもっともらしい値として算出する。こうして、異常値処理部５は、Ｒ－Ｒ間隔に不適切なデータが生じた時刻ｔ2とその直前の時刻ｔ1との間に、個数ＮだけＲ－Ｒ間隔のもっともらしい値を挿入することによりＲ－Ｒ間隔を補間することができる。
　図１０の（ａ）は心電図波形、および心電図波形から抽出したＲ－Ｒ間隔の時系列データの例を示す図である。図１０の（ａ）において、８０は心電図波形（単位［μＶ］）を示し、８１はＲ－Ｒ間隔（単位［ｍｓ］）を示している。この図１０の（ａ）に示した例では、１６時３１分４０秒付近の時刻で、生体の呼吸運動に由来してＲ－Ｒ間隔が大きく揺らいでいる。
　図１０の（ｂ）は図１０の（ａ）に示したＲ－Ｒ間隔の時系列データを平均化したデータ、および異常値処理部５が補間したデータの例を示す図である。図１０の（ｂ）の破線８２はＲ－Ｒ間隔の時系列データを第１の実施例で説明した手法で平均化したデータを示し、＊印８３は補間後のデータを示している。
　図１０の（ｂ）の例では、Ｒ－Ｒ間隔を平均化した値を求める際に、Ｒ－Ｒ間隔そのものの値を用い平均化している。また、誤った値が平均化に混入することを防ぐために、それまでの平均化した値に比べ一定の割合（この例では±３０％超）離れた値の場合には平均化に取り込まないこととしている。しかしそのため、本来の生体の状態に由来した揺らぎを示しているＲ－Ｒ間隔の値を不適切な値と見なして平均化に取り込まず、過去の値に固執する形となってしまい、それに基づいて無用な補間を続けるといった、誤った結果を招いていることが分かる。
　そこで、本実施例の平均化処理部３ａは、生体信号抽出部２が抽出した生体信号（Ｒ－Ｒ間隔）の平均化処理に際して、Ｒ－Ｒ間隔の逆数に基づく値に対して平均化処理を行い、この平均化処理で求めた値の逆数からＲ－Ｒ間隔の平均化データを算出することを、個々のデータ毎に行う。具体的には、Ｒ－Ｒ間隔の逆数に基づく値（以下、処理対象データＣと呼ぶ）として、Ｃ＝６００００／Ｒ－Ｒ間隔＝心拍数を用いる。ｉ番目の処理対象データをＣ（ｉ）、ｉ－１番目までの処理対象データを平均化した値をＣ´（ｉ－１）、所定の係数をｒとすれば、ｉ番目までの処理対象データを平均化した値Ｃ´（ｉ）は、式（１）と同様の次式によって求めることができる。
　Ｃ´（ｉ）＝ｒ×Ｃ（ｉ）＋（１－ｒ）×Ｃ´（ｉ－１）　・・・（２）
　平均化処理部３ａの逆数平均化処理部３０は、式（２）によりＣ´（ｉ）を算出する（図９ステップＳ１２）。そして、平均化処理部３ａの平均化データ算出部３１は、Ｒ－Ｒ間隔＝６００００／Ｃ´（ｉ）により、Ｒ－Ｒ間隔の平均化データを算出すればよい（図９ステップＳ１３）。
　ただし、平均化処理部３ａは、誤った値が平均化処理に混入することを防ぐために、ある時刻での処理対象データＣ（ｉ）が、直前の時刻までの処理対象データを平均化した値Ｃ´（ｉ－１）を中心とする所定の正常値範囲から外れている場合、処理対象データＣ（ｉ）を不適切であると判定して、平均化処理に取り込まないようにする。例えば処理対象データＣ（ｉ）が不適切と判定した場合、直前の時刻までの処理対象データを平均化した値Ｃ´（ｉ－１）をそのままＣ´（ｉ）とすればよい。ここでは、平均化データＣ´（ｉ－１）の±３０％の範囲を正常値範囲とする。
　図９のステップＳ１４，Ｓ１５，Ｓ１６，Ｓ１７，Ｓ１８，Ｓ１９の処理は、それぞれ図４のステップＳ４，Ｓ５，Ｓ６，Ｓ７，Ｓ８，Ｓ９と同様であるので、説明は省略する。
　図１１は図１０の（ａ）に示したＲ－Ｒ間隔の時系列データを本実施例の手法で平均化したデータ、および異常値処理部５が補間したデータの例を示す図である。図１１の○印および＊印９１はＲ－Ｒ間隔の時系列データおよび補間後のＲ－Ｒ間隔のデータを示し、破線９２はＲ－Ｒ間隔の平均化データを示している。
　最も変動が急な部分のＲ－Ｒ間隔の変動域がおよそ６５０ｍｓ→９００ｍｓ（変動量４２％）であるのに対し、Ｒ－Ｒ間隔の逆数に比例する心拍数の変動域はおよそ９２ｂｐｍ→６７ｂｐｍ（変動量２７％）となる。したがって、本実施例の手法によれば、Ｒ－Ｒ間隔が揺らいでいる箇所でも、それらの揺らぎは平均化に取り込まれており、補間処理したＲ－Ｒ間隔のデータが補間処理前のＲ－Ｒ間隔のデータから逸脱することがない。つまり、心拍数を用いることにより、変動のスケールが抑えられるので、平均化処理が安定化され、生体信号のデータ列に対する処理が正しく行なわれている。
　生体信号のデータの値が変動する場合、ある数値域においては、データの値そのものよりも、データの値の逆数に基づく値に対して平均化処理を行うことで、平均化データの変動幅を抑えることができ、平均化処理を安定化することができる。本実施例によれば、生体信号の時系列データにおいて、ノイズ等に起因する不適切なデータを取り除いて、もっともらしく修復することができ、それにより生体の状態のより正確な分析に結び付けることができる。
［第３の実施例］
　次に、本発明の第３の実施例について説明する。本実施例においても、生体信号処理装置の構成は第１の実施例と同様であるので、図３の符号を用いて説明する。本実施例は、図３の生体信号処理装置の微分部６と変化量低下判定部７の動作を説明するものである。
　図１２を参照して本実施例の生体信号処理装置の動作を説明する。図１２のステップＳ１～Ｓ６の処理は第１の実施例で説明したとおりである。
　次に、微分部６は、異常値処理部５によって処理された生体信号（ＲＳ振幅）の１階微分値と２階微分値を算出する（図１２ステップＳ２０）。
　ある時刻ｔkにおける生体信号（ＲＳ振幅）の補間された値をｆ（ｔk）とすると、１階微分値ｆ’（ｔk）は式（３）のように表すことができる。
　ｆ’（ｔk）＝｛ｆ（ｔk+1）－ｆ（ｔk）｝／（ｔk+1－ｔk）・・・（３）
　また、２階微分値ｆ”（ｔk）は式（４）のように表すことができる。
　ｆ”（ｔk）＝｛ｆ（ｔk+1）－２ｆ（ｔk）＋ｆ（ｔk-1）｝／（ｔk+1－ｔk）2
　　　　　　　　　　　　　　　　　　　　　　　　　　　　・・・（４）
　ＲＳ振幅の変化は呼吸運動に起因するため、息を止めているときなど、呼吸運動が停止している際は、１階微分値と２階微分値が共に０近傍の値をとる。そこで、変化量低下判定部７は、微分部６が算出した１階微分値と２階微分値が共に０を中心とする所定の範囲内の状態が所定時間以上継続した場合に、生体信号の変動が低い（呼吸運動が停止している）と判定し、１階微分値と２階微分値のうち少なくとも一方が０を中心とする所定の範囲から外れているか、あるいは１階微分値と２階微分値が共に所定の範囲内の状態の継続時間が所定時間未満の場合に、生体信号の変動が正常（呼吸運動が経時変化している）と判定する（図１２ステップＳ２１）。
　表示部１０は、変化量低下判定部７の判定結果を表示する（図１２ステップＳ２２）。これにより、生体信号の変動が低いことを警告することが可能になる。警告の方法としては、通知手段となる表示部１０に警告メッセージを表示させてもよいし、ＬＥＤ等の発光器からなる通知手段を別に設けて点滅させてもよいし、スピーカ等の通知手段を別に設けて音声で警告を発するようにしてもよい。
　図１３は、図５に示したＲＳ振幅の時系列データの１階微分値および２階微分値と、変化量低下判定部７の判定結果の例を示す図である。図１３の□印１１０は１階微分値を示し、△印１１１は２階微分値を示し、実線１１２は変化量低下判定部７の判定結果を示している。
　この図１３の例では、１階微分値と２階微分値が共に±８０［μＶ］の範囲内の状態が３秒以上継続した場合に、生体信号の変動が低いとして判定結果を“１”とし、１階微分値と２階微分値のうち少なくとも一方が±８０［μＶ］の範囲から外れているか、あるいは１階微分値と２階微分値が共に±８０［μＶ］の範囲内の状態の継続時間が３秒未満の場合に、生体信号の変動が正常として判定結果を“０”としている。
　こうして、本実施例では、呼吸運動が経時変化しているのか、停止しているのかを判定することができ、生体の呼吸運動を監視することができる。
　なお、本実施例では、図３に示した生体信号処理装置の微分部６と変化量低下判定部７の動作について説明しているが、図８に示した第２の実施例においても同様の動作を実現することができる。
［第４の実施例］
　次に、本発明の第４の実施例について説明する。本実施例においても、生体信号処理装置の構成および処理の流れは第１の実施例と同様であるので、図３、図４の符号を用いて説明する。本実施例では、ＲＳ振幅の時系列データを周波数解析する例について説明する。
　第１の実施例で説明したとおり、生体信号処理装置の周波数解析部９は、リサンプリング部８が取得した生体信号（ＲＳ振幅）の時系列データを周波数解析して生体信号のスペクトルを求め（図４ステップＳ８）、表示部１０は、このスペクトルを表示する（図４ステップＳ９）。
　図１４はＲＳ振幅を高速フーリエ変換により周波数解析したスペクトルの例を示す図である。横軸は周波数、縦軸は最大ピーク値が１となるように規格化されたスペクトルエネルギー密度である。このように、周波数分布を得たことにより、計測されたＲＳ振幅は、どの周波数帯域の成分により構成されているのか把握可能となり、呼吸運動のより特徴的な知見を得ることができる。
　図１４の例では、０．２秒おきにサンプリングされたＲＳ振幅の時系列データを３２個すなわち、計６．２秒間のデータからフーリエ変換の計算を行っている。高速フーリエ変換における周波数分解能は、次の式（５）で規定されるため、横軸のプロット間隔はおよそ０．１６Ｈｚである（文献「三上直樹，“はじめて学ぶディジタル・フィルタと高速フーリエ変換”，ＣＱ出版，ｐ．１３５－１３７，２００５年」参照）。
　ｄｆ＝１／Ｔ（＝１／６．４＝１／（０．２×３２）≒０．１６）　　  ・・・（５）
　式（５）が意味するところは、周波数分解能は計測時間によって一意に定まるということである。周波数分解能を０．１６Ｈｚよりも高くしたい場合には、リサンプリング部８のサンプリングレートが固定されている状況においては、計測時間を長くするより他ない。
　図１５は、図１４の場合と同じＲＳ振幅のデータを用いて最大エントロピー法により周波数解析したスペクトルの例を示す図である。最大エントロピー法は高速フーリエ変換とは異なり、周波数分解能は計測時間に無関係であるため、計測時間を長くせずとも高い周波数分解能が得られる。
　平時の呼吸の周期は３～４秒であり、高速フーリエ変換で周波数分解能を高めようとして計測時間を長くすると複数周期分のデータを含むこととなり、得られる周波数分布はそれらの統計量となってしまい、単回の呼吸活動の情報が埋もれてしまう。しかし、最大エントロピー法を用いれば高い分解能を短い計測時間で得られるため、統計量ではない単回の呼吸についての周波数情報を得ることができる。
　図１５においては、図１４に比べて周波数分解能が高まったために、呼吸周波数のピークが０．３５Ｈｚにあることが明瞭に観測でき、かつ、図１４では確認できない０．２３Ｈｚのピークも検出できている。
　最大エントロピー法による周波数の計算は、以下の方法により行う（文献「南茂夫，“科学計測のための波形データ処理”，ＣＱ出版，ｐ．１７３－１７４，１９８６年」参照）。最大エントロピー法には、Ｂｕｒｇ法とＹｕｌｅ－Ｗａｌｋｅｒ方があるが、ここではＢｕｒｇ法を例にとる。
　求めたいスペクトルエネルギー密度Ｓ（ω）は次の式で与えられる。
Figure JPOXMLDOC01-appb-M000001
　Δｔはサンプリングレートであり、図１５の例では０．２秒である。ωは角周波数であり、周波数ｆとω＝２πｆの関係をもつ。
　スペクトルエネルギー密度Ｓ（ω）を求めるためには、自己回帰モデルの係数ａmi、予測誤差の分散Ｐm、および用いるモデル次数ｍが分かればよい。モデル次数ｍは任意の整数をとるが、ここでは最大モデル次数を１６以上の値（文献「井上博，“循環器疾患と自律神経機能”，医学書院，ｐ．８５－８６，２０１０年」参照）から２０を選択し、次数ｍは１から２０をとるものとする。
　自己回帰モデルの係数ａmiを求めるためには、まず式（７）～式（９）からａmmを求めておく必要がある。
Figure JPOXMLDOC01-appb-M000002
　式（７）におけるＮはＲＳ振幅のデータの個数であり、ここでは３２である。また、式（７）～式（９）中の係数ｂmi，ｂ’miの初期値は次式のようになる。ｘiはＮ個のデータのうちのｉ番目のデータを現す。
Figure JPOXMLDOC01-appb-M000003
　求めたａmmから式（１３）、式（１４）の漸化関係式を用いて、自己回帰モデルの係数ａmiと予測誤差の分散Ｐmを求める。
Figure JPOXMLDOC01-appb-M000004
　式（１４）で用いるＰ0は、ＲＳ振幅のデータの平均値をｘaveとすると式（１５）で求められる。
Figure JPOXMLDOC01-appb-M000005
　次に、式（６）に代入すべき係数ａmiを決定するための統計量Ｑmを計算する。
Figure JPOXMLDOC01-appb-M000006
　式（１６）を用いて統計量ＱmをＱ1からＱ20まで計算し、それらのうち最小のＱmを与えるｍが、式（６）で使用するモデル次数ｍである。計算はｍ＝１から１ずつ増加させてＱ20まで計算する途中において統計量Ｑmの極小値が最初に出現した段階で計算を打ち切り、Ｑmが極小値となったときのｍを式（６）のモデル次数ｍとして使用してもよい。仮にＱ20まで計算しても統計量の極小値が現れなかった場合には、次数の候補の最大値（この場合は２０）を使用する。
　以上により、自己回帰モデルの係数ａmi、予測誤差の分散Ｐmおよびモデル次数ｍが入手できたので、式（６）により、周波数分布を求めることができる。図１５では、周波数分解能をΔｆ＝０．００１とした。すなわち、角周波数ωを０から０．００６２８ずつ増加させて、３．１４まで計算した。高速フーリエ変換の際のΔｆ＝０．１６と比べて高い分解能を得ることができた。
［第５の実施例］
　第１～第４の実施例では、判定の対象となる生体信号のデータＸ（ｉ）が、このデータＸ（ｉ）よりも過去の時刻の生体信号のデータを用いて算出された平均化データＸ´（ｉ－１）に基づく所定の正常値の範囲外の場合に、判定の対象となる生体信号のデータＸ（ｉ）を不適切と判定していたが、判定処理は上記の例に限るものではない。例えば、平均化データＸ´（ｉ－１）±α（αは規定値）の範囲を正常値の範囲としてもよい。
　また、第１～第４の実施例の異常値判定部４は、判定の対象となる生体信号のデータＸ（ｉ）とこのデータよりも過去の時刻の生体信号のデータ（直前の時刻までの生体信号のデータ）とを用いて算出した平均化データから得られる分散σ2を算出してもよい。そして、異常値判定部４は、この分散σ2が、過去の時刻の生体信号のデータを用いて算出した平均化データから得られる分散σp 2に基づく所定の正常値の範囲外の場合に、生体信号のデータＸ（ｉ）を不適切と判定するようにしてもよい。例えば２σp 2以下の範囲を正常値の範囲とし、分散σ2が２σp 2を上回る場合に、生体信号のデータＸ（ｉ）を不適切と判定すればよい。
［第６の実施例］
　第１～第５の実施例で説明した生体信号処理装置は、ＣＰＵ（Central Processing Unit）、記憶装置及びインターフェースを備えたコンピュータと、これらのハードウェア資源を制御するプログラムによって実現することができる。このコンピュータの構成例を図１６に示す。コンピュータは、ＣＰＵ４０と、記憶装置４１と、インターフェース装置（以下、Ｉ／Ｆと略する）４２とを備えている。Ｉ／Ｆ４２には、心電計１などが接続される。このようなコンピュータにおいて、本発明の生体信号処理方法を実現させるためのプログラムは、フレキシブルディスク、ＣＤ－ＲＯＭ、ＤＶＤ－ＲＯＭ、メモリカードなどの記録媒体に記録された状態で提供され、記憶装置４１に格納される。ＣＰＵ４０は、記憶装置４１に記憶されたプログラムに従って第１～第５の実施例で説明した処理を実行する。
　本発明は、心電図波形から得られる生体信号を解析する技術に適用することができる。
　１…心電計、２…生体信号抽出部、３，３ａ…平均化処理部、４…異常値判定部、５…異常値処理部、６…微分部、７…変化量低下判定部、８…リサンプリング部、９…周波数解析部、１０…表示部、３０…逆数平均化処理部、３１…平均化データ算出部。

Claims (9)

  1. 　生体の心電図波形から生体信号を抽出する第１のステップと、
    　この第１のステップで抽出した生体信号の時系列データを用いて平均化データを算出する第２のステップと、
    　前記第１のステップで抽出した生体信号のデータが適切かどうかを、このデータよりも過去の時刻の生体信号のデータを用いて算出した前記平均化データに基づいてデータ毎に判定する第３のステップと、
    　この第３のステップで不適切と判定した生体信号のデータの削除および補間のいずれかを行う第４のステップとを含むことを特徴とする生体信号処理方法。
  2. 　請求項１記載の生体信号処理方法において、
    　前記第２のステップは、
    　前記第１のステップで抽出した生体信号の時系列データのそれぞれの値の逆数に基づく値に対して平均化処理を行う第５のステップと、
    　この第５のステップの平均化処理で求めた値の逆数から前記平均化データを算出する第６のステップとを含むことを特徴とする生体信号処理方法。
  3. 　請求項１または２記載の生体信号処理方法において、
    　前記第３のステップは、
    　前記判定の対象となる生体信号のデータが、このデータよりも過去の時刻の生体信号のデータを用いて算出した前記平均化データに基づく所定の正常値の範囲外の場合に、前記判定の対象となる生体信号のデータを不適切と判定するか、もしくは前記判定の対象となる生体信号のデータとこのデータよりも過去の時刻の生体信号のデータとを用いて算出した前記平均化データから得られる分散が、前記過去の時刻の生体信号のデータを用いて算出した前記平均化データから得られる分散に基づく所定の正常値の範囲外の場合に、前記判定の対象となる生体信号のデータを不適切と判定するステップを含むことを特徴とする生体信号処理方法。
  4. 　請求項１乃至３のいずれか１項に記載の生体信号処理方法において、
    　前記第２のステップは、前記第１のステップで抽出した生体信号のデータの平均化処理に際して、このデータよりも過去の時刻のデータを用いて算出した前記平均化データの値を中心とする所定の範囲から外れている生体信号のデータ、および前記第４のステップで補間した生体信号のデータを、前記平均化データの算出に使用しないことを特徴とする生体信号処理方法。
  5. 　請求項１乃至４のいずれか１項に記載の生体信号処理方法において、
    　前記生体信号は、Ｒ波と１つ前のＲ波の時間間隔であるＲ－Ｒ間隔、およびＲ波のピーク値からＳ波のピーク値までのＲＳ振幅のいずれかであることを特徴とする生体信号処理方法。
  6. 　請求項１乃至５のいずれか１項に記載の生体信号処理方法において、
    　前記第４のステップで処理した後の生体信号の１階微分値と２階微分値とを算出する第７のステップと、
    　前記１階微分値と前記２階微分値とが共に０を中心とする所定の範囲内の状態が所定時間以上継続した場合に、生体信号の変動が低いことを警告する第８のステップとをさらに含むことを特徴とする生体信号処理方法。
  7. 　請求項１乃至６のいずれか１項に記載の生体信号処理方法において、
    　前記第４のステップで処理した後の生体信号を高速フーリエ変換および最大エントロピー法のいずれかにより周波数解析して、生体信号のスペクトルを求める第９のステップをさらに含むことを特徴とする生体信号処理方法。
  8. 　生体の心電図波形から生体信号を抽出するように構成された生体信号抽出部と、
    　この生体信号抽出部で抽出された生体信号の時系列データを用いて平均化データを算出するように構成された平均化処理部と、
    　前記生体信号抽出部で抽出された生体信号のデータが適切かどうかを、このデータよりも過去の時刻の生体信号のデータを用いて算出された前記平均化データに基づいてデータ毎に判定するように構成された異常値判定部と、
    　この異常値判定部で不適切と判定された生体信号のデータの削除および補間のいずれかを行うように構成された異常値処理部とを備えることを特徴とする生体信号処理装置。
  9. 　請求項８記載の生体信号処理装置において、
    　前記平均化処理部は、
    　前記生体信号抽出部で抽出された生体信号の時系列データのそれぞれの値の逆数に基づく値に対して平均化処理を行うように構成された逆数平均化処理部と、
    　この逆数平均化処理部の平均化処理で求めた値の逆数から前記平均化データを算出するように構成された平均化データ算出部とから構成されることを特徴とする生体信号処理装置。
PCT/JP2017/002483 2016-02-04 2017-01-25 生体信号処理方法および生体信号処理装置 WO2017135116A1 (ja)

Priority Applications (3)

Application Number Priority Date Filing Date Title
CN201780009657.6A CN108601546B (zh) 2016-02-04 2017-01-25 生物信号处理方法和生物信号处理设备
JP2017565497A JP6687645B2 (ja) 2016-02-04 2017-01-25 生体信号処理方法および生体信号処理装置
US16/075,116 US10918302B2 (en) 2016-02-04 2017-01-25 Biological signal processing method and biological signal processing apparatus

Applications Claiming Priority (2)

Application Number Priority Date Filing Date Title
JP2016019500 2016-02-04
JP2016-019500 2016-02-04

Publications (1)

Publication Number Publication Date
WO2017135116A1 true WO2017135116A1 (ja) 2017-08-10

Family

ID=59499544

Family Applications (1)

Application Number Title Priority Date Filing Date
PCT/JP2017/002483 WO2017135116A1 (ja) 2016-02-04 2017-01-25 生体信号処理方法および生体信号処理装置

Country Status (4)

Country Link
US (1) US10918302B2 (ja)
JP (1) JP6687645B2 (ja)
CN (1) CN108601546B (ja)
WO (1) WO2017135116A1 (ja)

Cited By (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
WO2019198691A1 (ja) * 2018-04-11 2019-10-17 シャープ株式会社 情報処理装置、およびウェアラブル端末
CN110432885A (zh) * 2019-09-11 2019-11-12 东北大学 一种光电容积脉搏波噪声去除方法
WO2023032281A1 (ja) * 2021-08-30 2023-03-09 ソニーグループ株式会社 情報処理装置、情報処理方法、及び、プログラム

Families Citing this family (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
KR102655669B1 (ko) * 2016-07-20 2024-04-05 삼성전자주식회사 생체신호의 특징 추출 장치 및 방법과, 생체정보 검출 장치
US11123004B2 (en) * 2018-12-20 2021-09-21 Biosense Webster (Israel) Ltd. Electrophysiological ripple mapping visualization method

Citations (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP2004337408A (ja) * 2003-05-16 2004-12-02 Seiko Instruments Inc 生体信号測定装置
JP2007535392A (ja) * 2004-05-01 2007-12-06 ビーエスピー バイオロジカル シグナル プロセッシング リミテッド 高周波ｑｒｓ群の分析のための装置および方法
JP2011098214A (ja) * 2004-03-24 2011-05-19 Dainippon Sumitomo Pharma Co Ltd センサを有する生体情報計測用衣服、生体情報計測システムおよび生体情報計測装置、および装置制御方法
JP2012065713A (ja) * 2010-09-21 2012-04-05 Gifu Univ 心電図データの異常心拍及びトレンドの除去方法、自律神経モニタ装置、及び敗血症発症警告装置
JP2013078543A (ja) * 2011-10-05 2013-05-02 Win Human Recorder Co Ltd 自律神経活動指標算出方法、自律神経活動指標算出装置とそれを用いたココロとカラダのバランスの視覚的表示システム
JP2013532573A (ja) * 2010-08-05 2013-08-19 レブ−エル ダイアグノスティクス オブ ハート ディジーズ リミテッド 心拍変動を測定する方法

Family Cites Families (15)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JPS526039B2 (ja) * 1971-12-17 1977-02-18
KR20030004387A (ko) 2001-03-06 2003-01-14 마이크로스톤 가부시키가이샤 신체 동작 검출 장치
JP3946108B2 (ja) * 2002-08-27 2007-07-18 パイオニア株式会社 心拍変動解析装置、心拍変動解析方法、および心拍変動解析用プログラム
US7699785B2 (en) * 2004-03-01 2010-04-20 Sleep System Laboratory Inc. Method for determining sleep stages
WO2008157544A1 (en) * 2007-06-18 2008-12-24 New York University Electronic identity card
CN201150533Y (zh) 2007-11-30 2008-11-19 沈阳东软医疗系统有限公司 一种血压测量装置
DK2346395T3 (en) * 2008-09-22 2018-06-06 Cheetah Medical Inc SYSTEM AND PROCEDURE FOR DETERMINING BLOOD FLOW
FR2946171B1 (fr) * 2009-05-29 2011-07-15 Groupe Des Ecoles De Telecommunications Get Ecole Nationale Superieure Des Telecommunications Enst Procede de quantification de l'evolution de pathologies impliquant des changements de volumes de corps, notamment de tumeurs
US9642534B2 (en) * 2009-09-11 2017-05-09 University Of Virginia Patent Foundation Systems and methods for determining location of an access needle in a subject
CN102846321B (zh) 2011-06-27 2015-03-25 联想(北京)有限公司 一种产生报警信息的方法和移动终端
JP6063775B2 (ja) * 2013-03-01 2017-01-18 東洋紡株式会社 居眠り予防方法、及び居眠り予防装置
WO2014147939A1 (ja) 2013-03-22 2014-09-25 パナソニック株式会社 生体信号計測システム、装置、方法およびそのプログラム
JP6243254B2 (ja) 2014-02-24 2017-12-06 日本電信電話株式会社 心拍検出方法および心拍検出装置
US9827992B2 (en) * 2014-03-19 2017-11-28 Ford Global Technologies, Llc Driver anomaly detection
CN104622446A (zh) * 2015-02-03 2015-05-20 南京理工大学 一种基于khm聚类算法的心率变异信号优化方法

Patent Citations (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP2004337408A (ja) * 2003-05-16 2004-12-02 Seiko Instruments Inc 生体信号測定装置
JP2011098214A (ja) * 2004-03-24 2011-05-19 Dainippon Sumitomo Pharma Co Ltd センサを有する生体情報計測用衣服、生体情報計測システムおよび生体情報計測装置、および装置制御方法
JP2007535392A (ja) * 2004-05-01 2007-12-06 ビーエスピー バイオロジカル シグナル プロセッシング リミテッド 高周波ｑｒｓ群の分析のための装置および方法
JP2013532573A (ja) * 2010-08-05 2013-08-19 レブ−エル ダイアグノスティクス オブ ハート ディジーズ リミテッド 心拍変動を測定する方法
JP2012065713A (ja) * 2010-09-21 2012-04-05 Gifu Univ 心電図データの異常心拍及びトレンドの除去方法、自律神経モニタ装置、及び敗血症発症警告装置
JP2013078543A (ja) * 2011-10-05 2013-05-02 Win Human Recorder Co Ltd 自律神経活動指標算出方法、自律神経活動指標算出装置とそれを用いたココロとカラダのバランスの視覚的表示システム

Cited By (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
WO2019198691A1 (ja) * 2018-04-11 2019-10-17 シャープ株式会社 情報処理装置、およびウェアラブル端末
JPWO2019198691A1 (ja) * 2018-04-11 2021-03-11 シャープ株式会社 情報処理装置、およびウェアラブル端末
CN110432885A (zh) * 2019-09-11 2019-11-12 东北大学 一种光电容积脉搏波噪声去除方法
WO2023032281A1 (ja) * 2021-08-30 2023-03-09 ソニーグループ株式会社 情報処理装置、情報処理方法、及び、プログラム

Also Published As

Publication number Publication date
JP6687645B2 (ja) 2020-04-28
CN108601546B (zh) 2021-03-26
CN108601546A (zh) 2018-09-28
US10918302B2 (en) 2021-02-16
US20190038164A1 (en) 2019-02-07
JPWO2017135116A1 (ja) 2018-08-16

Similar Documents

Publication Publication Date Title
WO2017135116A1 (ja) 生体信号処理方法および生体信号処理装置
JP5303802B2 (ja) 心電図から導出された無呼吸／低呼吸指数
EP3334337B1 (en) Monitoring of sleep phenomena
EP2953527B1 (en) Respiratory rate measurement
JP5475658B2 (ja) 痛み判定装置、痛み判定プログラムおよび痛み判定装置の制御方法
EP3381364B1 (en) Respiratory estimation method and device
Mirmohamadsadeghi et al. Respiratory rate estimation from the ECG using an instantaneous frequency tracking algorithm
JP2018011819A (ja) 生体信号処理方法および装置
JP5929020B2 (ja) 意識状態推定装置及びプログラム
JP6352826B2 (ja) Ｒ−ｒ間隔補間方法および心拍変動計測装置
JP6709116B2 (ja) 呼吸検出装置、呼吸検出方法および呼吸検出用プログラム
US9706945B2 (en) Respiration rate determination in impedance pneumography
JP4528583B2 (ja) 生体負荷検査装置
JP7146803B2 (ja) Ｅｍｇ信号から吸入を検出する及び神経呼吸ドライブの量を抽出する方法並びにシステム
KR20140114181A (ko) 심전도 신호에 기반하여 스트레스를 분석하고 추정하는 방법 및 장치
JP4488016B2 (ja) 心音計
US11311243B2 (en) Blood pressure data processing apparatus, blood pressure data processing method, and blood pressure data processing program
JP3314521B2 (ja) 心拍変動波形解析方法及び装置
JP2005102781A (ja) 睡眠時無呼吸型判別方法
EP2954842B1 (en) Biological information measuring apparatus and blood pressure analyzing method
JP6568022B2 (ja) 生体信号分析方法および装置
JP2018011753A (ja) 呼吸推定方法および装置
Buś et al. ECG-Derived Respiration-Complete Implementation and its Evaluation with Use of Clinical Data
US10993676B2 (en) Signal processing method and apparatus

Legal Events

Date Code Title Description
121 Ep: the epo has been informed by wipo that ep was designated in this application

Ref document number: 17747275

Country of ref document: EP

Kind code of ref document: A1

ENP Entry into the national phase

Ref document number: 2017565497

Country of ref document: JP

Kind code of ref document: A

NENP Non-entry into the national phase

Ref country code: DE

122 Ep: pct application non-entry in european phase

Ref document number: 17747275

Country of ref document: EP

Kind code of ref document: A1