How the New Discrete Period Transform Method Can Process Physiological Signals

问题:

How does the sliding discrete period transform (DPT) algorithm handle variations in heart rate and signal quality in real-time applications?

How the New Discrete Period Transform Method Can Process Physiological Signals

答案:

The DPT algorithm is designed to adapt to variations in heart rate and signal quality in real-time applications. It employs a sliding window approach, where a new window of data is analyzed with each new sample. This allows the algorithm to track changes in the signal over time and adjust its analysis accordingly. Additionally, the DPT algorithm uses a combination of autocorrelation and ensemble averaging to enhance the signal-to-noise ratio and reduce the impact of artifacts. This helps to ensure accurate and reliable results even in challenging conditions.

Introduction

Signals of physiological origin can be corrupted by noise and motion artifacts, and often share the same pass band as the signals themselves.¹ Biological signals are quasi-stationary and have periods and amplitudes that can change over time.² Simple filtering of data is not possible with such signals. One popular way for extracting information is to use another signal, which is temporally linked to the data, acting as the timeframe for ensemble averaging. Although ensemble averaging has been effectively applied to oximetry signals using an external cardiac trigger from an ECG source,³ in many instances an ECG source may be unavailable. In this article, a successful effort was made to process signals without an ECG trigger, yet yielding similar results.

Initially, an algorithm was developed to conduct a form of autocorrelation and ensemble averaging.⁴ However, it was soon discovered that ensemble averaging in the time domain was unnecessary, since all pertinent information could be found in the period domain data itself. Heart rate and blood oxygen saturation could be calculated directly from the results generated by the sliding discrete period transform (DPT).

This work was started with a review of the discrete Fourier transform (DFT) since it can generate the frequency spectrum of a signal that can then be used to determine its period.^5,6 Another goal of the research was to sample the data with a very high degree of resolution. Obtaining high resolution with the DFT requires collecting a large number of data samples. Because biological signals are quasi-stationary, collecting a large number of samples using the DFT often results in spectral smearing.⁷ What was needed was an algorithm that had high resolution, but would only require a small number of samples compared to the DFT. Since the intention was to use the algorithm on real-time data of undetermined length, a sliding form of the transform similar to a sliding DFT was used.

Methods

Algorithm Requirements

The initial objective was to find an algorithm that could be used to determine the underlying fundamental period of the data even if it was stochastic and nonstationary in nature. The initial algorithm requirements were:

Be able to determine the fundamental period of any biomedical signal such as ECG and SpO₂.
Have a rapid enough response time to track cardiac heart rate period and amplitude changes in real time.
Recover rapidly from signal interruptions or from excessive noise or motion artifacts.
Have sufficient computational speed so as not to be the limiting factor in determining the sampling rate.
Have low to moderate storage space requirements and the ability to be used in low power and portable devices.

Algorithm Development

Starting with the DFT, the objective was to find the period, so the frequency terms in the DFT equations were replaced with the period and instead of incrementing frequency as with the DFT, the period was incremented. Whereas the DFT has the frequency increase in a linear fashion such as (1f₀, 2f₀, 3f₀, …), where f₀ is the first harmonic, the DPT has the period increase in a linear fashion in multiples of the sampling period T₀. Despite the similarities in the equations of the two algorithms, the DFT cannot generate the same results as the DPT, as they are fundamentally different algorithms. The DFT and the DPT can be compared by analyzing the equations that describe their implementation. For sampling frequency f_S, the frequency bin k of the N-point DFT corresponds to frequency f_K = k × f_S / N Hz, and Equation 1 is the expression for the spectrum of the k^th frequency bin for the sample sequence XI... XI + N – 1.

Equation 1.

Where k = 0, 1, 2, ... N – 1

The i^th sample of the DFT is calculated as shown in Equation 2.

Equation 2.

The DFT basis functions have the same ordinate value as the previous N^th ordinate value for each harmonic shown in Figure 1. This occurs because all the harmonics in the DFT are commensurate such that the higher harmonics are exact multiples of the lower harmonics.

