116
edits
Changes
From OpenEMT.org
no edit summary
==Demodulation implementation==
**Demodulation implementation**
=======
The asynchronous demodulation process takes place in Matlab using an efficient matrix calculation technique. This method involves recording a number of samples of each signal of interest. This has the advantage of reduced solving time since all samples are available at the time of calculation. This is in contrast to real-time processing where samples are gathered and processed one-by-one in turn.
The signals of interest we wish to sample are:
1. The induced voltage on the 5-DOF sensor coil(s)
2. The composite coil current signal
Initially we consider the induced voltage on a single sensor that has been discretised through sampling:
<math display="block">x[n] = \sum_{i=1}^N V_i\sin\left(\frac{2 \pi f_i n}{f_s}+\varphi_i\right)</math>
where N is the number of frequency of interest (N = 8, one frequency per emitter coil), <math>f_i$ is the frequency of interest in hertz, <math>f_s$ is the sampling frequency of the DAQ, <math>V_i</math> is the amplitude and <math>\varphi_i</math> is an associated phase shift. Collecting and string <math>p$ points results in a vector <math>\mathbf{X}</math> containing <math>p$ samples:
<math display="block">\mathbf{X}=\left[x[0],x[1] ... x[p-1]\right]</math>
In the previous section asynchronous demodulation was discussed. The first step of the demodulation involves multiplying the sampled signal by two sinusoids <math>Y(n)=\sin(2\pi f_{ri} n)</math> and <math>X(n)=\cos(2\pi f_{ri} n)</math>, where <math>f_{ri}</math> is the simulated reference signal close in frequency to <math>f_i</math>. Using Euler's relation we combine <math>X(n)$ and <math>Y(n)</math> into a single discrete complex sinusoid:
<math display="block">\cos(2\pi f_{ri}n) + j\sin(2\pi f_{ri} n) = e^{\frac{2 \pi f_{ri} j n }{f_s}}</math>
where <math>j = \sqrt{-1}</math>. This complex exponential encapsulates the two quadrature signal components required for demodulation. Given the <math>N$ emitter coil frequencies, <math>N$ complex exponentials are required for the demodulation of each frequency. A <math>p \times N$ matrix of complex exponentials is created such that when pre-multiplied by <math>\mathbf{X}</math> results in a vector of voltage magnitude:
<math display="block"> \mathbf{E} = \begin{bmatrix}
\epsilon_i[0] & \dots & \epsilon_N[0] \newline
\vdots& \ddots & \vdots \newline
\epsilon_i[p-1] & \dots & \epsilon_N[p-1]
\end{bmatrix}$$
where
<math display="block">\epsilon_i[n]=e^{\frac{2 \pi f_{ri} j n }{f_s}}$$
The row vector of demodulated voltages is given by premultiplying <math>\mathbf{E}</math> by the row vector of <math>p$ samples:
<math display="block">\mathbf{Y} = [\widetilde{V}_1 ... \widetilde{V}_N] = 2\mathbf{X}\mathbf{E}</math>
<math>\mathbf{Y}</math> is a complex valued row-vector due to the presence of <math>j</math> in the complex exponential calculation. The absolute value of each entry in <math>\mathbf{Y}</math> is the amplitude received frequency of interest, <math>V_i</math>.
Finite impulse response (FIR filter)
------------------------------------
A finite impulse response filter is implemented as part of the demodulation step. This low pass FIR filter eliminates unwanted frequency components from the input sample stream <math>mathbf{X}</math>. Consider an FIR filter denoted by <math>f_i$ with $p$ coefficients (\textit{i.e.} equal to the number of input samples). The output of such a filter can be represented as:
<math display="block">p[n] = \sum_{i=0}^{p-1}f_ix[n-i]</math>
where <math>p[n]</math> is the filtered sample vector and <math>f_i</math> is a single filter tap coefficient. Each filter coefficient can be represented by a single matrix <math>F</math> as shown below:
<math display="block">F = [f_0, f_1, f_2 \dots f_{p-1}]</math>
The FIR filter can be applied to the demodulator by scaling each input sample of <math>X$ with the corresponding FIR coefficient given by <math>F$. This is achieved using an element by element multiplication:
<math display="block">\mathbf{Y} = 2(\mathbf{X}\circ \mathbf{F})\mathbf{E} </math>
where <math>\circ</math>represents the element-wise multiplication operator.
The amplitude of each frequency component can be determined by calculating the absolute value of <math>\mathbf{Y}</math>, which is a complex quantity. The phase of <math>\mathbf{Y}</math>an be found using the complex argument of <math>Y$
<math display="block">[V_1 V_2 \dots V_N] = |\mathbf{Y}|</math>
<math display="block">[\varphi_1, \varphi_2 \dots \varphi_N] = \arg(\mathbf{Y})</math>
The implementation of this calculation can be found in the Matlab code provided in the OSF file repository.
