LiU-ITN-TEK-A--17/057--SE
Implementation of vital sign
detection algorithms on a
high-performance digital
signal processor
Tobias Pettersson
LiU-ITN-TEK-A--17/057--SE
Implementation of vital sign
detection algorithms on a
high-performance digital
signal processor
Examensarbete utfört i Elektroteknik
vid Tekniska högskolan vid
Linköpings universitet
Tobias Pettersson
Handledare Adriana Serban
Examinator Qin-Zhong Ye
Upphovsrätt
Detta dokument hålls tillgängligt på Internet – eller dess framtida ersättare –
under en längre tid från publiceringsdatum under förutsättning att inga
extra-ordinära omständigheter uppstår.
Tillgång till dokumentet innebär tillstånd för var och en att läsa, ladda ner,
skriva ut enstaka kopior för enskilt bruk och att använda det oförändrat för
ickekommersiell forskning och för undervisning. Överföring av upphovsrätten
vid en senare tidpunkt kan inte upphäva detta tillstånd. All annan användning av
dokumentet kräver upphovsmannens medgivande. För att garantera äktheten,
säkerheten och tillgängligheten finns det lösningar av teknisk och administrativ
art.
Upphovsmannens ideella rätt innefattar rätt att bli nämnd som upphovsman i
den omfattning som god sed kräver vid användning av dokumentet på ovan
beskrivna sätt samt skydd mot att dokumentet ändras eller presenteras i sådan
form eller i sådant sammanhang som är kränkande för upphovsmannens litterära
eller konstnärliga anseende eller egenart.
För ytterligare information om Linköping University Electronic Press se
förlagets hemsida
http://www.ep.liu.se/Copyright
The publishers will keep this document online on the Internet - or its possible
replacement - for a considerable time from the date of publication barring
exceptional circumstances.
The online availability of the document implies a permanent permission for
anyone to read, to download, to print out single copies for your own use and to
use it unchanged for any non-commercial research and educational purpose.
Subsequent transfers of copyright cannot revoke this permission. All other uses
of the document are conditional on the consent of the copyright owner. The
publisher has taken technical and administrative measures to assure authenticity,
security and accessibility.
According to intellectual property law the author has the right to be
mentioned when his/her work is accessed as described above and to be protected
against infringement.
For additional information about the Linköping University Electronic Press
and its procedures for publication and for assurance of document integrity,
please refer to its WWW home page:
http://www.ep.liu.se/Abstract
This thesis investigates the possibility of detecting weak vital signs, such as heart-beat and respiration rate, through the implementation of quadrature demodula-tion and frequency spectrum algorithms on a high performance digital signal pro-cessor. This thesis has been part of an ongoing research project at Link ¨oping Uni-versity, with the aim to develop a sensor platform for wireless measurements of these vital signs. This sensor platform has been expected to consist of two major physical devices, a RF-radar front-end including a quadrature multi-port and a processing back-end which holds the implementation of the algorithms to detect the vital signs. The back-end consists of a data acquisition- and a processing-part which together forms the digital signal processor.
The results show that the implemented algorithms works in terms of being able to find artificial vital signs from quadrature signals. This result also confirms that the hardware solution proposed during this thesis, has been considered as viable for the aim of the project.
Acknowledgements
Thanks to Qin-Zhong Ye, Adriana Serban and Gustav Knutsson at the depart-ment of science and technology at Link ¨oping University for equipdepart-ment, guidance and help.
Special thanks to the student Henrik Kalv´er at Link ¨oping University for pro-viding motivational speeches and occasionally being a much needed sounding board that helped provide new angles towards solving problems.
Acronyms
Abbr. Description
ADC Analog to Digital Converter DMA Direct Memory Access DSP Digital Signal Processor FFT Fast Fourier Transform FIR Finite Impulse Response
FPGA Field Programmable Gate Array GUI Graphical User Interface
IDE Integrated development environment IIR Infinite Impulse Response
MATLAB Matrix Laboratory MCU Micro Controller Unit PCB Printed Circuit Board RF Radio Frequency
RMII Reduced Media-Independent Interface USB Universal Serial Bus
Contents
1 Introduction 1 1.1 Background . . . 1 1.2 Purpose . . . 2 1.2.1 Problem statement . . . 2 1.2.2 Delimitation . . . 2 1.3 Outline . . . 3 1.4 Method . . . 4 1.4.1 Hardware . . . 4 1.4.2 Software . . . 4 2 System overview 5 2.1 Concept . . . 52.1.1 The multi-port correlator . . . 6
2.1.2 Demodulation . . . 7
2.1.3 The sensor platform . . . 8
2.2 Task . . . 9
2.3 Hardware . . . 10
2.3.1 The Nucleo board . . . 10
2.3.2 The ADCs . . . 10
2.3.3 The Ethernet controller . . . 11
3 Data acquisition 12 3.1 Sampling . . . 12
3.2 Implementation . . . 13
3.2.1 Defining the sampling clock . . . 13
3.2.2 Dual-mode . . . 15
3.2.3 Configure the DMA . . . 16
3.2.4 Verification . . . 16
4 Processing 18 4.1 Concept . . . 18
4.2 The frequency spectrum . . . 18
4.2.2 Peak detection . . . 21
4.3 Demodulation with arc-tangent . . . 22
4.3.1 Adjusting the error . . . 23
4.3.2 Rotating into first octant . . . 25
4.3.3 The Taylor approximation in code . . . 26
4.3.4 The frequency spectrum with a phase-shifted signal . . . 28
5 Pre-processing 30 5.1 Concept . . . 30
5.2 Filtering . . . 30
5.2.1 Generating coefficients . . . 31
5.2.2 Implementing the filter . . . 32
5.3 Decimation . . . 35
5.3.1 Implementation . . . 35
5.4 Sliding window effect . . . 36
5.5 Data buffers . . . 36
5.5.1 The ADC-buffer . . . 36
5.5.2 The ring buffer . . . 37
5.5.3 Other buffers . . . 38 6 Communication link 39 6.1 Concept . . . 39 6.2 UART-communication . . . 39 6.2.1 Communication protocol . . . 40 6.3 Implementation . . . 43 6.3.1 Types . . . 44
7 Graphical User Interface 45 7.1 Concept . . . 45
7.2 MATLAB . . . 45
7.2.1 Describing the GUI . . . 45
7.2.2 Implementation . . . 46
7.2.3 Finalized GUI . . . 49
8 Test benches 50 8.1 Concept . . . 50
8.2 Signals generated in MATLAB . . . 50
8.2.1 Summary . . . 52
8.3 Time estimation . . . 53
9 Results 54 9.1 Concept and result outline . . . 54
9.2 Testing the algorithm . . . 54
9.2.1 ADC measures . . . 64
9.3 Measuring the algorithm . . . 66
9.3.1 The processing loop . . . 66
9.4 Summary . . . 67
10 Discussion 69 10.1 Tests of the algorithm . . . 69
10.1.1 Addition of noise . . . 69
10.1.2 The use of decimation . . . 70
10.1.3 Artificial vital signs . . . 70
10.1.4 Usage of the ADCs . . . 70
10.1.5 The GUI . . . 71
10.2 Time consumption of the algorithm . . . 71
10.2.1 The processing . . . 71
10.2.2 The communication link . . . 71
10.2.3 The data collection . . . 72
11 Conclusion 73 11.1 In general . . . 73
11.1.1 Hardware . . . 73
11.1.2 The algorithm . . . 73
11.1.3 Test and verification . . . 74
11.2 Further development . . . 74
List of Figures
2.1 Translation of a quadrature mixer into a multi-port [4]. . . 6
2.2 Blockdiagram describing the entire sensor platform. . . 8
2.3 Flowchart of the expected data through the DSP-kernel. . . 9
3.1 Sample code showing the configuration of the timer. . . 15
4.1 Sample code showing the FFT-implementation. . . 19
4.2 Reference spectrum calculated in MATLAB of the super-position signal. . . 20
4.3 Spectrum calculated using the DSP-FFT functions from ARM. . . . 21
4.4 Sample code showing the calculation of peaks in the output spec-trum. . . 21
4.5 Comparing the MATLAB atan2 to three different orders of Taylor expansions for arctan. . . 22
4.6 Showing the error at t=1 compared to the atan2. . . . 23
4.7 Comparing the adjusted Taylor expansion to the old one and to atan2. . . . 24
4.8 The implemented code for arctan approximation in the first octant. 24 4.9 The unit circle divided into four quadrants. . . 25
4.10 The implemented code to rotate any coordinate into the first octant. 26 4.11 The implemented code for approximation arctan with Taylor ex-pansion. . . 27
4.12 The reference frequency spectrum for the phase-shifted signal. . . 28
4.13 The frequency spectrum for the phase-shifted signal, calculated with the Taylor approximation in MATLAB. . . 29
4.14 The frequency spectrum for the phase-shifted signal, calculated on the DSP. . . 29
5.1 Frequency response of the implemented Butterworth filter. . . 32
5.2 The implemented filter in C-code. . . 32
5.3 The super position signal in the time domain. . . 33
5.4 The filtered result of the super position signal, in the time domain. 33 5.5 The filtered result from the DSP, in the time domain. . . 34
5.7 Sample code for the decimation of data. . . 35
5.8 Illustration of the data flow into the ADC-buffer. . . 37
5.9 A linear buffer in a ring buffer configuration. . . 38
6.1 Illustration of the packet structure sent between devices. . . 41
6.2 Skeleton for the uartSend-function. . . . 43
6.3 The if-statement to determine the amount of packets to send. . . . 44
7.1 The first computerized sketch of the GUI. . . 46
7.2 Setting up the figure. . . 46
7.3 The figure window, generated by the sample code. . . 47
7.4 Sample code for initializing the serial connect button in MATLAB. 47 7.5 Sample code for the callback connected to the serial connect but-ton. . . 48
7.6 The final GUI-window with all the elements added. . . 49
8.1 A simple illustration of a 32-bit word and its sub-components. . . 51
8.2 An illustration of how to measure the on time of the signal. . . 53
9.1 The algorithm response from test signal 1, without arc-tangent de-modulation. . . 55
9.2 The algorithm response from test signal 1, with arc-tangent enable. 55 9.3 The algorithm response from test signal 2, without any arc-tangent demodulation. . . 56
9.4 The algorithm response from test signal 2, with arc-tangent enable. 56 9.5 The algorithm response from test signal 3. The arc-tangent module is disabled. . . 57
9.6 The algorithm response from test signal 3, with arc-tangent enable. 57 9.7 The algorithm response from test signal 4, without arc-tangent. . . 58
9.8 The algorithm response from test signal 4, with arc-tangent enable. 58 9.9 The algorithm response from test signal 5, without arc-tangent. . . 59
9.10 The algorithm response from test signal 5, with arc-tangent enable. 59 9.11 The algorithm response from test signal 6, with the arc-tangent de-modulation disabled. . . 60
9.12 The algorithm response from test signal 6, with arc-tangent enable. 61 9.13 The algorithm response from test signal 8. Arc-tangent is disabled. 62 9.14 The algorithm response from test signal 8, with arc-tangent enable. 62 9.15 The algorithm response from test signal 10, without arc-tangent demodulation. . . 63
9.16 The algorithm response from test signal 10, with arc-tangent en-able. . . 63
9.17 The algorithm response to a 10 Hz sine wave, with arc-tangent dis-abled. . . 64
9.18 The algorithm response to a 1 Hz sine wave, with arc-tangent dis-abled. . . 65
9.19 The algorithm response to a 100 mHz sine wave, with arc-tangent disable. . . 65
List of Tables
3.1 The sampling characteristics for the ADCs . . . 14
3.2 ADC sampling and trigger frequencies . . . 17
4.1 Considered requirements for the FFT-implementation . . . 18
4.2 Modified constants for arctan approximation . . . 23
4.3 Rotation of coordinates in the unit circle . . . 25
5.1 Butterworth coefficients for the filter . . . 31
5.2 Butterworth coefficients sorted into a linear buffer . . . 32
6.1 Requirements on the communication protocol . . . 41
6.2 Types of messages sent through the UART-communication . . . 44
8.1 Summary of the different configurations for test signals . . . 52
9.1 Time consumed by the three main modules . . . 66
9.2 Time consumed by the two different fputc-functions . . . 66
9.3 Time needed per byte sent through the UART . . . 67
9.4 Time consumed by the entire processing loop . . . 67
9.5 Time used to copy data into the processing loop . . . 67
9.6 Total time needed to process new data . . . 68
1
Introduction
1.1
Background
Multi-port transmitters and receivers have been investigated at Link ¨oping Uni-versity, by the Communication Electronics research group, to be used in systems for high data rate transmission. Due to the nature of the correlator component, the multi-port technique has excellent phase-discriminating capabilities, which makes them appropriate for use in radar devices in special industrial, automo-tive and medical applications.
Among this variety of applications, the wireless detection of weak, low fre-quency vital signs such as heartbeat and respiration rate are of particular interest. Using the multi-port technique for such applications, benefits from low-power consumption, good linearity and low circuit complexity. However, as with any radar device or application, the Radio Frequency (RF) front-end should be com-plemented by a Digital Signal Processor (DSP) back-end in order to present the data to the user.
Currently there is a research project at Link ¨opings University, with the aim to design and optimize a sensor platform, to measure these vital signs. The sensor platform can be divided into three parts/subsystems with the following func-tions:
• Multi-port radio/radar system
• Acquisition board consisting of A/D conversion in real time • Digital signal processor for processing and windowing data
Furthermore this sensor platform should have a communication link towards a computer or server, in order to present the data to the user. This communication link should be included as a module in the DSP-subsystem.
1.2
Purpose
This master thesis aims at implementing the acquisition- and DSP-part of the described sensor platform. How these two parts are implemented and whether they are put together in one unit or not, is left for the thesis to answer. The two parts should be fast enough to measure and process the output ports from multi-port radar, while also being able to transfer this data to another location/device. If successful, the outcome of this thesis could be used to further aid the ongoing research project.
To also enhance the possible outcome of this thesis, a Graphical User Inter-face (GUI) or a similar interInter-face has been requested. This interInter-face should allow the data collected to be easily evaluated, using for instance a diagram or a plot. This interface should however be separate from the DSP-part but continuously receive new data, to keep itself updated. A suggestion has been made towards implementing this GUI in Matrix Laboratory (MATLAB).
1.2.1
Problem statement
This thesis will be built upon answering these questions, where the first one con-cerns what type of hardware to use. The last question is considered as secondary and will only be evaluated if time is still available.
• How should the data acquisition and the DSP-part be implemented? As separate units put together into a standalone device or as an already existing solution?
• How should the DSP-kernel with its different modules be implemented in the available hardware?
• What kind of communication link between the sensor platform and the computer should be used?
• (Secondary objective) How should the computer or server be utilized in order to present the data to the user? Should MATLAB as suggested be used, or is there another way such as using a C/C++ implementation?
1.2.2
Delimitation
The sensor platform should also include the RF-part with the multi-port corre-lator and the radio/radar system to become the intended platform. This part will not be covered in this thesis but is instead part of another students thesis. However, since this thesis and a big part of the acquisition system is heavily de-pendent on the levels of signals generated by the multi-port correlator, there will be a close development together with this other student.
1.3
Outline
This report will describe the work and procedures done by the author during this master thesis. This includes all information relevant to the results and to the conclusion of the thesis.
The Introduction-chapter will provide the reader with a basic understanding of the task upon which this thesis has been conducted, in terms of purpose and problem formulations. This chapter will also describe what this thesis contains in terms of outline.
The System overview-chapter will describe the approach made towards solving the task at hand. It will also provide brief theory for the entire sensor platform along with necessary theory for the data acquisition- and DSP-parts. This chapter will also provide information about the hardware used to hold the implementa-tion of the data acquisiimplementa-tion- and DSP-parts.
The Data acquisition-chapter describes how the collection of data has been im-plemented, together with tests and evaluations of how the implementation was working. From a development point of view, the Data acquisition-chapter was first in a chronological order, and also included the first tests made to verify whether the chosen hardware was enough in terms of collecting data or not.
The Processing-chapter contains the implementation of the two major algo-rithms used in this thesis, the phase-demodulation and the frequency spectrum calculation. This chapter was also the second in chronological order, and tested whether the performance of the hardware was enough.
The Pre-processing-chapter describes the implementation of a digital filter in order remove undesired components. This chapter also describes the decimation that followed the filtering in order to reduce the sampling speed of the system, as well as the implementation of the storage of data along the execution loop.
The Communication link-chapter will provide the reader with insight on how the communication between the DSP-part and a computer is implemented. This includes the design of a protocol for transmitting larger chunks of data, as well as the re-targeting of some standard C-library functions.
The Graphical User Interface-chapter describes the implementation and design of a GUI in MATLAB, where data will be presented towards the user in a MAT-LAB plot.
The Test benches-chapter contains the description of test made in MATLAB, to verify the implemented algorithms. The chapter also describes performance measures that should be done onto the algorithm in order to determine the time consumption.
The Result-chapter contains the results collected from the tests made according to the Test benches-chapter. These results are mainly presented as figures of the implemented GUI with the corresponding results.
The Discussion-chapter does just that. It discusses the results presented, but also the overall performance and drawbacks of the implemented system.
Finally the Conclusion-chapter summarizes this master thesis by answering the problem statement from above, along with some suggestions for further de-velopment.
1.4
Method
During this thesis several different sources of information have been used in or-der to gather the necessary knowledge needed to implement the task at hand. Where needed, these external sources for information will be cited to give the reader the possibility to verify decisions and approaches made. Alongside this information, decision regarding what kind of hardware and software to use has been made. Some of these decisions have been influenced by the needs of the research project.
1.4.1
Hardware
The hardware used in this project consists of an evaluation platform and test equipment. The evaluation platform is the STM32 Nucleo board from ST [1], upon which the developed algorithms has been deployed on, while the test equip-ment is separated into two devices. The first device is an 33120A Abritrary Wave-form Generator from Hewlett Packard (Agilent) [2] and the second device is an DSO-X 2012A Oscilloscope from Agilent Technologies [3].
1.4.2
Software
The software used in this project consists of a programming Integrated devel-opment environment (IDE) from ST, called STM Workbench, which is a modifi-cation to the freely distributed Eclipse IDE. Along with this workbench, follows also a compiler targeting the previously mentioned evaluation platform. This compiler is embedded into the STM Workbench IDE. ST also provides a firmware package for the evaluation platform, which has been used.
In order to test the algorithms, MATLAB from MathWorks has been used. MATLAB has also been used to develop and test the GUI, which is responsible for collecting the data sent from the evaluation platform. The use of MATLAB has been available through Link ¨oping University.
2
System overview
2.1
Concept
In order to measure vital signs, without any physical contact to the patient, a radar device will be implemented using the specified frequency of 5.8 GHz. The radar device will send out this signal towards the patient and then also collect the reflected signal. Upon hitting the patient and being back scattered, the origi-nal sigorigi-nal will be modulated with a phase that contains the relevant information about the vital signs. Eq. 2.1 shows the outgoing SLO and the phase-modulated SRF, both carrying the 5.8 GHz signal denoted as fLO.
SLO = A1cos(2π fLO)
SRF = A2cos(2π fLO+φ(t))
(2.1)
As can be seen in the equations, the reflected SRF contains an unknown phase called φ(t). This phase is expected to contain the vital signs gathered from the patient together with potential noise. Eq. 2.2 shows the expected contents of the unknown phase.
φ(t) =cos(2π fRR) +cos(2π fHR) +θ(t) (2.2)
Apart from θ(t), which is considered as noise, the two vital signs respira-tion and heartbeat rate, denoted fRRand fHR respectively, are hidden within this
phase. The frequency of the respiration rate and the heartbeat rate, could here be determined by translating the phase-signal from the time domain into the fre-quency domain. This translation could be performed by the use of Fourier trans-form.
2.1.1
The multi-port correlator
The receiver front-end using the multi-port technique is shown to the right in
Figure 2.1. This multi-port correlator is compared, for clarity of operation, to the
classical receiver front-end, to the left in Figure 2.1. The multi-port correlator con-sists of a Wilkinson power divider and three quadrature couplers, and is usually implemented with transmission lines on a Printed Circuit Board (PCB).
Figure 2.1: Translation of a quadrature mixer into a multi-port [4].
To perform the operation of a mixer, diodes and filters are connected at the four outputs of the multi-port correlator. Moreover, the use of the multi-port technique in conjunction with the diodes, will perform quadrature demodulation of the incoming SRF-signal by the means of the coherent SLO-signal.
The multi-port correlator will split the input signals SRF and SLO along ferent paths so that, at the output of the correlator, linear combinations with dif-ferent phase delays are generated. Through the squaring operation of the diodes followed by filtering of harmonics, the in-phase I and quadrature phase Q sig-nals are obtained, at the output of the multi-port front-end. Moreover, the sigsig-nals obtained are differential such that the in-phase signal consists of I−and I+while the quadrature phase signal consists of Q− and Q+. With the input signal SRF and SLO from Eq. 2.1, at the port P1 and P2 in Figure 2.1, the output at ports P3 and P4 can be written as Eq. 2.3, in phasor-form.
P3signal = A2 2 e i(2π fLO+φ(t)−180)+ A1 2 e i(2π fLO−360) P4signal = A2ei(2π fLO+φ(t)−270)+ A1ei(2π fLO−270) (2.3)
In a similar way, the output at ports P5 and P6 can be written as Eq. 2.4. P5signal = A2 2 e i(2π fLO+φ(t)−270)+ A1 2 e i(2π fLO−180) P6signal = A2 2 e i(2π fLO+φ(t)−180)+ A1 2 e i(2π fLO−270) (2.4)
After performing the squaring operation using the diodes along with filtering, the expected output signals are shown in Eq. 2.5.
I− = −A1A2 4 cos(φ(t)) I+ = A1A2 4 cos(φ(t)) Q− = −A1A2 4 sin(φ(t)) Q+ = A1A2 4 sin(φ(t)) (2.5)
The output signals can be further processed, either as differential signals or as single-ended signals. A differential to single-ended amplifier can be used to accomplish the signals in Eq. 2.6.
I = A1A2
2 cos(φ(t))
Q = A1A2
2 sin(φ(t))
(2.6)
At this point, the 5.8 GHz carrier has been removed through down conversion of the SRF-signal to the base band. This ideally means that only the phase, that
includes the information about the vital signs remain.
2.1.2
Demodulation
The quadrature signals at the output of the multi-port front-end needs further demodulation in order to extract the phase from the signal. This can be done by the use of the arc-tangent function as shown in Eq. 2.7.
φ(t) = arctanQ
I ⇔ φ(t) = arctan
sin(φ(t))
This would result in the unknown phase φ(t), being represented as the super-position signal described in Eq. 2.2. However, in order to get a valid φ(t), both the quadrature signals have to be collected from the output of the multi-port front-end at the same time.
The acquisition-part
The collection of this data is done by the acquisition part of the sensor platform. Since the acquisition part also is the connection to the DSP-part of the platform, this will be some sort of Analog to Digital Converter (ADC). Furthermore, since there are at least two outputs from the multi-port front-end, there has to be at least two separate ADCs in order to collect the data simultaneously.
The DSP-part
When the ADCs have collected a pair of samples from the output of the multi-port front-end, this will be sent to the DSP-part. The DSP-part will then demod-ulate the incoming samples according to Eq. 2.7 to get the desired φ(t). The demodulated phase will finally be translated into the frequency domain with the use of Fast Fourier Transform (FFT) in order to get a spectrum of the frequencies that the signal consists of. The result from the FFT, will be transferred through a communication link, to a computer or server where the data will be presented to the user in plot.
2.1.3
The sensor platform
To summarize the above sections and theory into a graphical perspective, the block diagram in Figure 2.2 was made. The block diagram briefly describes the sensor platform in terms of hardware.
2.2
Task
As mentioned, the task of this thesis is to implement the acquisition- and DSP-part of the sensor platform, together with an interface that can present the final data to the end user.
The acquisition part of the sensor platform should include two separate ADCs for sampling of the I and Q signals of the multi-port. These two signals have to be sampled simultaneously in real time, to not get skew in the data. Furthermore, these two ADCs have to have enough resolution and enough speed to give accu-rate samples. The DSP-part of the sensor platform should include the following features:
• Filtering and decimation
• Arc-tangent calculation between the I and Q signals to determine the phase • Fast Fourier Transform (FFT) to be able to determine the frequency of the
incoming signals
• Peak detection and windowing of data
• Communication link (Ethernet/USB) to computer or server
The flowchart in Figure 2.3 shows the expected flow of data from ADCs to user, with the previously mentioned features to the DSP-part.
2.3
Hardware
The first part of this thesis was put into deciding what kind of hardware that should be used in order to implement the acquisition- and DSP-part described in the previous section. The most important points included the need of at least two ADCs, either Universal Serial Bus (USB) or Ethernet connectivity for the commu-nication link and finally enough processing speed to perform the rather heavy FFT.
Together with the participants of the research project, the decision was made to use an existing evaluation platform such as the STM32 Nucleo-F767ZI. The discussion had previously been to use an Field Programmable Gate Array (FPGA) system as the DSP-algorithms, in order to efficiently build a pipe lined structure that could process data very fast. This was however rejected mainly because of the high price of such a device, but also because this particular evaluation platform included many of the peripherals needed.
2.3.1
The Nucleo board
The Nucleo development board contained a STM32F767ZI Micro Controller Unit (MCU), which is essentially an ARM 32-bit Cortex-M7 [5]. Furthermore, it con-tained both USB and Ethernet connectors which were possible to program and setup through the MCU. ST, who has developed the platform, also provided a programming IDE based on Eclipse, which made it possible to program and also debug the code with minor additional setup steps.
The MCU itself included several needed perks for this project. First of all, the processor core could run at 216 MHz and had dedicated DSP instructions, which would make the FFT calculations among other things more efficient. The amount of available memory on the device was up to 2 MB of flash together with 512 kB of SRAM. The MCU also had three separate ADCs which could be configured to work in parallel. These ADCs could also be configured to use the on-board Direct Memory Access (DMA) controller in order to reduce the number of calls to the processor.
Finally, since the board was an evaluation platform, it included an on board USB-to-serial controller acting as a debugging interface between the computer and the board. This allowed the board to be powered from the USB-port of a computer, regardless if the board was being debugged or not.
2.3.2
The ADCs
As previously mentioned, the STM32F767ZI MCU contained three separate ADCs. All of these had a maximum resolution of 12-bits, but could be configured for lower resolutions if necessary. The result from each conversion was stored into a 16-bit data register, from where it then could be collected either by the processor itself or by a DMA-unit. Each one of these ADCs allowed a maximum sample rate of 2.4 MSa/sec, which was considered more than enough for this application.
The ADCs could also be configure in dual or triple mode which allowed them to work in parallel. This kind of configuration would also increased the total amount of Sa/sec. Together with the DMA-unit it was also possible to generate user-defined interrupts such as on each conversion from the ADCs or when the DMA-unit has collected a certain amount of samples. By using the second type of user-defined interrupts it would be possible to buffer data into larger chunks of data that all could processed at once.
2.3.3
The Ethernet controller
One of the many features this Nucleo board had, was the Ethernet connection port which allowed the device to be connected to an Ethernet network. This essentially meant, that when properly configured, the device could send and re-ceive data from the internet or a local area network, which was likely to exist at or close to the installation of the future finished device.
The Ethernet controller supported 10M/100M communication speed, and was connected to the MCU through a Reduced Media-Independent Interface (RMII).
ST did not provide any own software when using the Ethernet device, but was
instead relying on the lightweight TCP/IP protocol suite, that was made available through LwIP-project [6]. This suite was included as a third-party software in the general firmware package that ST made available for the Nucleo board.
3
Data acquisition
3.1
Sampling
As previously mentioned, the vital signs that should be monitored was the heart-beat and the respiration rate. Both of these signals are relatively low in frequency. The heartbeat rate varies between different persons and depends on several factors. However, a general rule for the maximum heartbeat rate of a person is shown in Eq. 3.1.
HEARTBEATMAX =220−age, where age is years (3.1)
This sets the highest frequency for the heartbeat to about 3.6 Hz, when dis-regarding the age. On the contrary, the lowest frequency of the heartbeat goes towards zero [7].
The respiration rate of a person is also varying between persons, but in gen-eral an adult person at rest takes about 12 breaths per minute [8]. This rate will increase with increased activity, however this will never reach the fast pace of the heart. Therefore, the low rate of 12 breaths per minute can be used to set the lower limit of the frequency span, where the vital signs is expected to be. This span reaches from 0.2 Hz, from the respiration rate, to 3.6 Hz, from the heart rate. According to the Nyquist sampling criterion, the minimum sampling frequency needed, in order get a correct representation of the analog signal in the digital do-main, must be two times the highest frequency we want to represent. This would mean that, this system needs to gather at least 7.2 Sa/sec, to fulfill the criterion.
The three ADCs on the Nucleo board could sample at a maximum speed of 2.4 MSa/sec. Since this sampling speed was more than enough to fulfill the pre-viously mention criterion, the real challenge would be to reduce the sampling speed to a somewhat more reasonable amount. However, the low sampling speed stated by the criterion, was not ideal either. Therefore, in order to allow filtering
and decimation in the pre-processing steps, the sampling rate was decided to be 2048 Sa/sec.
3.2
Implementation
As previously mentioned, the MCU mounted on the Nucleo board ships with three ADCs that can run in parallel. This allows three samples to be gathered at the same time.
This project was designed with the intention of gathering two samples at the same time, which then also means that only two ADCs are in use. Gathering these two samples at the same time, is however crucial because of the arc-tangent demodulation.
The ADCs can be used together with a DMA in order to offload the processor when continuously sampling. The DMA could then be configured to alert the processor by triggering an interrupt, when a certain amount of samples has been collected. These samples could then be gathered by the processor and put into the algorithm. This way would also allow the processor to work with blocks of data rather than single samples, which is needed for the FFT-algorithm.
3.2.1
Defining the sampling clock
One of the most important steps when configuring the ADC, was to define the clock to which each ADC sample would be triggered. This would essentially determine the sample rate of the ADCs and would need to be low enough, in order to reach the desired sample rate.
The ADC-clock
The ADC-modules includes an ADC-clock which among other things, sets the time used to collect each sample. This means, the ADC-clock sets the sample time used during a single measure by the ADC. If the ADCs would be configured to continuously sample, this clock also sets the sample rate.
According to the reference manual of the STM32F767ZI [9], the ADC-clock was derived from the APB2 clock through a programmable pre-scaler. The APB2 clock was in turn derived from the system clock, which also could be configure through a programmable pre-scaler. This defined the ADC-clock as described in
Eq. 3.2. However, according to the datasheet of the STM32F767ZI, the maximum
ADC-clock was 36 MHz, which makes the equation only valid for ADC-clocks that are below this value.
fADC = fSYSCLK
APB2DIV IDER·ADCDIV IDER (3.2)
The conversion time needed for the ADC to make one measurement, was firstly determined by the resolution used by the ADC. This resolution was set
to the maximum of 12-bits as previously stated, which then resulted in 12 ADC-clock cycles, as each bit needed one ADC-clock cycle to process the value. Secondly, the ADCs could be configured with a sampling time, which would define how long the ADC would allow the input capacitor to charge. The lowest sampling time allowed by the processor was three ADC-clock cycles while the longest sampling time was 480 ADC-clock cycles. The sample time for one ADC conversion could then be calculated using Eq. 3.3
tconversion = tS·NR OF BITS
fADC (3.3)
If the ADCs would be configured for continuous sampling, the conversion time would equal to the time period for the sample rate. This means the sample rate would be the inverse of the conversion time, as in Eq. 3.4.
fS = 1
tconversion (3.4)
By using this relation together with the two previous equations, the maxi-mum and minimaxi-mum sample rate for a 12-bit ADC, using the ADC-clock can be summarized as in Table 3.1.
Table 3.1: The sampling characteristics for the ADCs
APB2DIV IDER ADCDIV IDER fADC fs
Max. sample rate ∼ ∼ 36 MHz 2.4 MSa/sec Min. sample rate 16 8 1.687 MHz 3.43 kSa/sec
Since the minimum sample rate according to Table 3.1 was only 3.43 kSa/sec, the ADCs can not be run in continuous mode and have to be triggered by an external source, like a timer, to reach the desired sample rate of 2048 Sa/sec.
Using a timer
The ADC-clock will still be used to define the conversion time for a single sample, since it is needed for the ADCs to work. Furthermore, since the lowest conversion time for a single sample still is faster than the desired conversion time, the ADC-clock will be configured according to the minimum sample rate in Table 3.1.
Triggering of the ADCs will be done by configuring a timer to run and gener-ate an interrupt with the desired frequency. The interrupt will then be forwarded to the ADCs, by configuring these to use an external event as trigger.
To configure the timer to work with an interrupt and the ADCs, only a few things needed to be done. Firstly, a timer that can be used to trigger the ADCs through the hardware, has to be used. Secondly, this timer has to be configured with a period that generates an interrupt with the desired frequency.
According to the reference manual for the STM32F767ZI, TI M2TRGO can be
used to trigger the ADCs. Furthermore, TIM2 which this refers to, has a 32-bit auto-reload counter which can be used to set the period without further scaling to the clock of the timer. Since this counter can support numbers of size equal to the system clock itself, the Eq. 3.5 shows how the period can be calculated.
Period= TI M2 CLK
DESIRED CLOCK SPEED (3.5)
Apart from setting the period and the counter mode to up, the only remaining thing to do was to configure the timer to generate a trigger signal when the period has been reached. Figure 3.1 shows the code used to configure the timer.
Figure 3.1: Sample code showing the configuration of the timer.
According to the figure, the period depends on a variable called
SystemCore-Clock, which is set during the configuration of clocks used by the device. In this
case the SystemCoreClock will be configured for the highest speed of the device, 216 MHz. The TIM2-module is connected to the APB1 bus, which in this case has a clock frequency of half the SystemCoreClock. The calculation for the timer period then matches that of Eq. 3.5.
3.2.2
Dual-mode
To measure the two separate input signals simultaneously, two separate ADCs will be used. Furthermore, these two ADCs must be configured to measure their input at the same time. This is done by configuring the ADCs in dual-mode, which is one of the available multi-mode configurations. The dual-mode should further be configured to use the regular simultaneous conversion, in order to achieve the desired configuration.
This dual-mode is essentially configuring the ADCs to work in a master-slave configuration, with ADC1 as master while ADC2 is slave. This also means that the master will automatically trigger the slave, when a conversion is initiated.
When a conversion is done, the two 12-bit values, one from each ADC, which are stored in a 16-bit vector (half-word) gets put together into a 32-bit vector. The 32-bit vector contains the ADC1 value in the lower 16-bits, while the ADC2 value is stored in the upper 16-bits. The resulting 32-bit vector can be read out from the conversion register connected to the master, ADC1.
3.2.3
Configure the DMA
As previously mentioned, one of the DMA-modules was used to offload the pro-cessor of the device. This would allow the propro-cessor to work without interruption for each conversion made by the ADCs. However, the DMA was configured to interrupt the processor after a specified amount of conversions had been done, which then allowed the processor to grab this data as a block.
According to the reference manual, only the DMA2-module was connected to the ADCs. Furthermore, out of the eight streams that the DMA2-module had, only the lower five had an actual channel to any of the ADCs. Since the two ADCs used were configured in dual-mode, only the master needed to be con-nected to the DMA. This resulted in the use of stream zero on DMA2, since this was connected to ADC1.
The DMA-buffer
Apart from configuring the DMA2-module to work with the ADCs, the interface towards the processor also had to be configured. This was done by defining a global memory space as a buffer, that the DMA could fill with data. The DMA was then configured to generate an interrupt for each half of the buffer. That is, an interrupt would be generated when the buffer was half full, and another inter-rupt would be generated when the buffer had been completely filled. In order to repeat this process each time the buffer has been filled, the DMA is configured to treat the provided buffer as circular. This option will force the DMA to start from the beginning of the buffer when it has reached the end.
3.2.4
Verification
To get the desired sampling rate for the data acquisition part of the algorithm, three peripherals has been linked together. The two ADCs that makes the mea-surement of the input signals are triggered by the timer TIM2, which has been configured for a frequency of 2048 Hz. When a conversion has been performed by the ADCs, the DMA2-module transfers this data into a global buffer without interference of the processor. The DMA2-module also generates two interrupts to alert the processor that a new chunk of data is available for processing.
To verify the functionality of the ADC, two measurements were done. The first one, estimates the actual trigger frequency generated by the timer, while the second one estimates the interrupt frequency of the DMA2-module. For the second measurement, the size of the global buffer has been set to two times the sampling frequency of 2048 Hz. This would in theory generate an interrupt
ev-ery second because of the relation to the sample rate. Table 3.2 summarizes the measurements compared to theory.
Table 3.2: ADC sampling and trigger frequencies
tINTERVAL f
Timer Calculated 488.28 µs 2048 Hz
Measured 488.3 µs 2047.9 Hz
DMA Calculated 1 s 1 Hz
4
Processing
4.1
Concept
As previously mentioned, the processing part of the project, refers to two major operations. First, the last part of the demodulation of the signal has to be done with the use of arc-tangent calculation, in order to extract the phase from the phase-shifted signal. Secondly, the frequency spectrum of the phase has to be cal-culated in order to determine which frequencies the phase consists of. When this has been done, a simple magnitude calculation to find the peaks of the frequency spectrum will remain before the signal can be plotted for the end-user.
The majority of these two operations has been implemented using the DSP-functions provided by the CMSIS library from ARM [10].
4.2
The frequency spectrum
The CMSIS library contains an example implementation of how to calculate the frequency spectrum using their predefined transform, magnitude and peak-detection functions. This example was used as a start point in the FFT implementation.
Table 4.1 contains the few requirements stated before starting the
implementa-tion of the FFT.
Table 4.1: Considered requirements for the FFT-implementation
Description Requirement
Number of points 1024 (minimum) Frequency resolution [Hz] 0.1 Hz
4.2.1
FFT
The CMSIS library function for the FFT, supports lengths from 16, 32, 64 etc up to a maximum of 4096 points. This made the first requirement fairly simple to reach. The library also provides predefined constants to setup the correct FFT-configuration. Furthermore, the library function operates on complex data and also returns complex data.
Figure 4.1 shows the initial FFT-implementation. The arm cfft f32-call is the call
to the CMSIS-library function for calculating the FFT-spectrum. The first input-argument to this function defines the setup for the FFT, as can be seen in the figure. The second argument is the data that should be transformed into the frequency domain, while the two last arguments defines how the function should expect this data.
Figure 4.1: Sample code showing the FFT-implementation.
To get the correct spectrum in a human readable format, the complex magni-tude also had to be calculated in order to convert the complex output vector from the FFT into real values. This is done through the function call to arm cmplx mag f32, which also only returns a vector of half the size of the input.
Frequency resolution
In order to present a valid spectrum for the incoming signals, the scaling had to be correct. As previously mentioned, the expected signals had a very low frequency (less than 10 Hz) which made it reasonable to only show a spectrum in that region.
The previously stated requirement for a frequency resolution of 0.1 Hz, then implies that each bin of the FFT-result contains an increase of 0.1 Hz frequency. In order to reach this, the sample-rate needs to be adjusted with regards to the length of the FFT. Eq. 4.1 shows the simple relation between sampling frequency and FFT length [11][12].
FFT resolution= Fs
N, where N= FFT SIZE (4.1)
This suggests that the sample rate needs to be one tenth of the length of the FFT to reach the required resolution. As mentioned earlier, the sampling goal for the
ADC:s were set to 2048 samples per second which made it reasonable to increase the length of the FFT to 2048. This would by default then mean that the FFT-resolution would be 1 Hz and also easier to reduce with decimation of the input samples after the filtering has been done.
Testing the FFT
To verify that the FFT worked as expected, a reference signal was created in MAT-LAB and then written to a header-file that could be compiled with the C-code. The signal generated was a super-position of the four frequencies 10 Hz, 40 Hz, 50 Hz and 500 Hz. The equation used to calculate the signal is seen in Eq. 4.2. The vector generated from this signal contained 2048 real samples, equal to the lenght of the FFT, and 2048 imaginary values simply set to zero.
superPosSign=0.5(sin(a f0) +...+sin(a f3)), where a=2πt (4.2)
To verify the output spectrum from the DSP, the same signal was also plotted using only MATLAB in order to get a reference spectrum. Figure 4.2 shows the reference spectrum from MATLAB. Here you can see a total of four peaks at the correct frequencies, and nothing else.
Frequency [Hz] 0 200 400 600 800 1000 1200 Magnitude [V] 0 0.1 0.2 0.3 0.4 0.5 0.6 FFT Spectrum (only FFT) X: 500 Y: 0.5 X: 50 Y: 0.5 X: 40 Y: 0.5 X: 10 Y: 0.5
Figure 4.2: Reference spectrum calculated in MATLAB of the super-position signal.
In comparison the calculated spectrum using the DSP can be seen in Figure
4.3. Comparing the two plots lead to the conclusion of the FFT-functionality in
Frequency [Hz] 0 200 400 600 800 1000 1200 Magnitude [V] 0 0.1 0.2 0.3 0.4 0.5 0.6 FFT spectrum (DSP) X: 500 Y: 0.5 X: 50 Y: 0.5 X: 40 Y: 0.5 X: 10 Y: 0.5
Figure 4.3: Spectrum calculated using the DSP-FFT functions from ARM.
4.2.2
Peak detection
From a user perspective, the full frequency spectrum is not interesting more than for debugging purposes. The user wants fast and accurate information about the heartbeat and respiration rate, and therefore there is a need of peak detection from this spectrum.
In Figure 4.1, the last function call to arm max f32 calculates the highest peak in the signal. However, there is at least two peaks that needs to be found.
Fig-ure 4.4 shows the code used to find multiple peaks in the spectrum. The
num-ber of peaks to find is specified at compile time of the code, through the macro
NR OF PEAKS.
4.3
Demodulation with arc-tangent
Before the frequency spectrum of the input signal can be calculated, the last part of the demodulation has to be done, by extracting the phase from the incom-ing signal. This was done by usincom-ing the trigonometric function arc-tangent. The CMSIS-library however, did not contain this function and therefore it had to be implemented.
This was done by implementing the fourth order Taylor series expansion of arc-tangent that can be seen in Eq. 4.3 [13]. However this series is only valid for t ≤ 1 because above this it will start to diverge. Also the accuracy for this expansion at t=1 is limited. arctan(t) = t−t 3 3 + t5 5 − t7 7 + t9 9 (4.3)
Figure 4.5 shows the accuracy issue at t = 1 for three different orders of the
Tay-lor expansions. The comparison was made towards the MATLAB implemented
atan2. t 0 0.2 0.4 0.6 0.8 1 1.2 phi [rad] 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9
1 Comparing atan2 to Taylor expansion
atan2 Taylor O(2) Taylor O(3) Taylor O(4)
Figure 4.5: Comparing the MATLAB atan2 to three different orders of Taylor expansions for arctan.
The error for each of the orders can be seen in Figure 4.6 which is a zoomed in version of the previous figure. It is obvious here that increasing order of the Tay-lor expansion gives a smaller error, which corresponds to the series being infinite for a correct result. However, an infinite series is not practical in the program-ming world and therefore the fourth order was used.
t 0.65 0.7 0.75 0.8 0.85 0.9 0.95 1 1.05 1.1 phi [rad] 0.6 0.65 0.7 0.75 0.8 0.85 0.9 0.95
1 Comparing atan2 to Taylor expansion [zoomed]
atan2 Taylor O(2) Taylor O(3) Taylor O(4) X: 1 Y: 0.8667 X: 1 Y: 0.8349 X: 1 Y: 0.7854 X: 1 Y: 0.7238
Figure 4.6: Showing the error at t=1 compared to the atan2.
The Eq. 4.3 was not practical to implement into the code in its current form, therefore it was re-written into Eq. 4.4.
arctan(t) =t(1−t2(1 3 −t 2(1 5−t 2(1 7− t2 9)))) (4.4) As can be seen this equation now only has a t2-term instead of a t7-term. This made it possible to calculate t2one time and then simply iterate it for each bracket. Furthermore, the negative sign can be applied when calculating the t2-term, mak-ing the code only contain addition [14].
4.3.1
Adjusting the error
Eq. 4.4 contains a total of five constants that can be pre-calculated in order to get
rid of their respective division. Also, if Figure 4.5 is studied, the error only starts to occur closer to t = 1 for all the different orders of the Taylor expansion. This indicates that the smaller constants could be adjust to reduce the error at t = 1.
Table 4.2 contains these new modified constants.
Table 4.2: Modified constants for arctan approximation
Order Value Approximation
1 1 1.000 (1/1) 2 1/3 0.33316 (1/3) 3 1/5 0.19048 (4/21) 4 1/7 0.09916 (12/121) 5 1/9 0.02723 (7/257)
As can be seen, the first two constants are basically the same as before, while the three following constants has been adjusted. The result of this adjustment can be seen in Figure 4.7.
t 0 0.2 0.4 0.6 0.8 1 1.2 phi [rad] 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9
1 Comparison of atan2 to Taylor expansion and Taylor approximation
atan2 TaylorApprox Taylor O(4)
Figure 4.7: Comparing the adjusted Taylor expansion to the old one and to atan2.
As can be seen in the figure there is no noticeable error compared to atan2 with the modified constants. And compared to the standard constants for the Taylor expansion, there is a rather large improvement.
Implementing in code
When implementing Eq. 4.4 into code, an iteration from the inner most bracket and out was made. Figure 4.8 show the equation written in C-code, starting by adding the inner most constant to the output angle phi.
Figure 4.8: The implemented code for arctan approximation in the first octant.
Notice that this code expects two inputs, x and y, in order to first calculate t. After that, t2can be calculated and the iteration can start in an un-rolled loop by adding the first constant.
4.3.2
Rotating into first octant
As previously mentioned the Taylor expansion in Eq. 4.4 was only valid for t≤1. In the unit circle, this is only true in the first octant where x ≥ y. However, the
input values could be anywhere on the circle, and therefore some rotation of the coordinates was needed.
Figure 4.9 displays a unit circle cut into the four quadrants with the first
quad-rant divided into two octants.
π 2 y π 3π 2 0 x π 4
Figure 4.9: The unit circle divided into four quadrants.
The coordinate rotation applied uses a maximum of three rotations to move any coordinate into the first octant. It is implemented by simple if-statements that check which of x and y is greatest and rotate according to that. The three statements check whether the coordinate is in quadrant three or four, quadrant two or if the coordinate is in octant two. Table 4.3 show the rotation performed for each of these three statements.
Table 4.3: Rotation of coordinates in the unit circle
Coordinate position Condition Rotation
3 or 4 [Quadrant] y<0 180◦
2 [Quadrant] x≤0 90◦ 2 [Octant] x ≤y 45◦
These three statements were applied into the C-code in a separate function called octantify, which can be seen in Figure 4.10. The resulting number of octants that the coordinate has been rotated with is returned in the parameter octant.
Figure 4.10: The implemented code to rotate any coordinate into the first octant.
During the last rotation in the first quadrant, that is from octant two into oc-tant one, there is a scaling that gets added to the coordinates during the rotation. The scaling factor is 1/√2, however this scaling is added to both x and y and will be eliminated during the calculation of t, in the first steps of arctan-calculation.
4.3.3
The Taylor approximation in code
The full Taylor approximation of arc-tangent, implemented in the C-code can be seen in Figure 4.11. The function starts with a simple test to see whether we have a simple case that does not need any calculation.
Figure 4.11: The implemented code for approximation arctan with Taylor expansion.
After rotation the coordinate, the actual rotation angle is saved in radians in the variable angleShift and later added to the calculated angle before the function returns.
Also notice the last if-statement of the code, which checks if the angle is larger than π. This is due to the arc-tangent function by definition is defined for [−π,
4.3.4
The frequency spectrum with a phase-shifted signal
The fully implemented arc-tangent was tested in a similar way to the previously mentioned FFT. Furthermore the frequency spectrum was calculated with the DPS-implemented FFT before plotting in MATLAB. As a reference, the corre-sponding spectrum were calculated in MATLAB, both by using atan2 and the approximated Taylor expansion.
As for input signal to the DSP, a phase-shifted signal containing the previously used super-position signal, was created in MATLAB and transferred in a header-file to the DSP. The base signal was generated in MATLAB using Eq. 4.5 with the input from Eq. 4.2, and then transformed into a complex I- and Q-signal using the Hilbert-transform function [15].
phaseShi f tSignal =sin(superPosSign) (4.5)
The reference output spectrum created in MATLAB can be seen in Figure 4.12, were four distinct peaks can be seen. These also corresponds to the four input frequencies in the super-position signal.
Frequency [Hz] 0 200 400 600 800 1000 1200 Magnitude [V] 0 0.2 0.4 0.6 0.8 1 1.2 FFT Spectrum (atan2) X: 500 Y: 0.8543 X: 50 Y: 1.028 X: 40 Y: 0.9189 X: 10 Y: 0.8807
Figure 4.12: The reference frequency spectrum for the phase-shifted signal.
The corresponding frequency spectrum generated in MATLAB but for the Taylor expansion can be seen in Figure 4.13. This spectrum also has the four fre-quencies in distinct peaks, as the reference.
Frequency [Hz] 0 200 400 600 800 1000 1200 Magnitude [V] 0 0.2 0.4 0.6 0.8 1 1.2 FFT Spectrum (TaylorApprox) X: 500 Y: 0.8543 X: 50 Y: 1.028 X: 40 Y: 0.9189 X: 10 Y: 0.8807
Figure 4.13: The frequency spectrum for the phase-shifted signal, calculated with the Taylor approximation in MATLAB.
Finally Figure 4.14 shows the frequency spectrum calculated on the DSP using the Taylor series expansion approximation for arc-tangent calculation.
Frequency [Hz] 0 200 400 600 800 1000 1200 Magnitude [V] 0 0.2 0.4 0.6 0.8 1 1.2 FFT spectrum (DSP) X: 500 Y: 0.8543 X: 50 Y: 1.028 X: 40 Y: 0.9189 X: 10 Y: 0.8807
Figure 4.14: The frequency spectrum for the phase-shifted signal, calculated on the DSP.
As for the frequency spectrum for the super-position signal, this spectrum is identical to the reference spectrum. The spectrum shows all the peaks in correct position, but compared to only the super position signal, all these spectrum’s have an increase in amplitude of each peak, which should be centered at 0.5 V.
5
Pre-processing
5.1
Concept
Before the data is processed by the arc-tangent demodulation and the FFT calcu-lation, potential unwanted components in the signal should be removed. This is done by filtering the signal through a digital filter. Furthermore, since the sampling rate of 2048 samples per second is quite high, the signal can be down-sampled by the use of decimation. But since the FFT-module required a certain amount of samples before a valid result could be generated, this decimated data had to be stored in wait for further processing.
This chapter describes the steps performed to prepare the data for processing in the FFT-module. This includes the description of storing data between steps, and the generation of the sliding window effect.
5.2
Filtering
As previously mentioned, the heartbeat and respiration rate are very low fre-quency signals, typically less than 5 Hz. This leads to the conclusion, that a low-pass filter should be used in order to remove higher frequency components that might affect the resulting frequency spectrum.
The CMSIS-library offered several type of base filters, however two of these was of particular interest, namely Finite Impulse Response (FIR) and Infinite Im-pulse Response (IIR). In general IIR-filters are more efficient and requires less memory, compared to FIR-filters, when used in DSP applications. This is due to FIR-filters having a frequency resolution similar to the FFT, but also because IIR-filters can be implemented in a recursive manner. In this particular applica-tion, the FIR-filter would have required a very long coefficients vector in order to achieve a good resolution at lower frequencies. Therefore, the filter implemented was decided to be of IIR-type with Butterworth coefficients to achieve a frequency response close to a box.
5.2.1
Generating coefficients
In order to implement a Butterworth filter, using the CMSIS-library, a Biquad
Cas-cade Direct Form II-transposed structure was used with coefficients generated in
MATLAB. These coefficients were then written to a separate header-file and in-cluded in the C-code implementation.
According to the CMSIS documentation of the Biquad Cascade Direct Form II-transposed structure, the Eq. 5.1 is applied for each stage of the filter [16].
y[n] = b0x[n] +d1
d1 =b1x[n] −a1y[n] +d2
d2 =b2x[n] −a2y[n]
(5.1)
And, as can be seen, each stage can as maximum represent a second order filter, so in order to achieve a filter of higher order, several stages are necessary. This gives
a1..a2 and b0..b2 as the coefficients that needs to be generated for a Butterworth
filter of second order. Furthermore, according to Eq. 5.1 the a1and a2coefficients
needs to be negated in order to give a correct result.
In MATLAB a filter of third order was generated using the function butter(). This resulted in a total of ten coefficients as can be seen in Table 5.1, and also that the filter then consists of two stages with one of the stages having a d2 =0.
Table 5.1: Butterworth coefficients for the filter
Coefficient Value a11 0.8706641 a12 0.0 a21 1.8534351 a22 −0.8712398 b10 0.0002878 b11 0.0002878 b12 0.0 b20 1.0 b21 2.0000105 b22 1.0000105
The filter was designed for a cut-off frequency of 45 Hz, to remove high fre-quency components but could possibly be set at an even lower cut-off frefre-quency if necessary.
Figure 5.1 shows the frequency response of the generated filter in MATLAB.
Frequency (Hz) 100 101 102 103 Resoponse (dB) -300 -250 -200 -150 -100 -50 0
Frequency response for Butterworth order 3 X: 45 Y: -3.01
Figure 5.1: Frequency response of the implemented Butterworth filter.
5.2.2
Implementing the filter
In the C-code, the implementation of the filter was similar to the implementation of the FFT. The CMSIS library offered a function to filter the data as can be seen in Figure 5.2 as the call to arm biquad cascade df2T f32-function.
Figure 5.2: The implemented filter in C-code.
The first argument of this function, the parameter butterworth f32, defines the structure of the filter. In the init-call, this structure is initialized to use, the number of stages specified in the macro N STAGES, the filter coefficients provided and also a state buffer to record the state d1and d2from Eq. 5.1. The filter coefficients
was provided in a linear buffer as shown in Table 5.2 and would follow the same pattern for any higher order filter.
Table 5.2: Butterworth coefficients sorted into a linear buffer
Index 0 1 2 3 4 5 6 7 8 9
Coefficient b10 b11 b12 a11 a12 b20 b21 b22 a21 a22
Testing the filter
To test the implementation, the super position signal generated (Eq. 4.2) for the FFT test was used. Furthermore, the filtered signal was returned both as a time
domain signal and as a frequency spectrum, using the FFT implementation.
Fig-ure 5.3 shows the original time domain signal, that is the super position signal
from Eq. 4.2. Time (t) 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 Voltage (t) -2 -1.5 -1 -0.5 0 0.5 1 1.5
2 Super position signal
Figure 5.3: The super position signal in the time domain.
As previously has been said, the super position signal contained a total of four frequencies, where one of these frequencies was significantly higher than the other ones. By testing the filter design in MATLAB, the filtered reference in
Figure 5.4 was obtained. As can be seen, the high frequency component of the
signal has been removed.
Time (t) 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 Voltage (t) -0.6 -0.4 -0.2 0 0.2 0.4
0.6 Super position signal filtered
Figure 5.4: The filtered result of the super position signal, in the time domain.
The result from the DSP-implemented filter in Figure 5.5, showed that the im-plementation was working correctly.
Time (t) 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 Voltage (V) -0.6 -0.4 -0.2 0 0.2 0.4 0.6 Filtered signal (DSP)
Figure 5.5: The filtered result from the DSP, in the time domain.
This result was also further validated by looking at the output frequency spec-trum from the DSP, in Figure 5.6. Where the high frequency peak at 500 Hz has been removed from the signal, and the 50 Hz peak has been weakened.
Frequency [Hz] 0 200 400 600 800 1000 1200 Magnitude [V] 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 FFT spectrum (DSP) X: 50 Y: 0.3459 X: 40 Y: 0.3501 X: 10 Y: 0.3531
Figure 5.6: The filtered frequency spectrum from the DSP.
It can be noted in the spectrum that the amplitudes of the peaks has been reduced by the filter compared to the non-filtered peaks seen under the
Process-ing-chapter for the super position signal. This was however not considered a