Figure 1. Fourier transform sine basis functions where the red is the first harmonic (1 Hz), blue is the second (2 Hz), and green is the third (3 Hz).

The term N in the DPT must be modified for each period since the periods are not simply multiples of each other, yet vary by one sample period as shown in Figure 2.

Figure 2. Period transform complex sinusoidal basis functions for three adjacent sin and three adjacent cos functions. This example assumes that these functions are spaced by a sample period of 10 ms.

Both of the sliding forms DFT and the DPT require an implementation that has circular or recurrence buffers that hold a fixed number of the most recent samples. One buffer is used when the input data is real, while two buffers are employed when the input data is complex. The i^th sample of the DPT transform can be written as shown in Equation 3.

Equation 3.

Where the term RBS is the recurrence buffer size, T_L is the length of the longest period, and T_N is the period of the current basis element being processed. Doing this allows the beginning and ending ordinate values for each basis period to have the same value. The period s extends from the minimum period to the maximum period chosen to cover the periods in the sampled data. The implementation utilizes a set of basis functions that are the incremental phase angles of the complex sinusoids shown in Figure 2.

The DPT is somewhat difficult to implement because the basis functions are comprised of sets of complex functions that are mostly incommensurate and differ one from another by the sample period. An efficient DPT transform uses the basis phase angles shown in Figure 3. This is the form that was used in the implementations herein.

Figure 3. Period transform basis phase angles showing the values of the complex phase angles for increasing number of sample periods per minute. The ascending curve shows the cosine phase angles and the descending curve shows the sine phase angles.

The phasors are easily derived using Equation 4, where s is the set of periods from the minimum chosen period to the maximum chosen period in steps of the sample period.

Equation 4.

Algorithm Implementation

The sliding DPT transform is implemented as an IIR filter and has a signal flow diagram of a comb filter followed by a resonator, as does the implementation of the sliding form of the DFT. The comb filter delay of N samples causes the transient response to be N–1 samples in length. Heart rate tuned comb filters have been used by others with some success.⁸ Since the components of the DPT complex basis functions or phase angles are not all harmonically related, the end points of these functions do not always form continuous functions in the sample space as the DFT does. However, implementing the DPT as a sliding transform wraps the basis functions, allowing the component basis functions to become continuous. As the data and basis functions slide by, and the correlation is computed, basis function continuity is maintained.

In a sliding window algorithm, a window of length N is slid over a data array of indeterminate length. In the case of the DPT, two recurrence buffers are maintained because the DPT can process both real and imaginary input data. If the input has only a real component, which is often the case, only one recurrence buffer is used. However, the results can still be complex depending on the phase relationship between the input and the basis functions. The results are stored in two ensemble buffers, each with a length of the chosen maximum period.

MATLAB Proof of Concept Model

Equation 4 was implemented in a MATLAB^® script. An example plot is shown in Figure 4 using sine and cosine functions as inputs with amplitudes of ±1 and periods of 45 ms, 79 ms, and 175 ms. The MATLAB script was limited to periods between 400 ms (200 periods/min) and 2 s (40 periods/min). A fixed total of 5000 data samples were processed for this example. Since the input data are sinusoidal waveforms with unity amplitudes, the amplitudes of each period are also unity. One can easily see the great resolution that this transform offers.

Figure 4. Amplitude spectrum showing the values from three sets of input sinusoidal data that are incommensurate relative to each other.

The results for a sinusoidal cosine wave with a period of 73 periods per minute and with amplitude of 4.5 are shown in Figure 5. In this example, a recurrence buffer with a length of 1500 data points was used. Note that there are some small errors, with an amplitude error of 0.366% and a period error of 0.234%. The size of these errors is generally acceptable for biomedical applications. These errors are of no consequence with peripheral capillary oxygen saturation (SpO₂) measurements since the SpO₂ is calculated ratiometrically with the ratio of ratios of the red and infrared spectrophotometric signals.^9,10 See Equation 6 and Equation 7.

Figure 5. The sliding period transform for a cosine waveform with a period of 73 periods per minute and amplitude of 4.5. The amplitude error is less than 0.37% and the period error is less than 0.24%.

Results

