BEN-GURION UNIVERISTY OF THE NEGEV FACULTY OF ENGINEERING SCIENCES DEPARTMENT OF ELECTRICAL AND COMPUTER ENGINEERING
Detection of One Lung Intubation Incidents
Thesis Submitted in Partial Fulfillment of the Requirements for the M.Sc. Degree
LIOR WEIZMAN August 2004
ii
BEN-GURION UNIVERISTY OF THE NEGEV FACULTY OF ENGINEERING SCIENCES DEPARTMENT OF ELECTRICAL AND COMPUTER ENGINEERING
Detection of One Lung Intubation Incidents
Thesis Submitted in Partial Fulfillment of the Requirements for the M.Sc. Degree
LIOR WEIZMAN
Supervised by: Dr. Joseph Tabrikian and Prof. Arnon Cohen
Author signature
:
Date:
Approved by the advisors
:
Date:
:
Date:
Approved by the Chairman of the Graduate Studies Committee
iii
Ben-Gurion University of the Negev Faculty of Engineering Science Electrical and Computer Engineering Department
Detection of One Lung Intubation Incidents
By: Lior Weizman Advisors: Dr. Joseph Tabrikian and Prof. Arnon Cohen
Thesis Submitted in Partial Fulfillment of the Requirements for the Degree of Master of Sciences in Engineering. August 2004
Abstract During general anesthesia, voluntary breathing is stopped and the patient needs to be ventilated. The most common and the only method that suits long-term period of anesthesia is ventilation using an endotracheal tube, which is placed in the trachea. The location of the tip of tube is critical: it should be placed, and maintained above the bifurcation. A correct position of the tube, in which both lungs are ventilated, is called Tracheal Intubation (TRI). If the tube is misplaced or shifted due to patient movements, cases of One Lung Intubation (OLI) may occur. Prolonged cases of OLI should be avoided since it may cause insufficient oxygenation and may damage the non-ventilated lung. OLI was found to be the most probable cause of
iv
desaturation and the fourth cause of malfunction during anesthesia, and there is currently no reliable mean to detect OLI situation. This work offers a new method to detect OLI by analyzing the lungs sound signals. The use of lungs sounds for monitoring and diagnosis of pulmonary function is well known, but earlier works on OLI detection by monitoring lungs sounds have not provided reliable methods under real constraints. The proposed method assumes a MIMO (Multiple Input Multiple Output) system, in which a multi-dimensional AR (Auto-Regressive) model relates between the input (lung signals) and the output (recorded sounds). The ML estimator for the unknown parameters of the model is derived, and it has been shown that the resulted solution is actually an extension of the SISO (Single Input Single Output) LPC (Linear Prediction Coding) estimation to the MIMO case. In the theoretical development of the solution, information theoretic criteria, AIC and MDL, are used to estimate both the number of sources and the AR model order in the selected convolutive mixture model. In the simulation, the developed algorithm successfully estimates both the AR model order and the number of sources. Because of the fact that the lungs are non-point sources as they were treated in the selected model, AIC and MDL were unable to estimate the number of ventilated lungs when applied to real breathing sound signals. Therefore, the method is modified to handle scattered sources case, which is more accurate to model the lung sources. A detector, based on the estimated eigenvalues of the source covariance matrix, is developed, in order to detect one lung ventilation situation. In order to evaluate the performance of the proposed algorithm for OLI detection, a database was collected in the main surgery room of Soroka University Medical Center Faculty of Health Sciences, Beer-Sheva, Israel. The database consists
v
of breathing sound signals recorded during the anesthesia stage at the beginning of the surgery. A total of 24 patients participated in the study, and their breathing sound signals were recorded in a surgery room in both situations: during correct ventilation, when the tip of the tube is placed above the carina, and during a situation of OLI when the tip of the tube is under the carina and only one lung is ventilated. Testing the algorithm on the collected database shows an OLI detection rate of more than 95% with probability of alarm of less than 5%.
vi
Acknowledgments
I wish to thank my supervisors, Dr. Joseph Tabrikian and Prof. Arnon Cohen for their valuable ideas, assistance, close supervision and support during the course of my research. I am also very grateful to Dr. Shai Teiman and Prof. Gabi Gurman from Soroka University Medical Center and Faculty of Health Sciences for their continuous guidance and support of the medical part of the work. Special thanks and appreciation is given to Dr. Alex Zlotnik for his major contribution in establishing the database, which included spending hours of his private time in the surgery room. Thanks are also extended to my colleagues, Yariv Amos, Rami Hagag, Dror Lederman, Yair Noam and Shai Yaakobovich for their help and assistance during my research and to Dr. Vladimir Lapidos for his assistance in maintaining the sampling system. I would also wish to express my appreciation to my family in Beer-Sheva, my aunts Bruria and Aliza and my uncles, Shalom and Shlomi for their help and support during all the years I’ve studied in Ben-Gurion University. I must say that I respect and do not take for granted their assistance. Finally, I wish to express my gratitude and love to my family: my brother, Zahi my sister Inbal and to my parents, without whom loving support this thesis might not have materialized.
vii
TABLE OF CONTENTS Chapter
Page
Glossary of Symbols and Abbreviations....................................................x List of Figures ......................................................................................... xii List of Tables ...........................................................................................xiv 1. Preface .................................................................................................1 1.1
Introduction..................................................................................................1
1.2
Goals..............................................................................................................3
2. Medical Background ..........................................................................4 2.1
Introduction..................................................................................................4
2.2
Breathing Physiology ..................................................................................6
2.2.1
The Ventilation Stage ............................................................................6
2.2.2
The Circulatory Stage ............................................................................6
2.3
Ventilation Under General Anesthesia .....................................................7
2.3.1
Intubation and One Lung Intubation (OLI) ...........................................7
2.3.2
Detection of One Lung Intubation .........................................................9
3. Summary of Previous Research on OLI Detection by Lungs Sounds Analysis.....................................................................................................13 3.1
Introduction................................................................................................13
3.2
Detection of OLI by LPC and PARCOR Features .................................13
viii
3.3
Detection of OLI Using Source Separation .............................................15
3.4
Motivation for Developing a New OLI Detection Method .....................17
4. Blind MIMO Identification and Information Theory as an Approach to Determine the Number of Sources ....................................18 4.1
Introduction................................................................................................18
4.2
Blind MIMO Identification Based on High Order Statistics .................18
4.3
Additional Methods of Blind MIMO Identification ...............................20
4.4
Akaike Information Criterion – (AIC) ....................................................21
4.5
Minimum Description Length – (MDL) ..................................................24
5. The Proposed Model for the Problem ..............................................25 5.1
Introduction................................................................................................25
5.2
Model Formulation ....................................................................................25
5.3
The Maximum-Likelihood Estimator of A and R...................................28
5.4
Implementing MDL and AIC Algorithms for Estimating the Number of
Signals and the Model order .................................................................................31
6. Simulation Results ............................................................................35 6.1
Introduction................................................................................................35
6.2
Results for Implementing the ML Estimator ..........................................35
6.3
Results for Implementing the AIC and the MDL Criterion ..................37
6.4
Simulation Results Under Non-Point Sources Assumption ...................42
ix
6.5
Simulation Conclusions .............................................................................46
7. Mismatches Between the Assumed and Real Models......................48 7.1
Introduction................................................................................................48
7.2
Description of the Mismatches Between Assumed and Real Models ....48
7.3
Modification of the Solution for the Real Model.....................................49
8. Implementation on Real Data ..........................................................52 8.1
Introduction................................................................................................52
8.2
Experimental Setup ...................................................................................52
8.3
Pre-Processing ............................................................................................53
8.4
Real Data Results .......................................................................................55
9. Summary and Conclusions...............................................................61 References ................................................................................................64 Appendix A: Database Description .........................................................69 Appendix A.1 Introduction ...................................................................................69 Appendix A.2 The Sampling System....................................................................69 Appendix A.3 The Data Acquisition.....................................................................70
Appendix B: Proof of Mathematical Identities.......................................73
x
Glossary of Symbols and Abbreviations Abbreviations
AIC
Akaike Information Criteria
AR
Auto-Regressive
ARMA
Auto-Regressive Moving-Average
LLI
Left Lung Intubation
LPC
Linear Prediction Coding
MDL
Minimum Description Length
MIMO
Multiple Input Multiple Output
ML
Maximum Likelihood
OLI
One Lung Intubation
PDF
Probability Density Function
RLI
Right Lung Intubation
SNR
Signal-to-Noise Ratio
SISO
Single Input Single Output
TRI
Tracheal Intubation
xi
Notations
x, y
scalar
x, y
column vector
X, Y
matrix
xT , XT
vector/matrix transposition
N
number of samples
M
AR model order
K
number of sources
L
number of sensors
p
number of free adjustment parameters
θ
parameters vector
f (Y / θ )
probability density function of Y given θ
xii
List of Figures Figure
Title
Page
Figure 2.1: Anatomy of respiratory system ..................................................................5 Figure 2.2: Tube location in One Lung Intubation ........................................................8 Figure 2.3: Absorption characteristics of Hb and HbO2 .............................................10 Figure 2.4: O2 Sat as function of
R IR
............................................................................11
Figure 3.1: Source separation block diagram .............................................................15 Figure 5.1: Block diagram of the MIMO AR model ..................................................26 Figure 6.1: NMSE of ML estimators as a function of N ..............................................36 Figure 6.2: Pe of estimating K as a function of N ........................................................38 Figure 6.3: Pe of estimating M as a function of N........................................................38 Figure 6.4: Log-likelihood, AIC penalty function and their sum as a function of M ..39 Figure 6.5: Log-likelihood, MDL penalty function and their sum as a function of M 40 Figure 6.6: Distribution of MDL and AIC decisions on K, N=300. (True K=2) .........41 Figure 6.8: NMSE of the estimated matrices under scattered sources assumption as a function of scattering parameters.........................................................................43 Figure 6.9: Pe of estimating K as a function of scattering level, ε for K=1................44 Figure 6.10: Pe of estimating K as a function of scattering level, ε for K=2. .............44 Figure 6.11: The eigenvalues of R versus the scattering level, ε , for K=1................45 Figure 6.12: The eigenvalues of R versus the scattering level, ε , for K=2.................45 Figure 7.1: The eigenvalues of R versus the scattering level, ε . ................................51 Figure 8.1: Experimental set-up – view from back (TRI shown) ................................53 Figure 8.2: The developed GUI for pre-processing operations ...................................54
xiii
Figure 8.3: Examples of the recorded signals – four microphones, both OLI and TRI. ..............................................................................................................................56 Figure 8.4: Second eigenvalue of Rˆ as a function of time. ........................................57 Figure 8.5: The histrogram of the second eigenvalue under TRI/OLI situations. .......58 Figure 8.6: The DET of a detector based on the second highest eigenvalue of the estimated R...........................................................................................................60 Figure A.1: Schematic diagram of the sampling system .............................................69 Figure A.2: Picture of the sampling system.................................................................70 Figure A.3: A patient’s breathings sounds are recorded during an experiment...........72
xiv
List of Tables Table
Title
Page
Table 3.1: Confusion matrix of classification by LPC features .................................14 Table 3.2: Confusion matrix of classification by PARACOR features .....................14 Table 3.3: Confusion matrix of classification by normalized energies ratio .............15
1
1. Preface 1.1
Introduction Ventilation of patients under general anesthesia or in intensive care, is
performed by an endotracheal tube, which is placed in the trachea. The location of the tip of tube is critical: it should be placed, and maintained above the bifurcation. A correct position of the tube, in which both lungs are ventilated, is called Tracheal Intubation (TRI). If the tube is misplaced or shifted due to patient movements, cases of One Lung Intubation (OLI) may occur. Prolonged cases of OLI should be avoided since it may cause insufficient oxygenation and may damage the non-ventilated lung. OLI was found to be the most probable cause of desaturation and the fourth cause of malfunction during anesthesia [1]. There are currently several means for OLI detection: •
The use of stethoscope is the simplest mean but it is non-automatic and requires high attention and its reliability is low.
•
The Capnograph, which measures the exhaled amount of CO2 was proved to be an unreliable method for detection of OLI [1].
•
Oxygen saturation, measured by pulse oximeter is one of the most widely used method today [2],[3] and the most reliable one, but its results are provided with latency of 2 to 5 minutes, which may be too long to prevent damage.
Although the above methods are used today to detect OLI, there is currently no reliable method, which provides an automated, real-time OLI detection.
2
The proposed method for OLI detection is based on analyzing the breathing sound signals. An algorithm for detection of the number of ventilated lungs from the recorded breathing sounds is developed according to “real world” constraints. In order to examine the suggested method and to check the possibility of developing a monitoring tool in the future, a database of recorded breathings was established. The database was composed of 24 patients and demonstrated both the OLI and TRI situations. Applying the method over the real data showed that it has the potential of becoming a useful non-invasive diagnostic system for OLI situations.
This work is arranged as follows: Chapter 1 is a preface with an introduction and a description of the goals of this work. Chapter 2 describes the medical and physiological background of the breathing and reviews the existing methods for detection of OLI situation. Chapter 3 reviews previous research on OLI detection by lungs sounds analysis. Chapter 4 gives a short background about blind MIMO identification and detection of number of signals using information theory. Chapter 5 describes the selected model for the problem and the developed estimators for the unknown parameters. Chapter 6 shows and discusses simulation results applied to the developed estimators from chapter 5. Chapter 7 presents the mismatches between the chosen model and the real one, and offers a modification of the solution in according to the “real world” constraints.
3
Chapter 8 discusses the implementation of the modified solution on real data and presents real data results. Chapter 9 summarizes the results and conclusions, and makes recommendations for future work.
1.2
Goals Previous research on OLI detection based on lungs sounds analysis showed
poor results when tested on human being breathing sound signals. The work’s main goal is to develop a reliable algorithm for detection of OLI based on respiratory sound signals, and to prove the concept of the method by testing it on a comprehensive database, which was established especially for this purpose. The database includes the breathing sound signals during real operation in a surgery room, therefore it demonstrates almost exactly the real world problem. We hope that this work will be the basis of developing a real-time monitoring system on the future.
4
2. Medical Background [1],[2] 2.1 Introduction The respiratory system is a gas exchanging system in which the body absorbs oxygen from the air and releases carbon dioxide out of the body. In order to perform the exchange that mentioned above, two systems are needed: The respiratory system and the circulatory system. The respiratory system is made of the following parts: A. The Nasal cavity – The upper chamber where air is introduced to the body, moistened and warmed. B. Pharynx – Located at the back of the mouth where air and food are passed downwards into the body. C. Larynx – voice box where the vocal cords are located and the airways begin. D. Trachea – rigid tube that connects the larynx with the bronchi (windpipe). E. Alveolus – tiny, thin-walled air sac at the end of the bronchiole branches where gas exchange occurs (plural – alveoli). F. Bronchus – The rigid tube that connects the larynx with a bronchial tree. G. Diaphragm –muscle at the base of the chest cavity that contracts and relaxes during breathing. H. Epiglottis – a flag of tissue that closes over the trachea when food is swollen. I. Pulmonary capillaries – small blood vessels that surround each alveolus. J. Pleural membranes – thin membranes that cover the lungs, separate them from other organs and from a fluid-filled chest cavity.
5
Figure 2.1: Anatomy of respiratory system [1]
6
2.2 Breathing Physiology [1] 2.2.1 The Ventilation Stage As mentioned before, the goal of the breathing is to bring oxygen to body tissues and to release oxygenation products. In inhalation, the diaphragm and intercostal muscles (those are the muscles between the ribs) contract and expand the chest cavity. Air flows in through the airways (from high pressure to low pressure) and inflates the lungs. In exhalation, the diaphragm and intercostal muscles relax and the chest cavity gets smaller. Air from lungs (high pressure) then flows out of the airways to the outside air (low pressure). As we breathe air it goes past the epiglottis and into the trachea. It continues down the trachea through the vocal cords in the larynx until it reaches the bronchi. From the bronchi, air passes into each lung. The air then follows narrower and narrower bronchioles until it reaches the alveoli, where gas exchange occurs.
2.2.2 The Circulatory Stage After the air reaches the alveoli, oxygen is transferred to blood. Within each air sac, the oxygen concentration is high, so oxygen diffuses across the alveolar membrane into the pulmonary capillary. Oxygen diffused from the alveoli to the capillaries binds to Hb and CO2 is released to the alveoli and out of the body. This exchange of gasses occurs rapidly (fractions of seconds). The carbon dioxide then leaves the alveolus when in exhalation and the oxygen-enriched blood returns to the left heart and from there pumped into the body.
7
2.3 Ventilation Under General Anesthesia [2] During anesthesia voluntary breathing is stopped and the patient needs to be ventilated. The patient can be ventilated by several methods: •
Using a special mask designed for ventilation and attached to the patient’s face.
•
Using a laryngeal mask, which placed inside the patient’s mouth and covers the entrance to the larynx.
•
Using a tube, which is placed in the trachea.
The last method (also called intubation) is the only method that suits a long-term period of anesthesia or a respiratory intensive care, and is practiced every day all over the world.
2.3.1 Intubation and One Lung Intubation (OLI) When there is a need to ventilate a patient, the tip of the tube should be above the carina (the point where the trachea splits to the right lung and the left lung). An intubation in which there is a correct placement of the tube is called Tracheal Intubation (TRI). If the end of the tube is placed under the carina, then only one bronchus is ventilated. This situation is called One Lung Intubation (OLI) or endobronchial intubation.
8 tube trachea carina
One Lung Intubation –
Tracheal Intubation -
endobronchial intubation
Correct position of the tube
(OLI)
(TRI)
Figure 2.2: Tube location in One Lung Intubation
Situation of endobronchial intubation (OLI) causes several problems: 1. Ventilation / Perfusion mismatch, which causes hypoxia (a situation where there is a decrease of oxygen supply to the brain even though there is adequate blood flow). 2. Ventilating only one lung may cause an atelectasis of the unventilated lung and over inflation of the ventilated lung, which receives more air than it can handle with. 3. Ventilating only one lung may cause changes in the partial pressure of the CO2, damage to the tissues, and even a cardiac arrest and death. One-Lung intubation can occur in the beginning of the anesthesia - right after the insertion of the tube, or during the surgery after the patient head is moved.
9
2.3.2 Detection of One Lung Intubation The most important phenomenon of OLI is decreasing in the oxygen level in the blood. Therefore, in many cases, monitoring the oxygen level in the blood warns of ventilation problems. There are several methods to estimate the oxygen level in the blood, such as measuring the partial oxygen pressure PaO2 and measuring the oxygen saturation in hemoglobin – O2 Sat. Today, O2 Sat is a standard measure in the OR. Methods to detect one lung intubation 1. The Stethoscope The most simple and available instrument of the physician is the stethoscope. We can detect OLI by listening to lungs and comparing between the two sides of the lungs. This action takes time and attention and requires experience, and therefore it is not done usually during the surgery. In a research that was done [3], from 79 cases of OLI, doctors used stethoscope in 6 cases to detect the situation, and none of them was detected. 2. Estimation of oxygen level in the blood a. The clinical method Decreasing in O2 saturation may be detected by observing the color of the skin and the mucosas (lips). When the HbO2 level decreases, a blue color appears on the skin and on the mucosas (like in a situation of strangulation). The accuracy of this method is very low. In addition, the light conditions in the surgery room not always enables to determine the color changes. Besides, the color changes are noticeable only when there is 5 gr Hemoglobin with CO2 (Unsaturated Hb).
10
b. The pulse oximeter [20] The idea behind this method is that different configuration of hemoglobin molecules such as Hb, HbO2 absorb light in different levels for different frequencies. As seen in the following figure, two different types of hemoglobin – the first one attached to oxygen (Oxyhemoglobin) and the other one does not (Deoxyhemoglobin), absorbs the light in different levels in different frequencies.
ז
Figure 2.3: Absorption characteristics of Hb and HbO2 [20]
Hb absorbs more light in the red range, while HbO2 absorbs more light in the infrared range. Most often wavelengths 940nm and 660nm are used. Blood that flows in pulses in the arteries is mainly oxygenated blood (HbO2), it is marked as AC, and the blood that flows in the veins is mainly unoxygented blood (HB) and it is marked as DC. Let us denote the while absorption components for the AC blood in the different wavelengths by AC660, 940, and the absorption components for the DC blood in the different wavelengths by calculating DC660,940.
11
R = The ratio: IR
AC660 DC660 AC940 DC940
allows us to calculate the percentage of
saturated oxygen in the blood, see in figure 2.4.
Figure 2.4: O2 Sat as function of
R IR
[2]
It can be seen in figure 2.4 that 85% of saturated oxygen in yields R IR
= 1 . This is usually the normalization that is done when the
instrument is manufactured according to experimentally results. The error that arises from this method is due to the fact that there are other types of Hemoglobin that we did not consider in our initial assumption. In addition, skin color and nail polish affect the instrument accuracy. In a research that was carried out on 2000 patients [3], 66 of 76 cases of OLI were detected using pulse oximeter (87% of detection) but the drawback of this method is high detection time which may lead to a critical situation when OLI in found.
12
Six patients were under 70% of oxygen saturation, and 30 patients were under 90% of oxygen saturation, and these are situations that are considered as critical situations. The detection time was 2-5 minutes from the OLI occurrence. 3. Capnography [2] The Capnograph measures the CO2 level in the air discharged from the lungs. Most of the Capnography instruments are based on the principle of absorbing IR waves. In a research that was done [3], from 54 cases of OLI only 2 cases were detected (4% of detection) using Capnograph. Therefore, the Capnograph is not a useful method to detect OLI.
It can be summarized that the existing instruments for OLI detection that are used today require experience and attention, their percentage of detection is low. In addition, there is no special instrument to detect OLI that assigned especially to this problem.
13
3.
Summary of Previous Research on OLI Detection by Lungs Sounds Analysis
3.1
Introduction Analysis of respiratory sound signals, in application to diagnosis and
monitoring of pulmonological malfunctioning, is well known in the literature [4][7],[35]-[36]. In order to perform OLI detection, a previous research has suggested methods that were based on recording the patient’s breathing sounds by microphones attached to the patient’s chest. A simple analysis of comparing the amplitude of the recorded sounds in the right and left sides has already been proposed but it was not reliable enough. The reason for this is that each one of the microphones records sounds generated by both of the lungs, and therefore even during OLI situation all of the microphones record the breathings. This chapter summarizes a previous research, which was intended to overcome this problem by suggesting methods to detect OLI from the recorded breathing sound signals.
3.2 Detection
of
OLI
by
LPC
and
PARCOR
Features An automated system for classification between 3 cases of intubation by breathing sounds has proposed and tested on dogs by [2]. The 3 cases of intubation were: 1. Tracheal intubation (TRI). 2. Right lung intubation (RLI). 3. Left lung intubation (LLI).
14
Two microphones were used, to record the sounds at four locations on the dog’s chest. In addition, a strain gauge was used to monitor inspiration/expiration/pause periods. About 30 seconds were recorded in each location. The recorded signals were pre-processed in order to cancel the PCG signal and a classification by extracting LPC [8] and PARCOR [9] features was then performed. In addition, classification by energies ratio was also tested. The results of the classification by LPC and PARACOR features extraction are presented by confusion matrices in Tables 3.1 and 3.2 respectively. Classification by normalized energies ratio was tested and its results are shown in Table 3.3.
Table 3.1: Confusion matrix of classification by LPC features [2] Type of intubation
RLI
LLI
TRI
RLI
68%
10%
22%
LLI
21%
43%
36%
TRI
20%
12%
68%
Table 3.2: Confusion matrix of classification by PARACOR features [2] Type of intubation
RLI
LLI
TRI
RLI
57%
14%
29%
LLI
17%
51%
32%
TRI
14%
20%
66%
15
Table 3.3: Confusion matrix of classification by normalized energies ratio [2] Type of intubation
RLI
LLI
TRI
RLI
89%
1%
10%
LLI
9%
90%
1%
TRI
6%
24%
71%
As it can be seen from the last table, detection of OLI using normalized energy ratio test gave good results. However, only 5 dogs were used in the database hence, the results are not reliable enough and no conclusion about human being can be derived.
3.3
Detection of OLI Using Source Separation In order to overcome the problem described in section 3.1 (The left
microphone records the breathing of the right lung and vice-versa), a source separation algorithm was proposed and tested in [7]. Source separation algorithms mostly assume a certain model of mixing between the generated signals. Estimation of the unknown parameters of the mixing model is done by the recorded signals in order to recover the original generated signals. Block diagram of a basic source separation system is shown in the following figure:
x1 xˆ1 y M = x ˆ = M → A → B → x xK xˆ K Figure 3.1: Source separation block diagram [7]
16
where x = { x1 , x2 ,..., xK } represents the vector of independent sources, A is a mixing T
matrix, and y = { y1 , y2 ,..., yL } is the measurement vector, constructed by y = Ax . T
The goal of the source separation algorithm is to find a conversion to estimate x from the observation vector, y. Assuming that this conversion exists, it can be expressed as the matrix B and the goal of the source separation algorithm can be defined as finding
ˆ -1 Ax = xˆ ≈ x (in cases where K ≠ L a separation matrix B , where By = BAx = A pseudoinverse of the matrix A is needed). The source separation algorithm used in [7] was JADE (Joint Approximate Diagonalization of Eigen Matrices) [10], [11]. This algorithm assumes that the elements of the matrix A are constant, and the source separation is obtained by the matrix B. In order to find the elements of the matrix B fourth order statistics (cumulants) and a contrast function that measures the diagonalization of the cumulants
matrix
are
used.
Cijkl [ y ] = Cum { yi , y j , yk , yl }
φJADE [ y ] =
∑
ijkl ≠ ijkk
and
The the
cumulant
matrix
is
defined
as:
contrast
function
is
defined
as:
C 2ijkl [ y ] . The matrix B, which minimizes the contrast function is
the selected matrix for separation. Minimization of the contrast function consists of whitening operation, cumulants computing and choosing sub-sequence of them, and joint diagonalization of the selected sequence. The joint diagonalization is performed using Jacobi method, which achieves maximal diagonalization of a set of matrices by performing rotation of the set, one axis each stage. The algorithm was tested on synthetic data, which simulated the generated lungs sounds, and showed 100% of OLI detection, under SNR of 3dB. However, implementing the algorithm on real breathing signals recorded in a surgery room,
17
showed only 73% of correct lung ventilation situation. The main reason for these relatively poor results is that the JADE algorithm assumes non-convolutive mixture of the sources, whereas in the reality the source signals are convolved with a filter that represents the chest cavity.
3.4
Motivation for Developing a New OLI Detection Method As explained in Chapter 2, the currently used methods for OLI detection
require experience and attention and their reliability is low. The methods based on lung sounds analysis that were recently published, gave indication that the application of sophisticated signal processing methods may yield an acceptable solution to the OLI detection problem. It was decided to develop a new algorithm, based on information theory, in order to detect the number of ventilated lungs from the recorded breathing sound signals. In addition, because of lack of high quality breathing signals database, a new database of recorded breathing sound signals was established (see appendix A).
18
4.
Blind MIMO Identification and Information Theory as an Approach to Determine the Number of Sources
4.1
Introduction It can be noticed that in order to detect OLI situation, there is no need to
estimate the breathing sound signals themselves, as proposed by [7] and explained in Section 3.3. In addition, there is no importance to determine if the OLI situation is because of right lung intubation (RLI) or left lung intubation (LLI), as proposed by [2] and explained in Section 3.2. Therefore, this work proposes a method for estimation exclusively the number of ventilated lungs The method assumes a MIMO (Multiple-Input-Multiple-Output) system, which relates the lungs and the sensors. In addition, an information theory approach is used to determine the number of sources (lungs) of the system. An overview of existing blind MIMO identification methods is given in this chapter, such as the use of information theory as a model selection criteria developed by Akaike (AIC) [12], and Schwartz and Rissanen (MDL) [14], [15].
4.2 Blind MIMO Identification Based on High Order Statistics The goal of blind MIMO identification is to identify an unknown system H(z) driven by unobservable inputs based on its outputs. However, such a problem has no unique solution if there is no assumption on the unknown system and the
19
unobservable inputs. In most of the research on this subject, the unknown system is assumed to be linear time-invariant (LTI) and restricted to models such as MA and ARMA, and the unobservable inputs are assumed to be white [22]-[24]. The principle behind higher order statistics methods, as already mentioned in Section 3.3, is to use third- or fourth-order statistics called cumulants. For example, the third order cumulant
of
a
stationary
signal,
x(t),
is
defined
as:
C 3 (t1 , t 2 ) = E[ x(t ) x(t + t1 ) x(t + t 2 )] . While an extensive amount of work has been done for estimation of the parameters of a scalar model using higher order statistics of the time series, less attention has been paid to multichannel models. One main drawback of HOS method is that certain independence assumptions between the input signals are required in order to make the use of cumulants a statistical tool to identify H(z). In [25], thirdorder cumulants are employed for estimating multichannel AR model of known orders. Inouye and Sako [26] also considered the identification problem of a linear MIMO system driven by colored inputs using the third-order cumulants of the outputs and clarified the equivalence classes in the blind identification problem. Inouye and Hirano [27] addressed the blind identification problem of linear MIMO systems driven by unobservable colored inputs using HOS, particularly the fourth-order cumulants of the outputs, where the unobservable inputs are mutually, temporally colored linear processes.
20
4.3
Additional
Methods
of
Blind
MIMO
Identification Besides of HOS-based methods of blind MIMO identification, additional methods have been developed. Mostly, those methods are based on assuming certain model and distribution of the sources and noise (usually Gaussian), and estimation of the unknown parameters. Pham and Tong [28] assumed a multichannel AR model without noise, and offered two iterative procedures to find the ML estimator of the unknown parameters of the model, while guaranteeing a stable solution. Hasan and Yahagi [29] proposed a method for the identification of multichannel AR processes contaminated by additive white observation noise, while the noise-compensated multichannel Yule-Walker method is used, applied by two iterative methods. Gorokhov and Loubaton [30] addressed the blind MIMO-FIR identification problem in the case where the number of inputs is strictly less than the number of outputs. Their method is based on the linear prediction approach introduced in [31]-[33], while the unknown parameters are estimated using the weighted least-squares approach. The practical value of the algorithm based on optimal weighting is low, because of an excessive numerical complexity. Hasan, Hossain and Haque [34] presented multichannel autoregressive estimation method using a finite set of noisy observations without a-priori knowledge about the noise power. Their method is based on solving alternatively a set of nonlinear equations using the iterative Newton-Raphson algorithm in order to estimate the noise variance, and a set of Yule-Walker linear equations in order to estimate the AR parameters.
21
It can be noticed that most of the methods of blind MIMO identification developed till now are implemented by iterative solutions, and may require many iterations in order to receive the desired solution. Unlike these methods, a noniterative solution for blind MIMO AR identification problem is developed in this work, lowering the complication time needed to receive the desired solution.
4.4
Akaike Information Criterion – (AIC) The information theoretic criteria for model selection, introduced by Akaike
[12],[13], Schwartz [14], and Rissanen [15] address the following general problem. Given a set of N observations Y = [ y (1),..., y ( N ) ] and a family of models, that is, a parameterized family of probability densities f (Y / θ ) , select the model that best fits the data. The goal of the Akaike Information Criterion (AIC) [12] is to minimize the Kullback-Leibler Divergence (KLD) between the real and estimated probability density functions. It is assumed that there exists a vector parameter vector θ , which is the parameters vector the real probability density function, marked as: g (y ) = f (y | θ ) . A family of probability density functions under the i-th hypothesis is defined by fi (y / θ i , H i ) . The KLD is defined as: g (y ) I ( g , fi ) = ∫ g (y)log dy = f i (y | θ i , H i ) y = ∫ g (y)log ( g (y) ) dy − ∫ g (y)log ( fi (y | θi , Hi ) )dy
y y 144 42444 3 14444244443 S% ( g , g )
S% ( g , f )
In order to minimize I(g,f), S% ( g , f ) has to be maximized. Let us denote:
(4.1)
22 S (θ ,θ ) = S% ( g , g ) = ∫ g (y ) log ( g (y ) ) dy = S% ( f (y | θ ), f (y | θ ) )
(4.2)
y
and
(
S (θ ,θˆ) = S% ( g , f ) = ∫ g (y ) log ( f i (y | θ i , H i ) )dy = S% f (y | θ ), f (y | θˆ) y
)
(4.3)
where θˆ denotes the Maximum-Likelihood estimator of θ . When θˆ is sufficiently close to θ , θˆ = θ + ∆θ , using taylor series, developed up to second order approximation we obtain: I ( g , f ) = S (θ , θ ) − S (θ , θ + ∆ θ ) ≅ ∂S (θ , θ 2 ) ∂ 2 S (θ , θ 2 ) 1 = = S (θ , θ ) − S (θ , θ ) + ∆ ⋅ ∆θ + ∆θ T ⋅ ⋅ θ ∂θ 2 2 ∂θ 2 2 θ 2 =θ θ 2 =θ ∂ 2 S (θ , θ 2 ) 1 T = − ∆θ ⋅ ∂θ 2 2 θ
(4.4)
⋅ ∆θ 2 =θ
where i,j-th element of the matrix
∂S (θ , θ ) ∂S (θ ,θ ) ∂ 2 S (θ , θ ) is defined as: = ∂θ ∂θ i , j ∂θ i ∂θ j
and θ k represents the k-th element of the vector θ . The maximum of S (θ ,θ 2 ) is reached at the point θ 2 = θ . Therefore, that the term ∂S (θ ,θ 2 ) ∂θ 2 θ
in (4.4) vanishes because of the zero value of the derivative in the 2 =θ
maximum point. In addition, we obtain that: ∂ 2 S (θ ,θ 2 ) ∂θ 2 2 θ
2 =θ
= ∫ g (y ) y
∂ 2 log f (y | θ ) ∂ 2 log f (y | θ ) y = d E ∂θ 2 ∂θ 2 θ =θ θ =θ
= − J (θ ) (4.5)
where J (θ ) is the Fisher information matrix. As a result we obtain: I (g, f ) =
(
1 1 ∆θ T J (θ )∆θ = θˆ − θ 2 2
) J(θ )(θˆ − θ ) T
(4.6)
23
The mathematical development of the criterion that is given below is based on [12]. Because of the fact that θˆ is ML estimator, the asymptotic distribution of N (θˆ − θ ) is approximated by a Gaussian distribution with zero-mean and covariance matrix J −1 .
(
N θˆ − θ
) J (θ )(θˆ − θ ) T
Therefore, the distribution of the quadratic form:
for sufficiently large N is approximated, under certain
regularity conditions, as chi-square distribution. The degree of freedom equals to the dimension of the restricted parameter space [16]. Thus, it holds that:
(
E∞ ( 2 NI ( g , f ) ) = N θˆ − θ
)
T
(
)
J (θ ) θˆ − θ + p
(4.7)
where E ∞ denotes the mean of the approximated distribution and p is the dimension of θ , or the number of parameters independently adjusted for the maximization of the likelihood. Therefore, if N N 2 ∑ log f (y (i ) | θ ) − ∑ log f (y (i ) | θˆ) i =1 i =1
(
is used to estimate of N θˆ − θ
(4.8)
) J (θ )(θˆ − θ ) it needs a correction for the downward T
bias introduced by replacing θ by θˆ . This correction is simply realized by adding p to (4.8). It can be noted that the term in (4.8) that includes θ can be ignored because it is constant and does not depend of the selected hypothesis. As a result, the AIC criterion is defined: AIC ( H i ) = −2 log f (y | θˆi ; H i ) + 2 pi
(4.9)
where pi is the number of independent unknown parameters in θˆ under hypothesis Hi .
24
In the proposed model, the unknown number of sources is expressed by the number of independent unknown parameters, pi . The number of sources that leads to minimization of AIC ( H i ) is the estimated number of sources of the system.
4.5
Minimum Description Length – (MDL) Inspired by Akaike’s pioneering work, Schwartz and Rissanen approached the
problem from quite a different points of view. Schwartz’s approach is based on Bayesian arguments, in which each competing model can be assigned a prior probability, and proposed to select the model that yields the maximum posterior probability. Rissanen’s approach is based on information theoretic arguments. Since each model can be used to encode the observed data, Rissanen proposed to select the model that yields the minimum code length. It turns out that in the large sample limit, both Schwartz’s and Rissanen’s approaches yield the same criterion, given by: MDL( H i ) = − log f (y | θˆi ; H i ) + 12 pi log N
(4.10)
where N is the number of samples, and pi is the number of independently adjusted parameters to get θˆ under hypothesis Hi. In similar to AIC, the selected number of sources that leads to minimization of MDL( H i ) is the estimated number of sources of the system. The performances of AIC and MDL strongly depend on the characteristics of the data (e.g., sample size, experimental design, type of random noise) and the models themselves (e.g., model equation, parameters). However, it has been shown [37],[38] that AIC is not consistent and tends to overestimate the number of unknown parameters, while the MDL criterion is consistent but tends to underestimation at low SNRs. These behaviors of the criteria will be examined in simulations further in this work.
25
5. 5.1
The Proposed Model for the Problem
Introduction Recovering the original generated breathing sounds from the recorded
breathing sounds is not essential in order to detect OLI. Therefore, we propose estimating the number of ventilated lungs from the recorded breathing sounds, which is a simpler task. In order to determine the number of ventilated lungs from analysis of the recorded sounds, an AR (Auto-Regressive) model that relates the lungs’ sounds and the microphones’ signals was assumed. This model was chosen because it is commonly used in applications of speech and audio processing and its computational complexity is relatively simple. In addition, a high order AR model can be a good approximation to the more general ARMA (Auto-Regressive Moving Average) model [17]. This Chapter formulates the proposed model for the problem and gives a mathematical development of the Maximum-Likelihood for the unknown parameters of the model.
5.2
Model Formulation Our goal is to estimate the number of sources from the received signals by the
sensors. Fig. 5.1 shows a block diagram of the proposed MIMO AR model, in which
x[n] represents the sources (lungs) of the system, and y[n] represents the sensors (microphones).
26
e[n]
AR Model
x[n]
C
K Channels
( L× K )
∑
L Channels
L Channels
y[n]
Sources’
Microphones’
Signals
Signals
y ( M ) [ n]
A
( L × LM )
delay
LM Channels
Figure 5.1: Block diagram of the MIMO AR model
Let K and L denote the number of sources (lungs) and sensors (microphones), respectively (K
(5.1)
The L ×1 measurements vector is defined as: y[ n] = [ y1[ n]
y2 [ n] ⋅ ⋅ ⋅
T y L [ n]] .
(5.2)
The relation between the source signals and the measurements is assumed to be given by a MIMO AR model:
y[ n] = Ay ( M ) [ n] + Cx[ n] + e[ n],
(5.3)
27
where y ( M ) [n] is an ML × 1 vector defined as follows:
y ( M ) [n] = y1( M ) [n] y (2M ) [n] ⋅ ⋅ ⋅ y (LM ) [n] T
T
T
T
(5.4)
and y (i M ) [n] is a M × 1 vector which contains the past values of the i-th sensor, yi[n], up to sample M:
y (i M ) [n] = [ yi [n − 1] yi [n − 2] ⋅ ⋅ ⋅ yi [n − M ]] . T
(5.5)
A is an L × ML AR coefficeints matrix defined as: T a11 ⋅ A= ⋅ ⋅ aT L1
⋅
⋅
⋅
⋅
⋅
⋅
a1TL ⋅ ⋅ ⋅ aTLL
,
(5.6)
where aij is an M × 1 vector, which relates the samples of the i-th sensor, yi[n], with the past values of the j-th sensor, yj[n-1],…, yj[n-M]. C is an L × K matrix whose i,jth element relates the samples of source j with sensor i. Finally, e[n] is an L × 1 vector representing the additive white Gaussian noise. It is assumed that the noise and source 2 signals are independent, zero-mean, Gaussian with covariance matrices σ I and I,
respectively. The last assumption is employed with no loss of generality, because the covariance of the sources is determined by the matrix C, as it can clearly be seen from (5.3). As a result, we obtain that the conditional distribution of y[n] | y ( M ) [n] is Gaussian: y[n]|y ( ) [n] ~ N ( Ay ( ) [n], R ), where R is defined as: M
M
R = CC T + σ 2 I .
(5.7)
It is noted that the unknown parameters: A, R, σ 2 , M and K must be estimated from a set of N measurements, y[1],…,y[N]. It is also assumed that all the initial conditions are zero, i.e. e[n],x[n]=0 for n<0.
28
5.3
The Maximum-Likelihood Estimator of A and R In order to determine the number of sources, K, we need first to estimate the
unknown matrices, A and R, from the N samples of the data: y[1],…,y[N]. For this purpose, the Maximum–Likelihood (ML) estimator is used. The ML estimator of the matrices A and R, is obtained by maximizing the logarithm of the conditional probability density function (pdf) of the output samples given the unknown matrices, which is (see Appendix B for proof): log f ( y[1],..., y[ N ] | R, A ) = −
NL N log(2π ) − log R − 2 2
1 N − ∑ (y[n] − Ay ( M ) [n])T R −1 (y[n] − Ay ( M ) [n]) . 2 n =1
(5.8)
Maximization of (5.8) with respect to the unknown parameters, A,R, is achieved via equating the corresponding partial derivatives to zero. Only the last term of (5.8) which is:
1 N − ∑ (y[n] − Ay( M ) [n])T R−1 (y[n] − Ay( M ) [n]) = 2 n=1 ( M )T (M ) −1 −1 T T −1 ( M ) T 1 N y [n]R y[n] − y [n]A R y [n] − y [n]R Ay [n] + = − ∑ 2 n=1 (Ay( M ) [n])T R−1 (Ay( M ) [n])
(5.9)
is relevant to calculate the derivative of the log-likelihood with respect to A. The derivative of a scalar a with respect to a matrix B is defined as: ∂∂ba11 ⋅ ∂a = ⋅ ∂B ⋅ ∂a ∂bM 1
⋅
⋅
⋅
⋅
⋅
⋅
b11 ⋅ ⋅ ⋅ , where B = ⋅ ⋅ ⋅ ∂a bM 1 ∂bMN ∂a ∂b1 N
⋅
⋅
⋅
⋅
⋅
⋅
b1N ⋅ ⋅ . ⋅ bMN
29
We shall use the following identities (for proof of non-trivial identities see appendix B): ∂ ∂ xT Ay = y T AT x = xy T , ∂A ∂A
[
(5.10)
]
∂ ( Ax) T C( Ax) = (C + C T )( Ax)x T . ∂A
(5.11)
If A is a square matrix and f(A) is a scalar function, then:
∂ log f ( A) 1 ∂f ( A) = ∂A f ( A) ∂A
.
(5.12)
If A is a square and invertible matrix, then: ∂A A
= A ( A −1 ) T
(5.13)
,
and we obtain:
∂ log( f (y[1],..., y[ N ] | R, A) = ∂A (5.14) 1 N −1 ( M )T ( M )T (M ) ( M )T −T −1 −T [n] − R y[n]y [n] + (R + R ) Ay [n]y [n], − ∑ −R y[n]y 2 n =1
∂ log( f (y[1],..., y[ N ] | R, A) = ∂R −1 (M ) ( M )T T T [n]AT N 1 1 N y[n]y [n] − Ay [n]y [n] − y[n]y −1 T = R−1 R (R ) − ∑ 2 2 n =1 + Ay ( M ) [n]y ( M )T [n]AT
.
(5.15)
In order to find the ML estimator of A and R, the above derivatives should be equated to zero. Since R is a covariance matrix then, R = R T and R −1 = R −T . Therefore we obtain two matrix equations with two unknown matrix variables, A and R −1 : −
1 N −1 −R y[n]y ( M )T [n] − R −1y[n]y ( M )T [n] + 2R −1Ay ( M ) [n]y ( M )T [n] = 0 ∑ 2 n =1 (M ) ( M )T T T [n]AT + N 1 N y[n]y [ n] − Ay [ n]y [ n] − y[ n]y R − ∑ =0 2 2 n =1 Ay ( M ) [n]y ( M )T [n]AT
(5.16)
(5.17)
30
Eq. (5.17) can be simplified to: N N N ˆ y ( M ) [n]yT [n] − y[n]y ( M )T [n] A ˆT + ˆ = ∑ y[n]y T [n] − A NR ∑ ∑ n =1 n =1 n =1
(5.18)
ˆ y ( M ) [n]y ( M )T [n] A ˆT . +A ∑ n =1 N
Extraction of A from (5.16) gives:
ˆ A ML
N N = ∑ y[n]y ( M )T [n] ∑ y ( M ) [n]y ( M )T [n] n =1 n =1
Note that (5.19) is derived while assuming the matrix
−1
(5.19)
.
N
∑ y(
M)
[n]y ( M )T [n] is
n =1
invertible.After few straightforward manipulations we obtain: ˆ = 1 R ML N
N
∑ y[n]yT [n] − Aˆ ML n =1
1 N
N
∑ y(
M)
(5.20)
[n]y T [n] .
n =1
The last equation may be expressed as: ˆ = 1 R ML N
∑ ( y[n] − Aˆ N
n =1
ML
)(
)
T
ˆ y ( M ) [ n] . y ( M ) [n] y[n] − A ML
(5.21)
The ML estimator is known to be a consistent estimator. Therefore, if the real AR MIMO system is stable, the developed ML estimated system would also be stable, given a sufficient number of samples. It can be shown that the developed ML estimator of the matrix A (5.19) is actually an extension of the SISO LPC estimation to the MIMO case. After substituting L=1 and K=1 in the MIMO AR model, i.e. in Eq. (5.19), and performing transpose operation on both sides of the equation, (5.19) turns to:
31
a1 y[n − 1] a y[n − 2] 2 ⋅ N ⋅ = • y n − y n − ⋅ ⋅ ⋅ y n − M [ 1] [ 2] [ ] [ ] ∑ ⋅ ⋅ n =1 ⋅ ⋅ aM y[n − M ]
−1
y[n − 1] y[n − 2] N ⋅ ∑ y[n] • ⋅ n =1 ⋅ y[n − M ]
N
by defining: ryy ( j ) = ∑ y[n] y[n − j ] and assuming g stationary signals, one obtains: n =1
ryy (1) a1 ryy (0) a r (1) ryy (0) 2 yy ⋅ ⋅ ⋅ = ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ aM ryy ( M − 1) ryy ( M − 2)
Multiplying
both
sides
ryy (1) ryy (0) r (1) ryy (0) yy ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ryy ( M − 1) ryy ( M − 2)
of
⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅
the
⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅
last
ryy ( M − 1) ryy ( M − 2) ⋅ ⋅ ⋅ ryy (0)
equation
−1
ryy (1) r (2) yy ⋅ ⋅ ⋅ ryy ( M )
with
the
matrix
ryy ( M − 1) ryy ( M − 2) ⋅ , leads to the well-known Yule⋅ ⋅ ryy (0)
Walker equations, which are the basis to SISO LPC solution [40].
5.4
Implementing MDL and AIC Algorithms for Estimating the Number of Signals and the Model order In order to estimate the number of sources, K, and the AR model order, M,
information theoretic criteria are applied in this Section, based on Wax and Kailath [18]. It can clearly be seen from (5.7) that the rank of the matrix CCT equals to the
32
number of source signals. Using a well-known spectral representation theorem from linear algebra, R can be expressed as: K
R = ∑ (λi − σ 2 )v i v iT + σ 2 I
(5.22)
i=1
where λ1 ,..., λK and v 1 ,..., v K are the highest eigenvalues and the corresponding eigenvectors, of R, respectively. (For proof see Appendix B). Note that K ∈ {0,1,..., L − 1} ranges over the set of all possible number of signals. Denoting by
Θ( K , M ) the parameter vector of the model, it follows that:
(
T ΘT ( K , M ) = λ1 ,..., λK , σ 2 , v1T ,..., v1T , a11 ,..., aTLL
)
(5.23)
T where a11 ,..., aTLL are L vectors with dimensions M × L each, composing the matrix
A, which consists of L2 M elements. ˆ , and p ,..., p , denote its Let the l1 ≥ l2 ≥ ⋅ ⋅ ⋅ ≥ lL denote the eigenvalues of R 1 L corresponding
eigenvectors.
According
to
[19],
the
ML
estimators
of
λ1 ,..., λK , σ 2 , v1T ,..., vTK are as follows: λˆi = li
i=1,…,K
(5.24a)
vˆ i = pi
i=1,…,K
(5.24b)
σˆ 2 =
1 L−K
L
∑l
i = K +1
i
(5.24c)
After substituting the ML estimators (5.24) in the log-likelihood function (5.8) and omitting terms that do not depend on the parameter vector Θ ( K , M ) with some straightforward manipulations, we obtain: L− K L K N 1 ˆ log f (y[1],..., y[ N ] | Θ ( K , M )) = const + log ∏ li L − K ∑ li i =1 2 i = K +1
(5.25)
33
(See Appendix B for derivation of (5.25)). The number of free adjustment parameters in Θ( K , M ) is obtained by counting the number of degrees of freedom of the space spanned by Θ( K , M ) . It follows that
Θ( K , M ) has
K {
eigenvalues
+ 1{ + σ
2
LK {
from eigenvectors
+
L2 M {
parameters. However, not all of
from matrix A
the parameters are independently adjusted; the eigenvectors are constrained to have a unit norm, a constraint that reduces K degrees of freedom. In addition, the eigenvectors also have to be orthogonal to each other; this constraint reduces 1 K ( K − 1) degrees of freedom. Therefore, the number of free adjustment parameters 2 of the model is: 1 1 K L + − K + L2 M + 1 . 2 2
(5.26)
The last term of (5.26) can be omitted because it is independent of K and M. The AIC for this problem is therefore given by: L−K L K 1 AIC ( K , M ) = N log ∏ li L − K ∑ li + 2 K ( L + 12 − 12 K ) + L2 M i =1 i = K +1
(5.27).
Similarly, the MDL criterion turns out to be: L− K L K N 1 1 MDL( K , M ) = log ∏ li L − K ∑ li + K ( L + 12 − 12 K ) + L2 M log N (5.28), i =1 2 2 i = K +1
ˆ It can be noted that M is expressed by the eigenvalues of R which is a function of ML ˆ (see (5.21)), a matrix with dimensions of M × ML . The number of signals and A ML
the AR model order are determined by the pair K and M for which either the AIC or the MDL is minimized.
34
L
It can be noted that the term
1 L− K
∑l
i = K +1
i
which appears in (5.27) and (5.28) is
actually the estimation of the noise level, σ 2 (according to (5.24c)). Therefore, if we discuss the case where the noise level is known, there will be two changes in the development of the two criteria: 1. Parameter vector of the model will not contain σ 2 , and (5.23) turns into:
(
T ΘT ( K , M ) = λ1 ,..., λK , v1T ,..., v1T , a11 ,..., aTLL
)
.
(5.29)
Therefore, the number of free adjustment parameters of the model is now (instead of (5.26)): 1 1 K L + − K + L2 M . 2 2
(5.30)
Note that this change from the case where the noise level is known actually does not affect the penalty function of the criteria. L 2. The term L −1K ∑ li i = K +1
L−K
now turns into: σ 2( L − K ) .
Therefore, the criteria for the case where σ 2 is known are given by: K AIC ( K , M ) = N log ∏ li σ 2( L − K ) + 2 K ( L + 12 − 12 K ) + L2 M , i =1
(
MDL( K , M ) =
)
(5.31)
K 1 N log ∏ li σ 2( L − K ) + K ( L + 12 − 12 K ) + L2 M log N , (5.32) 2 i =1 2
(
)
The pair K and M for which either the AIC or the MDL is minimized determines the number of signals and the AR order.
35
6. 6.1
Simulation Results
Introduction In order to evaluate the performance of the estimators as a function of the
parameters of the model, simulations with synthesis data were performed. The simulations checked the performance of the ML estimators of the matrices A and R, the estimated number of sources, K, and the AR model order, M. In addition, a comparison between the simulation results of the different criteria (MDL and AIC) was carried out. This Chapter ends with conclusions that arise from the simulation results.
6.2
Results for Implementing the ML Estimator The proposed MIMO AR system as defined in (5.3) was simulated and ML
estimators of R and A (Eqs. 5.22 and 5.23) were calculated. In order to quantify the performance of estimator, Normalized Mean Square Error (NMSE) for estimation of the matrices was defined as follows: L
NMSE ( A) =
LM
∑∑ (aij −aˆij )2 i =1 j =1 L LM
∑∑ a i =1 j =1
ij
2
L
, NMSE (R ) =
L
∑∑ (r i =1 j =1 L L
ij
−rˆij ) 2
∑∑ r i =1 j =1
, 2
ij
where aij and rij represent the i,j-th element of the matrices A and R, respectively. In all of the simulations the parameters of the system were selected to be (unless otherwise is indicated): K=2 (Number of sources) L=4 (Number of sensors)
36
M=5 (AR order)
σ 2 = 1 (The noise variance) 0.3 0.15 0.15 0.3 . The matrix C was chosen to be: C = 0.9 0.6 0.6 0.9
In Fig. 6.1 NMSE of the estimated matrices was drawn as function of the number of samples (N) that ranges between 100 to 10000. The NMSE was calculated according to averaging over 100 iterations for each N. NMSE of estimating A and R versus the number of samples
1
10
NMSE of estimated A NMSE of estimated R 0
Normalized Mean Square Error
10
-1
10
-2
10
-3
10
-4
10
2
10
3
10 N - number of samples
4
10
Figure 6.1: NMSE of ML estimators as a function of N A trend of decrease in the NMSE as the number of samples increases can clearly be seen in Fig. 6.1, as expected to be from ML estimator.
37
6.3
Results for Implementing the AIC and the MDL Criterion The behavior of the AIC and MDL criteria as a function of the number of
independent samples, N, was examined. Each one of the criteria was calculated for number of sources, K, in the range between 1 to 3 and model order, M, in the range between 1 to 7. Each N was tested for a total number of J=100 iterations, and probability of error of K or M for each N is computed as: Pe =
number of uncorrectly estimated K,M J
(6.1)
The next two figures depict the probability of error in estimating K and M, respectively, while the number of independent samples, N, ranges between 100 to 10000.
38 Pe of estimating number of sources, K 1 Pe of AIC Pe of MDL
0.9 0.8 0.7
Pe
0.6 0.5 0.4 0.3 0.2 0.1 0 2 10
3
4
10 Number of samles - N
10
Figure 6.2: Pe of estimating K as a function of N Pe of estimating AR order model, M 1 Pe of AIC Pe of MDL
0.9 0.8 0.7
Pe
0.6 0.5 0.4 0.3 0.2 0.1 0 2 10
3
10 Number of samles - N
Figure 6.3: Pe of estimating M as a function of N
4
10
39
In order to explain the necessity of adding the penalty function to the loglikelihood function (see (4.9), (4.10), (5.27), (5.28)), the components of AIC and MDL, which are the log-likelihood and the penalty function, were drawn as a function of tested model order, M. In addition, the sum of these values, which compose the corresponding criterion, was also drawn. The number of independent samples, N, was set to 1000, and the number of sources, K, was set to 2. It can clearly be seen that the penalty function allows both criteria to correctly estimate the model order (M=5 - in the minimum point of each criterion). 4
2500 3.5
The MDL AICcriterion
x 10
Log-likelihood function penalty function the MDL criterion AIC
3 2000 2.5 1500 2 1.5 1000 1 500 0.5 0
10
2
2 4
6 3 8 10 4 12 5 14 M - the AR tested model
16 6
18
20 7
Figure 6.4: Log-likelihood, AIC penalty function and their sum as a function of M
40
MDL
1400
Log-likelihood function penalty function MDL
1200 1000 800 600 400 200 0
1
2
3
4 5 M - the AR tested model
6
7
Figure 6.5: Log-likelihood, MDL penalty function and their sum as a function of M
Additional important performance measure is the histogram of the decisions of the AIC and the MDL under low number of samples conditions. Therefore, we set N=300 and the examined the criteria decisions for 100 iteration, while the iteration are different from each other by the random noise. As a results, the histogram of the decisions of the criteria for K and M, is drawn in the next two figures:
41 Histrogram of estimated number of sources
% of estimated number of sources by the criterion
100
AIC MDL
90 80 70 60 50 40 30 20 10 0
1
2 3 Tested number of sources - K
Figure 6.6: Distribution of MDL and AIC decisions on K, N=300. (True K=2) Histrogram of estimated AR model order
100
AIC MDL
% of estimated AR model by the criterion
90 80 70 60 50 40 30 20 10 0
1
2
3
4 5 6 7 8 Tested AR order model - M
9
10
Figure 6.7: Distribution of MDL and AIC decisions on M, N=300. (True K=1)
42
6.4 Simulation Results Under Non-Point Sources Assumption As explained in Chapter 2, each lung is composed of several independent point sources, and therefore should be treated as a scattered sources. In order to examine the algorithm performances under this assumption, scattered and independent sources were synthesized as follows. The matrix C now represents two scattered sources, in which each source is scattered into 4 different and independent sources. Therefore, the columns of C were chosen to be: C = [c1 c1 + ε∆c11 c1 + 2ε∆c12
c1 + 3ε∆c13
c2
c 2 + 2ε∆c 21 c 2 + 3ε∆c 22
where c1 is orthogonal to c 2 , and the vectors
{ci
∆ci1
∆c i 2
c 2 + ε∆c 23 ]
∆ci 3 } is an
orthonormal group. ε is a constant which actually represents the scattering level. The vector of sources x[n] represents 8 independent sources. Therefore, the product Cx[n] represents 2 scattered sources, where each one is constructed by 4 independent sources, and the scattering level is controlled by ε . The ML estimator was used to estimate the matrices A and R under the assumption of scattered sources and the matrix C as defined above. The NMSE of the estimated matrices, A and R, was examined as a function of the scattering level, ε , where the number of independent samples, N, was constant and equal to 1000. In Fig. 6.8 the results of this simulation is shown, and it can clearly be noticed that the NMSE of the estimated matrices is not affected by the scattering level, therefore, even under the assumption of scattered sources, the ML estimator can be used as an estimator for the unknown matrices.
43 NMSE of estimating A and R , N=1000
0
10
Normalized Mean Square Error
NMSE of estimated A NMSE of estimated R
-1
10
-2
10
-3
10
0
0.1
0.2
0.3
0.4
0.5 ε
0.6
0.7
0.8
0.9
1
Figure 6.8: NMSE of the estimated matrices under scattered sources assumption as a function of scattering parameters.
In addition, under the same conditions, the AIC and MDL criteria were implemented, in order to estimate the number of sources, K, as a function of the scattering level, ε . In Figs. 6.9 and 6.10 the probability of error, Pe, defined in (6.1), is drawn as a function of ε for K=1 and K=2, respectively. In addition, the eigenvalues of the estimated matrix R are also drawn as a function of ε for K=1 and K=2 in Figs. 6.11 and 6.12, respectively.
44 Pe of estimating number of sources, K 1 Pe of AIC Pe of MDL
0.9 0.8 0.7
Pe
0.6 0.5 0.4 0.3 0.2 0.1 0
0
0.1
0.2
0.3
0.4
0.5 ε
0.6
0.7
0.8
0.9
1
Figure 6.9: Pe of estimating K as a function of scattering level, ε for K=1
Pe of estimating number of sources, K 1 Pe of AIC Pe of MDL
0.9 0.8 0.7
Pe
0.6 0.5 0.4 0.3 0.2 0.1 0
0
0.1
0.2
0.3
0.4
0.5 ε
0.6
0.7
0.8
0.9
1
Figure 6.10: Pe of estimating K as a function of scattering level, ε for K=2.
45 The eigenvalues of R, K=1
12 λ1 λ2 λ3 λ4
10
8
6
4
2
0
0
0.1
0.2
0.3
0.4
0.5 ε
0.6
0.7
0.8
0.9
1
Figure 6.11: The eigenvalues of R versus the scattering level, ε , for K=1 The eigenvalues of R, K=2
18 λ1 λ2 λ3 λ4
16 14 12 10 8 6 4 2 0
0
0.1
0.2
0.3
0.4
0.5 ε
0.6
0.7
0.8
0.9
1
Figure 6.12: The eigenvalues of R versus the scattering level, ε , for K=2.
46
The explanation for decreasing in the probability of error for K=2 when ε is in the range between 0.5 to 0.7
(Fig. 6.10) can clearly be found by observing the
eigenvalues for this case. It can be seen from Fig. 6.12 that in this case, the lowest eigenvalues ( λ1 and λ2 ) are getting closer to each other, and therefore both criteria assums that they are representing the noise level and decide about the correct number of sources, K=2.
6.5
Simulation Conclusions From observation of the simulations in the previous section, several
conclusions may be obtained: 1. It can be noted that the NMSE of the ML estimator decreases as the number of samples, N, grows (Fig. 6.1). This fact shows that the developed ML estimators for the matrices A and R are consistent, as expected from the ML estimator to be. 2. From Figs. 6.2-6.3 it can be noted that given sufficient number of samples, both criteria, MDL and AIC, successfully estimate the number of sources, K, and the AR model order, M. 3. Generally, it can be said (according to Figs. 6.2 and 6.3) that for low number of samples the estimation of the AR model order, M, is better than the estimation of the number of sources, K. This phenomenon probably occurs because M is the order of a modeled system, an AR MIMO model, in contrary to K, which is estimated according to a non-modeled matrix, R. 4. From Figs. 6.6-6.7 it is noted that when the AIC wrongly estimates K or M, it tends to over-estimate it, while the MDL tends to sub-estimation. This
47
phenomenon fits the theoretical behavior of the estimator, as published in [38]. 5. In Fig. 6.8 it is seen that even when the sources are scattered, the performance of the ML estimator of the matrices A and R is not highly affected, even in high scattering level. Therefore, the ML estimator developed for point sources can be used also in cases where the sources are scattered. 6. In Figs. 6.9-6.10 it is seen that AIC and MDL do not accurately estimate the number of sources when the sources are scattered above a certain level ( ε > 0.2 ). According to the MDL implementation for the model (5.28) it can be seen that the criterion estimates the noise variance from the lower eigenvalues. Therefore, when the eigenvalues are distant from each other, the criterion decides that they are product of source signal and not noise. Therefore, these methods cannot be used in order to estimate the number of sources in the selected model, if the sources are scattered (as in the case of lungs sources). 7. It can however be seen from Figs. 6.11-6.12 that even when the sources are highly scattered, the eigenvalues of the source signals sub-space are clearly separated from the eigenvalues of the noise sub-space.
48
7.
Mismatches Between the Assumed and Real Models
7.1
Introduction The use of AIC and MDL in the previous chapters is restricted to cases where
the sources are point sources. It is obvious that the point sources assumption is not satisfied in the real world because the anatomy of the human body is not parallel to the theoretical assumed model. A discussion about the mismatches between the assumed and the real model is presented in this chapter, and modification of the algorithm in order to adopt it for the “real world” problem is offered and tested.
7.2
Description
of
the
Mismatches
Between
Assumed and Real Models It was expected that the proposed model for the problem would be different from the real model. First of all, the lungs are not point sources as they are treated in the selected model. In fact, as explained in Section 2.1, each lung is constructed of alveoli and this is the reason why each lung is composed of many distributed and independent point sources. Therefore, instead of the results in the theoretic and simulated model, where the simulated energy of two lungs was expressed over two eigenvalues of R, in the “real world” the energy of the point sources composing the lungs is distributed over all of the eigenvalues of R as seen in Figs. 6.11-6.12. As a result, it turns out that in the real world problem the implemented model order selection criteria, MDL and AIC, are unable to estimate the real number of ventilated lungs, similarly to the simulation results for non-point sources (Figs. 6.9-6.10).
49
Indeed, modeling the human body by AR model may not be accurate. Still, high AR model order can sometimes be a good approximation even for some nonlinear models. Hence, it is believed that this mismatch is not the main reason for the fact that the criteria are unable to estimate the real number of sources. A modification of the solution according to the real anatomy of a body (scattered sources) had to be done, as explained below.
7.3
Modification of the Solution for the Real Model It can be noted that the implementation of AIC and MDL to MIMO-AR
model, as appears in Section 5.4, was developed for general case (with the constraint that the number of sources is smaller than the number of sensors). As explained above, the lungs are non-point sources and therefore AIC and MDL were unable to correctly estimate number of sources when applied to real breathing sound signals. In addition, in the OLI special case, the real number of sources can only be one or two. Therefore, a modification of the solution for the OLI detection problem is done as follows. We shall examine the likelihood ratio, which is: f ( y[1],..., y[ N ] | R, A; H1 ) N ˆ | H − N log R ˆ |H log = log R 2 1 y y R A f [1],..., [ N ] | , ; H 2 2 ( 2 )
(
(
ˆ |H where R 1
)
(
ˆ |H and R 2
)
)
(
)
(7.1)
represents the determinant of the estimated covariance
matrix, R, under OLI and TRI situations, respectively. ((7.1) can clearly be obtained from the proofs of (5.8) and (5.25) in Appendix B). Using identity (B.13) and assuming the noise variance, σ 2 , is known, and under point sources assumption, (7.1) can be simplified into:
50
( ) ( )
L−2
l1l2 σ 2 f ( y[1],..., y[ N ] | R, A; H1 ) N l N log = log 22 = log L −1 2 σ 2 l1 σ f ( y[1],..., y[ N ] | R, A; H 2 ) 2
(7.2)
2 ˆ ( l ≥ l ), and L is the number of where {li }i =1 are the two highest eigenvalues of R 1 2
ˆ is sensors. As can clearly be seen from (7.2), the second highest eigenvalue of R actually the detector of OLI situation, under the assumption that σ 2 is known, and that the sources are point sources. Indeed, in the OLI case the lungs are not point sources, and therefore the energy of the sources is distributed over all the eigenvalues ˆ and is not restricted to the eigenvalues that represent the source signals. of R ˆ as a detector for Even though, selecting the second highest eigenvalue of R OLI situation proved itself as a reliable detector. The reason for this, is that during OLI situation only one lung is ventilated and therefore, the sources are less scattered than in case of TRI where both lungs are ventilated. This fact causes the energy of ˆ . Therefore, in these point sources to be less distributed over the eigenvalues of R ˆ value is smaller than in TRI case, a fact OLI case the second highest eigenvalue of R which enables us to choose it as an OLI detector. Figs. 6.11 and 6.12 are given in Fig. 7.1 side by side in the same scale in order to illustrate this phenomenon.
51
The eigenvalues of R, K =1
16
λ1 λ2 λ3 λ4
14
12
10
10
8
8
6
6
4
4
2
2
0
0.2
λ1 λ2 λ3 λ4
14
12
0
The eigenvalues of R, K =2
16
0.4
0.6 ε
0.8
1
0
0
0.2
0.4
0.6 ε
Figure 7.1: The eigenvalues of R versus the scattering level, ε .
0.8
1
52
8. 8.1
Implementation on Real Data
Introduction After performing the modification of the solution for the real world problem,
implementation on real data, obtained by recording of real patients in a surgery room, was performed. Because of the noisy environment of the surgery room, preprocessing of the recorded sounds was essential prior to implementing the proposed method. In this chapter the experimental setup of the patient’s breathings recording is described, followed by the description of the pre-processing performed on the recorded breathings sounds. The last section of this chapter shows the results of an OLI detector based on the second highest eigenvalue of R, as described in Section 7.3.
8.2
Experimental Setup In order to examine the proposed method and to check the possibility of
developing a monitoring tool in the future, a database of lung sounds during OLI and TRI was established. The database was composed of 24 patients which were recorded in a surgery room in the two described situations: during correct ventilation, when the tip of the tube is placed above the carina (TRI), and during a situation of OLI when the tip of the tube is under the carina and only one lung is ventilated. In the experiments, the microphones were attached to the patients’ back, as shown in Fig. 8.1, recorded the breathing sounds of the patients in both situations. The ventilations were performed manually and not mechanically, in order to achieve higher SNR in the recorded sounds, and the real position of the tube was validated
53
each time by fiber-optic. The experiments were performed in the main surgery room of Soroka medical center – Beer-Sheva Israel, during the anesthesia part at the beginning of the surgery.
Figure 8.1: Experimental set-up – view from back (TRI shown)
8.3
Pre-Processing Because of the noisy environment of the surgery room and in order to fit the
breathing sounds to the model assumptions, pre-processing of the breathing sounds was performed. A Graphical User Interface (GUI) that enables the user to perform pre-processing jobs on the recoded sounds was developed especially for this purpose. This tool enables the user to observe the recorded signals from each microphone, to transfer the signal through a butterworth bandpass filter, to resample and to normalize
54
the signal according to the user’s defined parameters. A screen capture of this GUI is shown in figure 8.2.
Figure 8.2: The developed GUI for pre-processing operations
Using this tool, few pre-processing operations were performed on the recorded signals: 1. In order to attenuate noises which are not in the spectral range of the breathing signals, such as monitor beeps, doctors’ discussions etc., the recorded signals
55
were band-pass filtered by a Butterworth filter with a bandwidth of 100Hz600Hz. 2. The recorded sounds were sampled at 4kHz. Because of the cut-off frequency of 600Hz, down-sampling operation with a factor of 0.3 was performed. 3. The signal amplitude from each microphone depends on the location of the microphone on the patient’s body, the anatomy of the recorded patient and on the gain of the sampling system. Therefore, it turned out that each channel’s output had a different signal amplitude, and the noise variance, σ 2 , was different on each one of the sensors, in contrary to our assumption in the model. In order to overcome this problem, normalization of each channel according to the noise level on the channel was performed as follows: First, the noise variance of the i-th channel
(denoted as σ i2 ) was estimated
according to time segments that contain only noise. Next, the i-th channel was multiplied by the factor:
1000
σ i2
. Performing this operation over all of the 4
channels, resulted an identical noise variance to all channels: σ 2 = 1000 .
8.4
Real Data Results The recorded breathing signals contain both situations of OLI and TRI. Our
goal in implementing the algorithm on the real data is to detect the OLI cases from the recorded signals. In order to do so, the pre-processed data were divided into windows of 2000 samples each, with 80% of overlapping. Trying to use information criteria in order to estimate the AR model resulted an infinite model order when applied on real breathing signals. The reason is, that the real breathing signals contain an additive noise, that turns the AR system into ARMA system (see [39] for
56
additional explanation about this phenomenon). Therefore, an arbitrary AR order of 15 was set. The unknown matrices A and R were estimated for each window, using the ML estimated developed in Section 5.3. The second highest eigenvalue of R was extracted. Fig. 8.3 shows a few breathing cycles of both OLI and TRI situations, recorded by the four microphones after pre-processing. As it can be seen from this figure, determination between OLI and TRI cases by only the amplitude of the recorded sounds is not a simple task, as explained in Section 3.1. Upper left microphone
4
x 10
One Lung
1
Two Lungs
0 -1 5
10
15
20
Time [Sec.]
25
30
35
40
30
35
40
30
35
40
30
35
40
Lower left microphone
4
x 10
Two Lungs
One Lung
2 1 0 -1 -2
5
10
15
20
25
Time [Sec.]
Upper right microphone
4
x 10
Two Lungs
One Lung
5
0 -5 5
10
15
20
Time [Sec.]
25
Lower right microphone
4
x 10
4
Two Lungs
One Lung
2 0 -2 5
10
15
20
25
Time [Sec.]
Figure 8.3: Examples of the recorded signals – four microphones, both OLI and TRI.
57
ˆ as a function of time, as a result Fig. 8.4 shows the second highest eigenvalue of R ML
of processing the experiment shown in Fig. 8.1. As it can clearly be seen from Fig. 8.4, OLI and TRI cases can clearly be discriminated, by the value of the second ˆ highest eigenvalue of R ML in every breathing cycle. The results of the proposed
algorithm were consistent over the 24 experiment performed on patients as explained in Section 8.2. 6
5.5
x 10
Second highest eigenvalue of R
5
Two Lungs
4.5 4 3.5 3 2.5 2
One Lung
1.5 1 0.5 0
5
10
15
20
25
30
35
40
Time [Sec.] Figure 8.4: Second eigenvalue of Rˆ as a function of time.
In order to quantify the results, the values of the second highest eigenvalues were normalized in all of the 24 experiments as follows. First, the variance of second eigenvalue’s values in the regions of TRI was computed (and noted by vl , where l
58
indicated the experiment number). Next, the values of the second eigenvalue of the lth experiment were multiplied by:
108 . Performing this operation over all of the vl
experiment resulted that the TRI regions in the second eigenvalue’s values have a unitary variance of 108 . The histograms of the second highest eigenvalue under TRI/OLI situations are shown in Fig. 8.5. The left bars represent the histogram of the normalized values when OLI situation occurs, and the right bars represents the values in TRI situations.
Figure 8.5: The histrogram of the second eigenvalue under TRI/OLI situations.
Because of the fact that the database is not comprehensive enough to be used for both training and validation of the proposed method, estimation of the
59
performance of the system was performed as follows. Twenty different experiments were used to extract the histograms in order to train the system. The rest of the 4 experiments were used to validate the system and were tested according to the extracted statistics in the training process. This process was repeated 6 times, each time a different group of four validation experiment was used. As a result, a validation was performed using a total number of 24 experiments, in a “patient independent” mode. There are two types of errors in OLI detection: Pmiss – the probability of a true OLI to be wrongly detected as TRI, and false alarm, PFA, the probability of TRI to be detected as OLI. The Detection Error Tradeoff (DET) curve is a common mean to display these errors. The DET curve provides information about the device’s performance, where each point on the curve shows the PFA and Pmiss for a given threshold. The threshold of a real monitoring system should be calculated according to the requested sensitivity of the system, while taking into consideration the allowed Pmiss of the system. Fig. 8.6 shows the DET curve of the proposed decision system, which was computed according to the 6 iterations that were described above. The Equal Error Rate (EER) point is defined as the point on the DET curve where Pmiss= PFA, is 4.8%. In this work’s application, more importance should be given to Pmiss rather than to PFA. Therefore, it is assumed that in a practical system the selected activity point on the DET curve will be where Pmiss=2% and PFA=9%.
60
System Performance
0
10
Pmiss
DET curve EER
-1
10
-2
10
-2
10
-1
10 PFA
0
10
Figure 8.6: The DET of a detector based on the second highest eigenvalue of the estimated R.
61
9.
Summary and Conclusions
Summary This thesis includes both theoretical and practical achievements. In the theoretical aspect, MIMO AR system identification was performed according to second order statistics. In addition, the number of sources and the model order of the system were also estimated. In fact, this is the first time that information criteria are used to estimate the combination of the number of sources and the model order in a convolutive mixture. Theoretical development of implementation of the AIC and the MDL criteria for combined estimation the number of sources and the model order in MIMO AR model is presented in Chapter 5. Simulation results of this theoretical algorithm are shown in Chapter 6. From the practical point of view, this thesis endeavored to take OLI detection research one-step forward by presenting the concept of a new automatic system for OLI detection. In Chapter 7 the results from the theoretical part are modified according to the anatomy of the human body, and algorithm for detection of OLI by monitoring lungs sound is presented. In Chapter 8 it is shown that assuming a MIMO AR model with distributed sources and selecting the second highest eigenvalue of the residual covariance matrix as a feature, proves itself as a reliable method for detection of OLI for the database used in this research.
Performance In order to examine the algorithm performance, a database of recorded breathing sound signals of patients during OLI and TRI situations was established. To our knowledge, this is the first time that such a comprehensive database of real OLI
62
breathing signals was established. After implementing the developed algorithm on the recorded breathing signals of 24 patients, the results were summarized and shown in Figs. 8.5 and 8.6. From these figures it can be seen that using the offered method an OLI correct detection probability of 95.2% with Pmiss of 4.8% and PFA of 4.8% can be achieved. Note that there may not be enough patients recorded in the database so it may not demonstrate the real probability density function of the proposed classifier (the second highest eigenvalue of the residual covariance matrix) under OLI and TRI situations. I is expected that the more patients will be added to the database, the less concave the DET curve will be, indicating inferior performance of the detector.
Recommendation for future work This work proves the concept of detection of OLI situation using breathing sound signals analysis. The proposed method is meant to be the basis for developing a real-time monitoring system for OLI detection. As explained in Section 8.2 the database was based on experiments in which the ventilations were performed manually. It is supposed that implementing the method into a real-time machine in order to detect OLI while the patient is mechanically ventilated would provide even better results. The reason for this being is that while the patient’s lungs are mechanically ventilated, the variance of the amplitude of the breathing signals is much smaller than it is when using manual ventilation (as performed in our experiments). Because of the fact that pre-processing operations according to the surgery conditions has to be done, it is assumed that automatic training of the system before every surgery in order to enable it to set the optimal gain for each microphone will be essential part of the future system. In addition, our estimators may need to be
63
modified in order to have computation complexity that enables them to be computed in a real-time system. From the theoretical point of view, the MIMO AR identification developed in this work estimates the AR coefficients, the matrix A, and the covariance matrix, R. It can be noted that the mixing matrix, C, is not estimated because it is not an essential need in order to estimate the number of sources. An interesting idea is to estimate the mixing matrix, C, in order to develop a new source separation method. In addition, in this work, the stability of the solution is not guaranteed. Stabilizing the solution or offering a new method that ensures a stable solution, may lead to a better performance.
64
References [1] http://www.howstuffworks.com [2] G. Sod-Moriah, “Ventilation Monitoring During Anesthesia and Respiratory Intensive Care – M.Sc. Thesis Report,” Ben-Gurion University of the Negev,
Department of Electrical and Computer Engineering, Beer-Sheva, Israel, October 1995. [3] R. K. Webb, J. H. van der Walt, W. B. Runciman, J. A. Williamson et. al. , “Which Monitor? An Analysis of 2000 Incident Reports,” Anesthesia and Intensive
Care, vol. 21, pp. 529-542, 1993. [4] A. Cohen and D. Landsberg, “Analysis and Automatic Classification of Breathing Sounds,” IEEE Trans. vol. BME-31, pp. 585-590, 1984. [5] A. Cohen and A. Berstein, “Acoustic Transmission of the Respiratory System using Speech Stimulation,” IEEE Trans., vol. BME-38, pp. 126-132, 1991. [6] G. Sod-Moriah, A. Cohen and G. Gurman, “Detection of One-lung Intubation Incidents in General Anesthesia and Intensive Care,” Proc. of the 13th Int. Conf.
BIOSIGNAL 96, pp. 282-284, Brno, Czech Republic, 1996. [7] A. David, T. Lazmi, “Detection of One-Lung Intubation – Final Project Report,”
Ben-Gurion University of the Negev, Department of Electrical and Computer Engineering, Beer-Sheva, Israel, 2003. [8] J. Makhoul, “Linear Prediction: A Tutorial Review,” Proc. of the IEEE, vol. 63 no. 4, pp. 561-580, April 1975. [9] A. Cohen, “Biomedical Signal Processing,” CRC Press, Florida 1986. [10] J. F. Cardoso, “Blind Signal Separation: Statistical Principles,” Proc. of the
IEEE, vol 9, no. 10, pp. 2009-2025, 1998.
65
[11] J.F. Cardoso, A. Souloumiac, “Blind Beamforming for Non-Gaussian Signals,”
Proc. Inst. Elect. Eng., vol. 140, no. 6, pp. 362-370, 1993. [12] H. Akiake, “A New Look at the Statistical Model Identification,” IEEE Trans.
Automat. Contr. vol. AC-19, pp. 716-723, 1974. [13] H. Akiake, “Information Theory and an Extension of the Maximum Likelihood Princilple,” Proc. 2nd Int. Symp. Inform. Theory, Suppl. Problems of Control an
Inform. Theory, pp. 267-281, 1973 [14] G. Schwartz, “Estimating the Dimension of a Model,” Ann. Stat., Vol. 6, pp. 461-464, 1978. [15] J. Rissanen, “Modeling by Shortest Data Description,” Automatica, vol. 14, pp. 465-471, 1978. [16] P.J. Huber, “The Behavior of Maximum Likelihood Estimates Under Nonstandard Conditions,” Proc. 5th Berkeley Symp. Mathematical Statistics and
Probability, vol. 1, pp. 221-233, 1967. [17] B. Porat, “Digital Processing of Random Signals, Theory and Methods,” Prentice Hall, Englewood Cliffs, NJ, 1994. [18] M. Wax and T. Kailath, “Detection of Signals by Information Theoretic Criteria,” IEEE Trans. Acoust. Speech, Signal Processing, vol. ASSP-33, pp. 387392, 1985. [19] T. W. Anderson, “Asymptotic Theory for Principal Component Analysis,” Ann.
J. Math. Stat., vol. 34, pp. 122-148, 1963. [20] http://www.oximeter.org [21] G. Kowalewski, “Einf¨urung in die Determinantentheorie,” pp. 91-93, Chelsea, New York, 1948.
66
[22] G. B. Giannakis, Y. Inouye, and J. M. Mendel, “Cummulant Based Identification of Multichannel Moving-Average Processes,” IEEE Trans. Automat. Contr., vol. 34, pp.783-787, 1989. [23] Y. Inouye and T. Umeda, “Estimation of Multivariate ARMA Processes Using Cumulants," IEICE Trans. on Fundamentals, Communications and Computer
Sciences, vol. E77-A, no. 5, pp. 748-759, 1994. [24] A. Swami, G. B. Giannakis, and S. Shamsunder, “A Unified Approach to Modeling Multichannel ARMA Processes Using Cummulants,” IEEE Trans. Signal
Processing, vol. 42, pp. 898-913, 1994. [25] M. R. Raghuveer, “Multichannel Bispectrum Estimation,” in Proc. 3rd ASSP
Wrokshop on Spectrum Estimation and Modeling, pp. 21-24, 1986. [26] Y. Inouye and B. Sako, “Identifiability of Multichannel Linear Systems Driven by Colored Inputs and its Application to Blind Signal Separation,” in Proc. ISCAS-
94, vol.5, pp. 57-60, 1994. [27] Y. Inouye, and K. Hirano, “Cumulant-Based Blind Identification of Linear Multi-Input-Multi-Output Systems Driven by Colored Inputs,” IEEE Trans. Signal
Processing vol. 45, no. 6, pp. 1543-1552, 1997. [28] D. T. Pham and D.Q. Tong, “Maximum Likelihood Estimation for a Multivariate Autoregressive Model,” IEEE Trans. Signal Processing, vol. 42, no. 11, pp. 30613072, 1994. [29] M. K. Hasan, and T. Yahagi, “An Iterative Method for Identification of Multichannel Autoregressive Processes with Additive Observation Noise,” IEICE
Trans. Fund., vol. E79-A, no. 5, pp. 674-680, 1996.
67
[30] A. Gorokhov and P. Loubaton, “Blind Identification of MIMO-FIR Aystems: A Generalized Linear Prediction Approach,” Signal Processing, vol. 73, pp. 105-124, 1999. [31] K. Abed-Meraim, P. Duhamel, D. Gesbert, P. Loubaton, S. Mayrargue, E. Moulines, D. Slock, “Prediction Error Methods for Time-Domain Blind Identification of Multichannel FIR Filters,” In Proc. 1995 IEEE ICASSP, Detroit, MI, pp. 1968– 1971, 1995.
[32] K. Abed-Meraim, P. Loubaton, E. Moulines, “A Subspace Algorithm for Certain Blind Identification Problems,” IEEE Trans. Inform. Theory, vol. 43, pp. 499–511, 1997. [33] K. Abed-Meraim, E. Moulines, P. Loubaton, “Prediction Error Methods for Second-Order Blind Identification,” IEEE Trans. Signal Processing, vol. 45, pp. 694– 705, 1997. [34] M. K. Hasan, M. J. Hossain and M.A. Haque, “Parameter Estimation of Multichannel Autoregressive Processes in Noise,” Signal Processing, vol. 83, pp. 603-610, 2003. [35] G. R. Wodicka, K. N. Stevens, H. L. Golub, E. G. Cravalho and D. C. Shannon “A Model of Acoustic Transmission in the Respiratory System,” IEEE Trans. on
Biomedical Engineering, vol. 36, no. 9, pp. 925-934, 1989. [36] V. K. Iyer, P. A. Rammoorthy and Y. Ploysongsang, “Autoregressive Modeling of Lung Sounds: Characterization of Source and Transmission,” IEEE Tras. on
Biomedical Engineering, vol. 36, no. 11, pp. 1133-1137, 1989. [37] K. M. Wong, Q. T. Zhang, J. P. Reilly and P. C. Yip, “On Information Theoretic Criteria for Determining the Number of Signals in High Resolution Array Processing,” IEEE Trans. on Signal Processing, vol. 98, no. 11, pp. 1959-1971, 1990.
68
[38] A. P. Liavas and P. A. Regalia, “On the Behavior of Information Theoretic Criteria for Model Order Selection,” IEEE Trans. on Signal Processing, vol. 49, no. 8, pp. 1689-1695, 2001. [38] J. G. Proakis, and D. G. Manolakis, Digital Signal Processing: Principles,
Algorithms and Applications, Macmillan (1996). [40] E. C. Ifeachor and B. W. Jervis, Digital Signal Processing, A Practical
Approach, pp. 710-714, Prentice Hall, 2002.
69
Appendix A: Database Description Appendix A.1 Introduction In order to test the proposed algorithm, a suitable database had to be collected or acquired. The database consists of recording the breathing sounds of 24 patients who were operated in the main surgery room at Soroka University Medical Center, Beer-Sheva, Israel. Detailed description of the database is given in the following sections.
Appendix A.2 The Sampling System The sampling system consisted of three part: 4 microphones, amplifier and an antialiasing filter. The sampling was performed by a Laptop computer which was equipped with A/D converter connected as shown in the following figure:
Microphones
Laptop Computer
Amplifier
A/D Converter
Figure A.1: Schematic diagram of the sampling system
Because of the safety reasons, the system has to floating. Hence, the laptop researchable battery was used as the power source. Fig. A.12 depicts of the sampling system in the surgery room:
70
Figure A.2: Picture of the sampling system
The microphones are piezo-electric microphones with internal Low-Pass Filter which suits the spectral range of the breathing sound signals.
The amplifier has a controlled gain for each channel. The laptop computer is equipped with Keithley A/D converter, 12bit resolution, and is able to receive up to 8 single-ended channels. The computer is also equipped with network adapter, which enables the user to transmit the sampled data directly via the Internet.
Appendix A.3 The Data Acquisition The sampling process was performed in the main surgery room of the hospital. After few experiments in which the microphones were located on the chest
71
of the patient, it was decided to move the location of the microphones to the back of the patient, because of several reasons: •
When the microphones are located on the back of the patient, the patient (in most cases) actually lies on them and his body blocks most of the surrounding voices, such as doctors discussions, monitor sounds etc.
•
The sound transmission from the lungs to the back is better than it is to the chest.
•
The chest cavity of the patient is expanded during ventilation, and as a result, the chest of the patient rises up and down during ventilation. When microphones are attached to the chest of patient, the movement of the chest vibrates the microphones and this vibration results in noises in the recorded signals. Familiar phenomenon does not occur if the microphones are attached to the patient’s back.
And indeed, the quality of the recorded breathing sounds was highly improved after the microphones were transferred to the back of the patients. Each phase of recording began the moment the ventilation tube inserted to the trachea. In order to obtain recording samples of OLI and TRI situations, the patient was intubated using the “tube withdrawal” method. Using this method, the tube is inserted deeply until OLI situation occurs and recorded. Afterwards, the tube is drawn slowly, until it reaches its correct location, and the TRI situation is recorded. The location of the tip of the tube was validated and written after each movement of it, by stethoscope and the final location validated by fiber optic. Every phase of recording lasts 2-4 minutes. A total of 24 patients comprised the study. All of the patients were informed and gave their consent to participate in the study. The Ethics
72
Committee of the Ministry of Health (the Helsinki Committee) approved the study. A picture of a patient during the experiment is shown in Fig. A.3.
Microphones’ cables Figure A.3: A patient’s breathings sounds are recorded during an experiment
In Fig. A.13 the microphones’ cables are seen in both sides of the patient. They lead to the microphones attached to the back of patient and records the breathing, while the anesthesiologist validates location of the tip of the tube by fiber-optic.
73
Appendix B: Proof of Mathematical Identities Proof of 5.8 Claim: Given the model defined in (5.1)-(5.7) the pdf of the output samples given the unknown matrices is: log f ( y[1],..., y[ N ] | R , A ) = − −
NL N log(2π ) − log R − 2 2
1 N (y[ n] − Ay ( M ) [n])T R −1 (y[ n] − Ay ( M ) [n]) ∑ 2 n =1
Proof: According to conditional probability definition, the following equation is obtained:
f (y[1],..., y[ N ] | R, A) = f (y[ N ] | y[1],..., y[ N − 1]; R, A) ⋅ f (y[1],..., y[ N − 1] | R, A) (B.1) It is can clearly be seen from (5.3)-(5.5) that y[ N ] | y[ N − M ],..., y[ N − 1] is independent of y[1],…,y[N-M-1]. Therefore, using the notation in (5.2), (B.1) is simplified to: f (y[1],..., y[ N ] | R, A ) = f (y[ N ] | y[ N − M ],..., y[ N − 1]; R , A) ⋅ ⋅ f (y[1],..., y[ N − 1] | R, A) = = f (y[ N ] | y
(M )
(B.2)
[ N ]; R, A ) ⋅ f (y[1],..., y[ N − 1] | R , A)
Repeating this procedure N times leads to: f (y[1],..., y[ N ] | R, A) = f (y[ N ] | y ( M ) [ N ]; R , A) ⋅ f (y[ N − 1] | y ( M ) [ N − 1]; R, A) ⋅ ... ⋅ f (y[1] | y ( M ) [1]; R , A) ⋅ f (y[0] | R, A)
(B.3)
y[0]=0 with probability 1 (zero initial conditions). Therefore, (B.3) becomes: f (y[1], ..., y[ N ] | R, A) = f (y[ N ] | y ( M ) [ N ]; R , A) ⋅ f (y[ N − 1] | y ( M ) [ N − 1]; R, A) ⋅ ... ⋅ f (y[1] | y ( M ) [1]; R , A)
(B.4)
Substituting Gaussian pdf into (B.4) leads to: N
f (y[1],..., y[ N ] | R, A ) = ∏ i =1
1
( 2π )
L 2
R
1 2
e
(
)
(
)
T 1 (M ) [ n ] R −1 y [ n ]− Ay ( M ) [ n ] − 2 y [ n ]− Ay
Performing log on both sides of (B.5) leads to (5.8).
(B.5)
74
Proof of 5.10 Claim: Given a matrix A with dimensions M × N and vectors x and y with dimensions M × 1 and N × 1 respectively which are independent with the elements of
A then: ∂ xT Ay = ∂ y T AT x = xy T ∂A ∂A Proof: We begin with proof of the right side equality: ∂f ∂a 11 ⋅ ∂ ∂ T ∂ [f ]= ⋅ x Ay = ∑∑ xi y j aij = ∂A i j ∂A ∂A 142 4 43 4 ⋅ f ∂f ∂a M1
[
x1 y1 ⋅ = ⋅ ⋅ x M y1
⋅
⋅
⋅
⋅
⋅
⋅
]
⋅
⋅
⋅
⋅
⋅
⋅
∂f ∂a1N ⋅
⋅ = ⋅ ∂f ∂a MN
x1 y N ⋅ ⋅ = xy T ⋅ x M y N
Applying the trace operation on xT Ay (which is a scalar) leads to proof left side equality by identical way.
75
Proof of 5.11 a11 ⋅ Claim: Given a matrix A defined as: A = ⋅ ⋅ a1L
as:
⋅
c11 ⋅ C= ⋅ ⋅ c1L
⋅
⋅
⋅
⋅
c1L ⋅ ⋅ ⋅ c LLL
⋅
x = [x1 ⋅ ⋅ ⋅ x N ] then: T
and
a
[
⋅
⋅
⋅
⋅
⋅
⋅
column
a1N ⋅ ⋅ , a matrix C defined ⋅ a LN
vector
x
defined
as:
]
∂ ( Ax) T C( Ax) = (C + C T )( Ax) x T ∂A
Proof: It can be clearly seen that:
(Ax )T
N = ∑ a1 j x j j =1
L
⋅ ⋅ ⋅
∑a j =1
N
(Ax )T C = ∑∑ c1i aij x j i=1
N
Lj
L
N
∑∑ c
⋅ ⋅ ⋅
j =1
x j , and therefore,
i =1 j =1
N
j =1
k =1
(Ax )T CAx = ∑ ∑ cli ∑ aij x j ∑ alk xk , L
l =1
L
i =1
a x j and as a results we obtain that:
Li ij
N
which is actually a scalar. Derivation of
this scalar with respect to the matrix A is a matrix, in which every element of it is a derivation of the scalar with respect to the relevant element in the A, as explained in section 5.3. The resulted scalar can be decomposed into:
N
j =1
k =1
(Ax )T CAx = ∑ ∑ cli ∑ aij x j ∑ alk xk = L
l =1
L
i=1
N
L N N N N L ∑ cli ∑ aij x j ∑ alk xk + ∑ ck1i ∑ aij x j ∑ ak1k xk ∑ j =1 l =1 i =1 k =1 i =1 j =1 k =1 4 314243 1444424 443 l ≠ k1 144244 B B B 12 2 4441142444443 14 L
B1
76
Let us derivate this scalar with respect to ak1k 2 , the (k1, k2) element of A. We will derive the element B1 first. B12 is constant with respect to ak1k 2 , therefore, we obtain: L N ∂B1 = ∑ ∑ alk xk clk1 xk 2 . ∂ak1k 2 l =1 k =1 l ≠ k1
B2 can be decomposed into:
N L N B2 = ∑ a k1k x k + a k1k 2 x x2 ∑ c k1i ∑ aij x j + aik 2 x k 2 = kk =≠1k i =1 jj =≠1k 2 2 N L N = ∑ a k1k x k + a k1k 2 x x2 ∑ c k1i ∑ aij x j + c k1k1 a k1k 2 x k2 j =1 kk =≠1k ii =≠1k 2 1 After performing the multiplications B2 will be composed of 6 elements. Derivation of each element with respect to ak1k 2 will give us:
L N N N ∂B2 2 = c k1k1 x k2 ∑ a k1k x k + x k 2 ∑ c k1i ∑ aij x j + c k1k1 ∑ a k1 j x j + 2a k1k2 c k1k 2 x k 2 ∂a k1k 2 k =1 i =1 j =1 j =1 k ≠k2 j ≠ k2 i ≠ k1 And as a result we obtain that: ∂B1 ∂B1 ∂ (Ax ) CAx = + = ∂ak1k 2 ∂ak1k 2 ∂ak1k 2 T
L N N N N 2 = ∑ ∑ alk xk clk1 xk 2 + ck1k1 xk 2 ∑ ak1k xk + xk 2 ∑ ck1i ∑ aij x j + ck1k1 ∑ ak1 j x j + 2ak1k 2 ck1k 2 xk 2 i =1 l =1 k =1 k =1 j =1 j =1 i ≠ k1 l ≠ k1 k ≠ k2 j ≠k2 14444442444444 3 L
B3
B3 can be written as: B3 = ∑ ∑ alk x k c k1l x k2 + c k1k1 x k2 ∑ a k1k x k , and as a result we l =1 k =1 k =1 L
l ≠ k1
obtain that:
N
N
k ≠ k2
77
T L N ∂ (Ax ) CAx N 2 == ∑ c lk 1 + c k1l ∑ a lk x k x k 2 + 2c k1k1 x k 2 ∑ a k1k x k + 2a k1k 2 c k1k 2 x k 2 = ∂a k1k 2 l =1 k =1 k =1
(
)
l ≠ k1
k ≠k2
(B.6)
L N L N N = ∑ c lk 1 + c k1l ∑ a lk x k x k 2 + 2c k1k1 x k 2 ∑ a k1k x k = ∑ c lk1 + c k1l ∑ a lk x k x k 2 l =1 k =1 l =1 k =1 k =1 l ≠k
(
)
(
)
1
Now we will show that the element (k1,k2) in the expression: (C + C T )( Ax)x T is identical to the element (B.6) above:
(C + C ) T
2c11 c + c 21 12 ⋅ = ⋅ ⋅ c L1 + c1L
N x1 ∑ a1i x i iN=1 x a 2i xi 1∑ i =1 ⋅ Axx T = ⋅ ⋅ N x a Li x i ∑ i =1
c12 + c 21
⋅
⋅
⋅
2c 22
⋅
⋅
⋅
⋅
⋅
⋅
⋅
N
x 2 ∑ a 1i x i
c L1 + c1L c L 2 + c 2 L ⋅ , and ⋅ ⋅ 2c LL
⋅
⋅
⋅
⋅
⋅
⋅
i =1 N
x 2 ∑ a 2i xi i =1
⋅
N x L ∑ a 1i x i i =1 N x L ∑ a 2i xi i =1 ⋅ . ⋅ ⋅ N x L ∑ a Li x i i =1
The (k1,k2) in the expression: (C + C T )( Ax)x T is the product of the k1-th row in
(C + C ) with the k2–th column in the matrix Axx T
[c
k11
+ c1k1
c 2 k1 + c k1 2
⋅ ⋅ ⋅ c Lk1 + c k1L
]
T
which is:
N x k2 ∑ a1i xi iN=1 x a 2i xi L k2 ∑ N i =1 = c + c a x x ⋅ • ∑ lk k k2 lk1 k1l ∑ l =1 k =1 ⋅ ⋅ N x k ∑ a Li xi 2 i =1
(
)
78
We got an identical expression to (B.6) which is the (k1,k2) element of the derivation of (Ax ) CAx with respect to A, therefore: T
[
]
∂ ( Ax) T C( Ax) = (C + C T )( Ax) x T ∂A
79
Proof of 5.13 Claim: Given a square and invertible matrix A with dimensions NxN then: ∂A ∂A
= A ( A −1 ) T
∂A ∂a11 ⋅ ∂A Proof: It is obvious that: = ⋅ ∂A ⋅ ∂A ∂a N 1
⋅
⋅
⋅
⋅
⋅
⋅
∂A ∂a1N ⋅ ⋅ . ⋅ ∂A ∂a11
We will define kij an the algebraic complementary of element aij in the matrix A so that: N
N
j =1
i =1
A = ∑ aij k ij = ∑ a ij k ij
(B.7)
The adjoint matrix is defined as: ⋅
k11 ⋅ adj ( A) = K T = ⋅ ⋅ k1N
From (B.7) it is obtained that
∂A ∂aij
=
⋅
⋅
⋅
⋅
⋅
k N1 ⋅ ⋅ ⋅ k NN
(B.8)
∂ ∂aij
N ∑ aij k ij = k ij , and using (B.8) it is j =1
⋅
⋅
obtained that: k11 ⋅ ∂A = ⋅ ∂A ⋅ k N 1
⋅
⋅
⋅
⋅
k1N ⋅ T ⋅ = ( adj ( A) ) = K ⋅ k NN
(B.9)
80
According to [21], for every square matrix: adj ( A) = A A −1
From (B.9) and (B.10) it is obtained that: ∂A ∂A
(
= A ( A −1 )
)
T
= A ( A −1 ) T
(B.10).
81
Proof of 5.22 Definitions: C is matrix with degree of K. T λ1 ,..., λK are the eigenvalues of CC in descending order ( λ1 ≥ λ2 ≥ ⋅⋅⋅ ≥ λK ).
R is a covariance matrix defined as: R = CC T + σ 2 I .
l1 ,..., l L are the eigenvalues of R in descending order ( l1 ≥ l2 ≥ ⋅ ⋅ ⋅ ≥ lL ).
v1 ,..., v L are the eigenvectors of R, where v i corresponds to li . Claim: R can be expressed as: K
R = ∑ (λi − σ 2 )v i v iT + σ 2 I i =1
Proof:
Because of the fact that R is a covariance matrix, it can be decomposed as R = VΛV T , where Λ is a diagonal matrix of the eigenvalues of R, V is a matrix
consists of { v i }L as columns and VV T = I . i =1 Based on [18] it follows that:
li = λi ,
i=1,…,K
li = σ 2 ,
i=K+1,…,L
Therefore, R can be developed as follows: L
L
R = VΛVT − σ 2I + σ 2I = ∑ li vi vTi − σ 2I + σ 2I = ∑ li vi vTi − σ 2 V VT + σ 2I = i =1
L
L
∑ l v v − ∑σ i =1
i
K
i
T i
i =1
2
vi viT + σ 2I = ∑ (li − σ 2 ) vi vTi + σ 2I =
= ∑ (λi − σ 2 ) vi vTi + i =1
i =1
L
i =1
L
∑
i = K +1
K
(σ 2 − σ 2 ) vi vTi + σ 2I = ∑ (λi − σ 2 )vi vTi + σ 2I i =1
82
Proof of 5.25 Claim: Substituting (5.24) in (5.8) and omitting terms independent of Θ( K , M ) as defined in (5.23), leads to (5.25). Proof: Substituting (5.24) in (5.8) leads equation (B.11):
log f ( y[1],..., y[ N ] | R, A ) = −
NL N ˆ − log(2π ) − log R 2 2
1 N ˆ [n])T R ˆ [n]) ˆ −1 (y[n] − Au − ∑ (y[n] − Au 2 n =1
(B.11)
The sum term in (B.11) is a scalar, and therefore the trace operation can be performed on it:
N ˆ [n])T R ˆ [n])} = ˆ −1 ( y[n] − Au tr ∑ {(y[n] − Au n =1 N ˆ −1 ˆ [n])( y[n] − Au ˆ [n])T } = tr NR ˆ −1R ˆ = tr ( NI ) = NL tr ∑ {R (y[n] − Au n =1
(
)
(B.12)
The second equality in (B.12) arises by substituting (5.21). The determinant of matrix is the product of its eigenvalues. Therefore, recalling
{li }i=1 L
the eigenvalues of
ˆ in descending order ( l ≥ l ≥ ⋅ ⋅ ⋅ ≥ l ), R ˆ can be R 1 2 L
simplified into: L
ˆ = R ∏ li
(B.13).
i =1
ˆ are used to estimate As explained in section 5.4, the smallest L-K eigenvalues of R
σ 2 , therefore, (B.13) turns into: K 1 L ˆ = R l ∏ i L − K ∑ li i = K +1 i =1
L−K
(B.14)
83
Therefore, substituting (B.12), (B.14) and Θ( K , M ) into (B.11) leads to:
log f (y[1],..., y[ N ] | Θ( K , M )) = L−K L K NL N NL 1 = log(2π ) + log ∏ li L − K ∑ li + i =1 { 2 243 2 i = K +1 14 2 (*)
(B.15)
(*)
The elements marked with (*) in (B.15) are independent of K and M, and omitting them from (B.15) leads to (5.25).