Sliding Window DPT Use in Pulse Oximetry

For sliding window algorithms to operate properly when applied to pulse oximetry, two recurrence arrays are required: one for the red history and another for the infrared history. To complete the sliding transform picture, the updated contents of the recurrence buffers, which have the length of the period bin being processed, are rotated by the basis function corresponding to that period. The length of this buffer determines the overall resolution and once enough data has entered the process to fill these buffers, the transform results reach a stable limit and only change in amplitude or period as the input data changes. For the data processing reported, the recurrence buffers held the last 10 s of data.

The raw data was collected by researchers at Analog Devices and the software used to process the data was a sliding DPT in MATLAB script. Figure 6 shows the raw data taken from one subject along with band-pass filtering of 1 Hz to 4 Hz and using a flat smoothing moving average filter with a total width of 200 ms. The spectrums that are shown in Figure 7 are those after the recurrence buffers were filled and the spectrums reached stable amplitudes. The DPT will continue to track any changes in the raw data as new data is sampled and the spectrums will be updated accordingly.

Figure 6. Raw, filtered, and smoothed photoplethysmographic pulse data from a subject using a MAX30101 PPG AFE device. The top waveforms show the raw infrared and red signals while the bottom waveforms show the filtered and smoothed data.

Figure 7. This graph shows the red and infrared spectrums using the sliding window DPT. The larger of the two peaks is the infrared spectrum and the smaller is the red spectrum.

To estimate SpO₂, the commonly known equation referred to as the ratio of ratios was first used. The AC component entries used the peak values from the spectrum plot shown in Figure 7 and the DC component entries used the average values of the unfiltered signals shown in Figure 6.

Equation 5.

Equation 6.

Equation 7.

The SpO₂ and heart rate collected from a Masimo oximeter using its SET algorithm were compared with the data taken at the same time with an ADI MAX30101 pulse oximeter sensor. The data from a subject was chosen at random and the results plotted in Figure 8 and Figure 9.

Figure 8. Comparison of DPT processed photoplethysmographic data.

Figure 9. Comparison of processed heart rate data from the MAX30101 oximeter processed with the discrete period transform and a Masimo oximeter.

It is common medical practice to evaluate values produced by two different instruments that are measuring the same parameter. One of the instruments is assumed to produce correct results and used as the standard.

Bland and Altman developed a method based on the agreement between two quantitative measurements.^11,12 They did this by studying the mean difference and constructing the limits of agreement. The Bland-Altman plot analysis is a simple way to evaluate a bias between the mean differences and to estimate an agreement interval. If two medical instruments are put to this test, with one considered the standard, the results should remain within two standard deviations or 95% for the second device to be considered as good as the standard for clinical use.

Using the Bland-Altman method, as opposed to correlation analyses that identifies the relationship between two variables, a statistical method is applied that addresses the difference between two variables.

The accuracy and precision of the DPT algorithm was evaluated by collecting data from 26 healthy adult subjects using the MAX30101 pulse oximeter sensor and comparing that with a Masimo oximeter containing its latest Signal Extraction Technology®.¹³ The cohort consisted of 15 male and 11 female subjects with ages between 20 and 40 years. Our study was designed to compare the results of the two oximeters on each of the individuals in the cohort and not any differences between men and women. Note that SpO₂ between the sexes does vary somewhat. One study has shown that, for young healthy adult individuals, the mean SpO₂ for males was 97.1±1.2%, while for females the mean was 98.6±1.0%.¹⁴

The results using the Bland-Altman criteria are depicted in Figure 10 and Figure 11, where each circle represents the Bland-Altman results for an individual subject. All of the SpO₂ comparisons met the Bland-Altman criteria.

Figure 10. The percentage differences of SpO₂ between a Masimo oximeter and an ADI oximeter using the DPT algorithm. The Bland-Altman criteria was met.

Figure 11. Heart rate differences in beats per minute between a Masimo oximeter and an ADI oximeter using the DPT algorithm. The Bland-Altman algorithm was met in all cases but one. An arrow shows where the analysis falls outside of two standard deviations.

In Figure 11, an arrow points to a value in the comparison of heart rate that falls outside of two standard deviations. The heart rate vs. time plot for this subject is shown in Figure 12, where the Masimo oximeter had a standard deviation of 1.7892 and the MAX30101 oximeter using the DPT algorithm had a standard deviation of 0.8935. In this case, it is difficult to identify which instrument was more accurate, but the standard deviations give an indication.

Figure 12. Plot of heart rate vs. time for a Masimo oximeter and an ADI oximeter. For the 25 s period, the Masimo oximeter had a standard deviation of 1.7892 and the MAX30101 had a standard deviation of 0.8935. The stepped waveform is the signal from the Masimo oximeter and the smooth signal is from the ADI oximeter running the DPT algorithm.

A Protype Oximeter System Using the SDPT Algorithm

Finally, a prototype oximeter was designed with an Arm® microprocessor running with a bare metal operating system. The Raspberry Pi Zero was used as the computer platform with the MAX30102 integrated circuit as the sensor. The operating system and the sliding window DPT were implemented in the standard C language. Figure 13 shows the prototype. The entire oximeter was powered by a USB 3.0 connection. Two digital-toanalog converters sent data, as determined by the attending software, via a ribbon cable to a Tektronix DPO-4034 oscilloscope where they were plotted. The plots were then sent to a desktop computer using a network connection. Figure 15 shows the results from a single subject taken over a period of approximately 9 s, post a 10 s period to fill the recurrence buffer.

Figure 13. The prototype Raspberry Pi Zero-based pulse oximeter. The MAX30102 SpO₂ sensor is in the finger clip shown in the upper left corner of the picture.

The red and infrared DC signals were extracted from the raw signals using a first-order low-pass IIR filter, while the AC signals were extracted using a first-order high-pass IIR filter. See Figure 14. The time constant for these filters was set at approximately 1 s. The data was sampled at 100 SPS with the interrupt from the MAX30102 as the timing signal. The device’s output was in a 12-bit fixed point digital format for both the red and the infrared signals.

Figure 14. The infinite impulse response (IIR) filters used to extract the AC signals and the DC signals from the raw spectrographic data.

Figure 15. PPG waveforms from the Raspberry Pi Zero oximeter prototype showing the red pulses on the top and infrared pulses on the bottom. The heart rate is approximately 58 bpm. Inverted waveforms are shown to more accurately represent the actual arterial pressure in the finger.

Once the red and infrared AC signals were extracted by the filters, they were processed by the DPT without any further signal preprocessing. The first harmonics of the spectrographic signals produce peaks as shown in Figure 16. The location of the data peaks on the abscissa determines the heart rate while the amplitudes of the red and infrared data peaks were used to directly calculate the SpO₂ using the ratio of ratios equations.

Figure 16. The spectrums generated by Raspberry Pi Zero using the sliding window discrete period transform with an SpO₂ value of 97% and a heart rate of 58 bpm. The cursor b (center vertical blue line) shows where the heart period is measured as 1.03 s. The rectangular signal at the top left points to where the 400 ms period is located on the abscissa and the rectangular signal at the top right points to where 2000 ms period is located on the abscissa.

Discussion

The raw optical signals generated by oximeters contain large steady DC components and small oscillating AC components that are approximately 1% of the DC signal. These oscillating components contain the pulsatile activity in the capillaries. Any motion or other artifacts can easily override these signals and prevent an accurate reading. Over the years, a great deal of time has been spent on methods to separate the signals from any artifacts. Often the methods have proven to be very complex and difficult to implement.^16,17

It was for these reasons that this research was undertaken. The DPT algorithm resolves many of these challenges by employing a transform that requires a small number of samples yet achieves accurate measurements. Making the measurements in the period domain and spacing each bin by the sample period provides the required resolution. The period and amplitude information from the DPT could then be used to directly calculate heart rate and blood oxygen saturation without returning to the time domain.

Conclusion

Period domain analysis utilizing an incremental DPT algorithm is an effective and efficient way to process periodic biomedical signals for spectral content. It provides the capabilities of frequency domain analysis, with advantages in implementation. ADI’s MAX30101 integrated circuit sensor running the DPT algorithm was shown to be accurate enough to replace a Masimo oximeter in the practice of medicine.

Acknowledgements

We would like to especially thank Dr. Amit Nimunkar of the University of Wisconsin Department of Biomedical Engineering for his constructive suggestions and help proofreading this article. We would also like to thank Dr. Everett Smith for his valuable help and encouragement and Analog Devices for providing the raw subject data used in this project.

References

¹ Amal Jubran. “Pulse Oximetry.” Critical Care, Vol. 3, No. 2, February 1999.
² Han-Wook Lee, Ju-Won Lee, Won-Geun Jun, and Gun-Ki Lee. “The Periodic Moving Average Filter for Removing Motion Artifacts from PPG Signals.” International Journal of Control, Automation, and Systems, Vol. 5, No. 6 December 2007.
³Brendan Conlon, James A. Devine, and James A. Dittmar. “ECG Synchronized Pulse Oximeter.” U.S. Patent 4,960,126, October 1990.
⁴ James Reuss and Dennis Bahr. “Period Domain Analysis in Fetal Pulse Oximetry.” Proceedings of the Second Joint 24th Annual Conference and the Annual Fall Meeting of the Biomedical Engineering Society, 2002.
⁵ Eric Jacobsen and Richard Lyons. “The Sliding DFT.” IEEE Signal Processing Magazine, Vol. 20, No. 2, March 2003.
⁶ Eric Jacobsen and Richard Lyons. “An Update to the Sliding DFT.” IEEE Signal Processing Magazine, Vol. 21, No. 1, January 2004.
⁷ Lawrence R. Rabiner and Bernard Gold. Theory and Application of Digital Signal Processing. Prentice-Hall, January 1975.
⁸ Ludvik Alkhoury, Ji-Won Choi, Chizhong Wang, Arjun Rajasekar, Sayandeep Acharya, Sean Mahoney, Barry S. Shender, Leonid Hrebien, and Mose Kam. “Heart-Rate Tuned Comb Filters For Processing Photoplethysmogram (PPG) Signals in Pulse Oximetry.” Journal of Clinical Monitoring and Computing, Vol. 35, No. 4, August 2021.
⁹“Recommended Configurations and Operating Profiles for MAX30101/MAX30102 EV Kits.” Maxim Integrated, March 2018.
¹⁰ Sang-Soo Oak and Praveen Aroul. “How to Design Peripheral Oxygen Saturation (SpO2) and Optical Heart Rate Monitoring (OHRM) Systems Using the AFE4403.” Texas Instruments, March 2015.
¹¹ Douglas Altman and J. Martin Bland. “Measurement in Medicine: The Analysis of Method Comparison Studies.” Journal of the Royal Statistical Society: Series D (The Statistician), Vol. 32, No. 3, September 1983.
¹²J. Martin Bland and Douglas G. Altman. “Statistical Methods for Assessing Agreement Between Two Methods of Clinical Measurement.” The Lancet, February 1986.
¹³ Julian M. Goldman, Michael T. Petterson, Robert J. Kopotic, and Steven J. Barker. “Masimo Signal Extraction Pulse Oximetry.” Journal of Clinical Monitoring and Computing, Vol. 16, No. 7, 2000.
¹⁴ Sagi Levental, Elie Picard, Francis Mimouni, Leon Joseph, Tal Y Samuel, Reuben Bromiker, Dror Mandel, Nissim Arish, and Shmuel Goldberg. “Sex-Linked Difference in Blood Oxygen Saturation.” The Clinical Respiratory Journal, Vol. 12, No. 5, May 2018.
¹⁵ Sam Koblenski. “Everyday DSP for Programmers: DC and Impulsive Noise Removal.” November 2015.
¹⁶ Surekha Palreddy. “Signal Processing Algorithms.” Design of Pulse Oximeters, first edition, October 1997.
¹⁷ Terry L. Rusch, Ravi Sankar and John E. Scharf. “Signal Processing Methods for Pulse Oximetry.” Computers in Biology and Medicine, Vol. 26, No. 2, March 1996.