# Real-Time Lightweight Chaotic Encryption for 5G IoT Enabled Lip-Reading
Driven Secure Hearing-Aid

Ahsan Adeel  Department of Computing Science and Mathematics, University of
Stirling, FK9 4LA, UK  Jawad Ahmad  School of Engineering and Built
Environment, Glasgow Caledonian University, Glasgow, G4 0BA, UK
<EMAIL_ADDRESS>Amir Hussain  Department of Computing Science and
Mathematics, University of Stirling, FK9 4LA, UK

###### Abstract

Hearing-loss is the third most common chronic health condition. More than 10
million people in the United Kingdom (UK) and approximately 360 million in the
world are suffering from a debilitating hearing loss, with the number
estimated to rise to 14.5m by 2031. Existing audio-only hearing-aids are known
to perform poorly in noisy situations where overwhelming noise is present.
Next-generation audio-visual (lip-reading driven) hearing-aids stand as a
major enabler to realise more intelligible audio. However, high data rate, low
latency, low computational complexity, and privacy are some of the major
bottlenecks to the successful deployment of such advanced hearing-aids. To
address these challenges, we envision an integration of 5G Cloud-Radio Access
Network (C-RAN), Internet of Things (IoT), and strong privacy algorithms to
fully benefit from the possibilities these technologies have to offer. The
envisioned 5G IoT enabled secure audio-visual (AV) hearing-aid transmits the
encrypted compressed AV information and receives encrypted enhanced
reconstructed speech in real-time which fully addresses cybersecurity attacks
such as location privacy and eavesdropping. For security implementation, a
real-time lightweight AV encryption is utilized. For speech enhancement, the
received AV information in the cloud is used to filter noisy audio using both
deep learning and analytical acoustic modelling (filtering based approach). To
offload the computational complexity and real-time optimization issues, the
framework runs deep learning and big data optimization processes in the
background on the cloud. Specifically, in this work, three key contributions
are reported: (1) 5G IoT enabled secure audio-visual hearing-aid framework
that aims to achieve a round-trip latency up to 5ms with 100 Mbps datarate (2)
Real-time lightweight audio-visual encryption based on piece-wise linear
chaotic map (PWLCM), chebyshev map, secure hash, and a novel substitution box
(S-Box) algorithms (3) Lip-reading driven deep learning approach for speech
enhancement in the cloud. The effectiveness and security of the proposed
secure AV hearing-aid is extensively evaluated using widely known security
metrics such as correlation coefficient, entropy, contrast, energy, number of
pixel change rate (NPCR) and unified average change intensity (UACI). Lip-
reading driven deep learning approach for speech enhancement is evaluated
under four different dynamic real-world commercially-motivated scenarios
(cafe, street junction, public transport, and pedestrian area) using benchmark
Grid and ChiME3 corpora. The critical analysis in terms of both speech
enhancement and AV encryption demonstrate the potential of the envisioned
technology in acquiring high quality speech reconstruction and secure mobile
AV hearing aid communication.

##  1 Introduction

Hearing impairment is a hidden disability with no painful symptoms. People
with serious hearing-loss find themselves socially isolated and depressed with
more negative consequences including headaches, muscle tension, increased
stress, insecurity, and sadness. Hearing aids are the most widely used devices
for the majority of hearing losses compensation. The global hearing aid
industry, estimated around US $6 billion, is expected to grow at 6 percent
annually through 2020, according to the market research firm
MarketsandMarkets. However, existing hearing aids often work poorly for speech
in noise. Current devices make the signal more audible but remains deficient
to restore intelligibility i.e., no improvement in SNR. Thus, the existing
audio-only hearing-aids are not robust to reverberation and intelligibility
wins at the cost of higher cognitive load in noisy environment [1][2].

Despite decades of research, only few speech enhancement algorithms could
reliably increase the intelligibility of speech in noise, especially in
extreme noisy conditions such as cocktail party. A limited number of research
developments in the field of speech enhancement have been implemented into
commercially available hearing-aids. For example, spectral subtraction can be
very effective in stationary conditions, but the processed speech remains
unintelligible. In case of multiple microphones availability, beamforming
algorithms can possibly lead to improvements in speech intelligibility.
However, such approaches are hard to employ in an unpredictable noisy
situations. Consequently, existing audio-only hearing aid approaches achieve
benefit just by simply amplifying signal and offering very little advantage
for speech in high levels of noise [3].

Human performance in noisy environment is known to be dependent upon both
aural and visual cues, which are combined by sophisticated multi-level
integration strategies to improve intelligibility. The multimodal nature of
the speech is well established in literature, and it is well understood how
speech is produced by the vibration of vocal folds and configuration of the
articulatory organs. The correlation between the visible properties of the
articulatory organs (e.g., lips, teeth, tongue) and speech reception has been
previously shown in numerous behavioural studies [4][5][6][7]. Therefore, a
clear visibility of some articulatory organs could be effectively utilized to
extract a clean speech signal out of a noisy audio background. The biggest
advantage of using visual cues to extract clean audio features is their
inherent noise immunity [8].

However, embracing the multimodal nature of speech presents both opportunities
and challenges for hearing assitive technology. The real-time implementation
of AV hearing-aid demands for high datarate, low latency, low computational
complexity, and high security. To address these challenges, we propose a novel
integration of 5G, IoT, and strong privacy algorithms. The proposed AV
hearing-aid framework is envisioned to address challenges such as
cybersecurity attacks (location privacy, eavesdropping), interference between
medical IoT devices (that can cause hearing-aids to operate incorrectly with
potentially life threatening consequences), low-cost wireless technology
design, low power consumption, limited battery, and high datarate requirement.
Inspired by our previous work [9], the novel wireless hearing-aid framework
offloads the computational complexity and real-time optimization issues by
running deep learning and big data optimization processes in the background on
the cloud. The hearing-aid transmits the encrypted compressed audio/visual
information and receives encrypted enhanced reconstructed speech in real-time.
The hearing-aid connects to an indoor 5G wireless access point and
back/fronthaul core network that serves as the communication infrastructure of
the system.

The rest of the paper is organized as follows: Section 2 presents the
envisioned 5G IoT enabled secure AV-hearing aid framework. Section 3 presents
the proposed real-time lightweight chaotic encryption algorithms including
applied transformations, image, and audio encryption flow diagrams. Section 4
explains the proposed speech enhancement framework including designed enhanced
visually derived Wiener filter (EVWF) and long-short term memory (LSTM) based
lip-reading regression model. In Section 5, the used AV datasets and feature
extraction methodologies are presented. Section 6 presents the performance
evaluation of the proposed AV encryption and speech enhancement algorithms.
Finally, Section 6 concludes this work.

##  2 5G IoT enabled Secure Audio-Visual Hearing-Aid Framework

Modern digital hearing aids are marvels of sophisticated engineering. To hear
modern audio, a low-latency and high datarate wireless solution is needed that
would enable in-ear hearing devices to connect seamlessly [10]. In this
article, a novel 5G IoT enabled secure audiovisual hearing-aid framework is
proposed to acquire desired high quality processed speech in noisy
environments. An example of state-of-the-art 5G IoT architecture is shown in
Figure 1. The IoT enabled devices, supporting a wide variety of applications,
to connect to the Internet, utilizing gateway for connectivity. For gateway
design, different access technologies such as WiFi and 4G LTE could be used.
However, they both are incapable of supporting thousands of connected IoT
devices. WiFi suffers from packet collision and limited quality of service
(QoS), whereas 4G LTE suffers from high delay and high packet loss for large
number of users [11]. In addition, the wireless systems operating in
unlicensed frequency bands require additional network equipment, resulting
extra operation and capital expenditures. The unlicensed solutions are also
prone to congestion with exponential increase in IoT deployment. In contrast,
the next generation 5G wireless networks [12][13] are capable of providing
higher datarates, enhanced mobile coverage, improved user experience at
relatively lower cost, and dense connectivity [14]. Furthermore, to address
aggravating detrimental greenhouse (CO2) gas emissions due to ultra-dense 5G
wireless networks and increased network’s energy consumption, 5G C-RAN is a
widely accepted solution that enables improved environmental sustainability,
OPEX, resource management, and energy efficiency [11].

Figure 1: 5G IoT Architecture Figure 2: Envisioned 5G IoT enabled Secure
Audio-Visual Hearing-Aid Figure 3: Proposed Lip-Reading Driven Deep Learning
Approach for Speech Enhancement

The envisioned 5G IoT enabled secure audiovisual hearing-aid is depicted in
Figure 2. It is to be noted that the computational complexity and real-time
processing issues due to deep learning and big data optimization processes are
addressed by running them in the background on the cloud. The mobile hearing-
aid only transmits the encrypted compressed AV information and receives
encrypted enhanced reconstructed speech in real-time. For end-to-end
communication, an indoor 5G wireless small sized cell and back/fronthaul core
network are proposed as the communication infrastructure of the system [15].
For IoT gateway, we propose the use of an efficient IoT gateway over 5G
wireless system (i.e. an integration of Massive MIMO, mmWave, 5G-Cloud-RAN
(C-RAN), and indoor small sized cells for line of sight (LOS) signalling),
developed and tested in [11]. The developed IoT gateway in [11] promises the
uplink latency of 10ms and 5ms with and without compression respectively. In
addition, it ensures minimum interference between medical IoT devices with low
power consumption. The novelty of these gateways lie in efficient uplink IoT
traffic classification and optimal uplink data (traffic) compression
strategies. This helps in relaxation of uplink traffic burden and results in
efficient utilization of uplink wireless resources. More details on 5G-CRAN
and front/backhaul connectivity are comprehensively presented in [11].

For real-time lightweight audio-visual encryption, PWLCM, Chebyshev map,
secure hash algorithm and a novel S-Box algorithms are utilized. In the
literature, conventional encryption approaches such as advanced encryption
standard (AES) and Rivest–Shamir–Adleman (RSA)/Elliptic Curve (signing) are
suitable for high processing power systems but incompatible with embedded low
power sensor networks. Therefore, lightweight cryptography can potentially
address real-time encryption challenges [16]. In our proposed scheme, the
encrypted audio and video signals are exploited in the cloud by the designed
novel lip-reading driven speech enhancement system, depicted in Figure 3. The
proposed speech enhancement approach leverages the complementary strengths of
both deep learning and analytical acoustic modelling (filtering based
approach) that operates at two levels. In the first level, a novel deep
learning based lip-reading regression model is employed. In the second level,
lip-reading approximated clean-audio features are exploited, using an EVWF,
for estimating the clean audio power spectrum. Finally, the Wiener filter is
applied to the magnitude spectrum of the noisy input audio signal, followed by
the inverse fast Fourier transform (IFFT), overlap, and combining processes to
produce enhanced magnitude spectrum. More details are presented in Section 4.
The proposed AV speech enhancement framework finally transmits the enhanced
encrypted speech to the mobile hearing-aid.

##  3 Proposed Real-Time Lightweight Chaotic Encryption

In the proposed scheme, PWLCM, Chebyshev map, SHA and a novel S-Box algorithms
are effectively used for real-time light-weight encryption. The applied
transformations are briefly explained in the subsequent sections.

###  3.1 Applied Transformations

####  3.1.1 PWLCM

As outlined in shannon novel paper [17], a good encryption scheme is composed
of two stages: (i) Confusion and (ii) Diffusion. In confusion stage, a
correlation between key and ciphertext is made complex. Diffusion means that a
minor change in plaintext should change the corresponding ciphertext
significantly. The proposed algorithm uses PWLCM in confusion process. The
PWLCM can be written as:

$y_{n+1}=f(y_{n},\lambda)=\begin{cases}\frac{y_{n}}{\lambda},&$if
$y_{n}\in[0,\lambda]\\\ \frac{1-y_{n}}{1-\lambda},&$if
$y_{n}\in(\lambda,0.5]\\\ F(1-y_{n}),&$if $y_{n}\in(0.5,1],\\\ \end{cases}$
(1)

where, $y_{n}$ are pseudo-random chaotic values, $y_{n}\in(0,1)$, and
$\lambda$ is the control parameter. Both $\lambda$ and $y_{0}$ serve as an
initial condition and called as key for chaotic pseudo-random number
generation.

####  3.1.2 Chebyshev Map

In [18], Huang et al., proposed a novel key generator method using Chebyshev
map. The Chebyshev map can be defined mathematically as [18, 19]:

$T_{k}(x)=cos(k\times arc\mbox{ }cos(x)),$ (2)

where k = 0, 1, 2, …,N and $x\in[-1,1]$. Huang suggested $k=4$ for less
computation and better use of Chebyshev. In the proposed scheme, k = 4 is used
which defines the Chebyshev function as:

$f(x_{i})=8x^{4}_{i-1}-8x^{2}_{i-1}+1,i=1,2,....N$ (3)

####  3.1.3 Logistic-Sine Map

To overcome the drawbacks of one-dimensional (1D) Logistic map, Zhou et al. in
[20] proposed a novel method of chaotic maps combination. For a larger chaotic
map, the authors combined two exiting 1D Logistic and Sine maps. Logistic-sine
map has many advantages over traditional maps and is mathematically defined as
[20]:

$z_{n+1}=(rz_{n}(1-z_{n})+(4-r)\frac{sin(\pi z_{n})}{4})mod(1),$ (4)

####  3.1.4 Secure Hash Algorithm (SHA)

SHA generates a fix length value known as hash code by applying some function
to the plaintext message. In the literature, SHA has different variants
depending on the size of the output e.g, SHA-1, SHA-256 and SHA-512 for 128,
256, and 512 bits outputs, respectively. In the proposed scheme, we used
SHA-512 such that $H(m)$ = $h$(512 bits). Secret Key in the proposed scheme is
dependent on SHA-512. A minor change in the plaintext generates a completely
different hash and different initial key parameters.

####  3.1.5 Affine Transformation

Affine transformation is a one to one mapping that transforms a unique
plaintext into a unique symbol. The following affine transformation is used in
the proposed scheme:

$AT(w)=\begin{bmatrix}1&1&1&1&1&0&0&0\\\ 0&1&1&1&1&1&0&0\\\ 0&0&1&1&1&1&1&0\\\
0&0&0&1&1&1&1&1\\\ 1&0&0&0&1&1&1&1\\\ 1&1&0&0&0&1&1&1\\\ 1&1&1&0&0&0&1&1\\\
1&1&1&1&0&0&0&1\\\ \end{bmatrix}\times\begin{bmatrix}w_{7}\\\ w_{6}\\\
w_{5}\\\ w_{4}\\\ w_{3}\\\ w_{2}\\\ w_{1}\\\ w_{0}\\\
\end{bmatrix}\oplus\begin{bmatrix}0\\\ 1\\\ 1\\\ 0\\\ 0\\\ 0\\\ 1\\\ 1\\\
\end{bmatrix}$ (5)

where $w_{i}$ are coefficients of $w$ i.e multiplicative inverse modulo
($w^{8}+w^{4}+w^{2}+w^{1}+1)$.

###  3.2 Image Encryption Scheme

The flow chart for the grayscale images using PWLCM, Chebyshev, SHA-512 and
affine transformation is shown in Figure 4. The detailed steps of our proposed
cryptosystem are as follows:

Figure 4: Proposed image encryption algorithm.

  * •

Step 1: Convert colour image $I_{c}$ of size $A\times B$ to gray-scale image
$I_{g}$ and save result in $\psi$.

  * •

Step 2: Apply SHA-512 on gray-scale plaintext image $\psi$ and save
hexadecimal hash value in variable $\theta$.

  * •

Step 3: Select first and last 12 hash values and save in $\kappa_{1}$ and
$\kappa_{2}$.

  * •

Step 4: Convert hexadecimal values saved in $\kappa_{1}$ and $\alpha_{1}$ to
decimal values and store result in $\kappa_{1}$, and $\beta_{2}$,
respectively.

  * •

Step 5: Generate SHA-based initial conditions for PWLCM and Chebyshev using
following equations:

$y_{0}=\dfrac{\kappa_{1}}{2^{48}}$ (6) $x_{0}=\dfrac{\kappa_{2}}{2^{48}}$ (7)

  * •

Step 6: Iterate PWLCM $A$ times and store chaotic values in $\alpha$. Randomly
permute rows of gray-scale image $I_{g}$ using the sequence $\alpha$ and save
values in $I_{rp}$.

  * •

Step 7: Iterate Chebyshev map $B$ times and store chaotic values in $\beta$.
Randomly permute columns of $I_{rp}$ using the sequence $\alpha$ and save
values in $I_{permuted}$

  * •

Step 8: Iterate Logistic-sine map $A\times B$ times and store random values in
$\gamma$.

  * •

Step 9: Apply following operations on $\gamma$:

$R_{1}=\mbox{Mod}(\gamma\times 10^{14},256),$ (8) $R_{2}=floor(R_{1}).$ (9)

  * •

Step 10: Rearrange row-vector $R_{2}$ in matrix form $R$ and Bit-wise XOR
random matrix $R$ with $I_{permuted}$ to get $\phi$.

  * •

Step 11: Apply affine transformation on $\phi$ and store values as a
ciphertext image $C$.

For decryption, encryption steps are followed in the reverse order.

###  3.3 Audio Encryption Scheme

The flow chart for the audio signal using PWLCM, Chebyshev, SHA-512 and affine
transformation is shown in Figure 5. The detailed steps of our proposed
cryptosystem are as follows:

Figure 5: Proposed audio encryption algorithm.

  * •

Step 1: Convert 1D audio signal into 2D of size $A\times B$ and save result in
$\psi$.

  * •

Step 2: Apply SHA-512 on $\psi$ and save hexadecimal hash value in variable
$\theta$.

  * •

Step 3: Select first and last 12 hash values and save in $\kappa_{1}$ and
$\kappa_{2}$.

  * •

Step 4: Convert hexadecimal values saved in $\kappa_{1}$ and $\alpha_{1}$ to
decimal values and store results in $\kappa_{1}$, and $\beta_{2}$,
respectively.

  * •

Step 5: Generate SHA-based initial conditions for PWLCM and Chebyshev using
equations (6-7)

  * •

Step 6: Iterate PWLCM $A$ times and store chaotic values in $\alpha$. Randomly
permute rows of audio vector $A_{g}$ using the sequence $\alpha$ and save
values in $I_{rp}$.

  * •

Step 7: Iterate Chebyshev map $B$ times and store chaotic values in $\beta$.
Randomly permute columns of $A_{rp}$ using the sequence $\alpha$ and save
values in $A_{permuted}$

  * •

Step 8: Iterate Logistic-sine map $A\times B$ times and store random values in
$\gamma$.

  * •

Step 9: Apply operations given in equations (8-9) to $\gamma$.

  * •

Step 10: Rearrange row-vector $R_{2}$ in matrix form $R$ and Bit-wise XOR
random matrix $R$ with $A_{permuted}$ to get $\phi$.

  * •

Step 11: Apply affine transformation on $\phi$ and store values as a
ciphertext audio $C$.

For decryption, encryption steps are followed in the reverse order.

##  4 Speech Enhancement Framework

The state-of-the-art VWF and designed EVWF are depicted in Figure 6 (a) and
(b) respectively. The authors in [8] presented a hidden Markov model-Gaussian
mixture model (HMM/GMM) based two-level state-of-the-art VWF for speech
enhancement. However, the use of HMM/GMM models for the estimation of clean
audio features from visual features and cubic spline interpolation for the
approximation of high dimensional clean audio power spectrum from the
estimated low dimensional audio features are not optimal choices. The HMM/GMM
model suffers from poor generalization and cubic spline interpolation method
fails to estimate the missing power spectral values that leads to a poor audio
power spectrum estimation. In contrast, the designed EVWF addressed the
limitations of state-of-the-art VWF [8] by employing an inverse filter-bank
transformation (i.e. a pseudoinverse of the approximated audio features) for
audio power spectrum estimation as compared to the cubic spline interpolation
method. In addition, the use of LSTM addressed the generalization and accurate
clean speech coefficient estimation issues. The designed EVWF also eliminates
the need for voice activity detection (VAD) and noise estimation. The designed
LSTM based lip-reading model consists of input layer, two LSTM layers, and
output dense layer. In the designed LSTM model, prior visual features were
feeded into the stacked LSTM layers to exploit the existing temporal
correlation. The lower LSTM layer used 250 cells for encoding the input visual
information and passed its hidden state to the second LSTM layer, which has
300 cells. The output of the second LSTM layer was then feeded into the fully
connected (dense) layer which has total 23 neurons with linear activation
function. The designed LSTM model was trained with the objective to minimise
the mean squared error (MSE) between the predicted and the actual audio
features using stochastic gradient decent algorithm and RMSProp optimiser.
More details on EVWF mathematical formulation and LSTM model are
comprehensively presented in our previous work [21].

(a) State-of-the-art visually-derived Wiener filtering [8].

(b) Proposed enhanced visually-derived Wiener filtering

Figure 6: State-of-the-art visually-derived Wiener filtering and proposed
enhanced visually-derived Wiener filtering

##  5 Dataset and Audio-Visual Feature Extraction

For AV encryption and speech enhancement, Grid [22] and ChiME3 [23] corpora
are used. For noisy utterances generation, the clean videos from Grid corpus
were mixed with ChiME3 noises for different SNR levels ranging from -12 to
12dB. For preprocessing, sentence alignment and prior visual frames were used
to stop the model from learning redundant information and improve mapping
between visual and audio features (exploiting their temporal information).

###  5.1 Audio-Visual feature extraction

For audio feature extraction, the input audio signal was sampled at 50kHz and
segmented into N 16ms frames with 800 samples per frame and 62.5% increment
rate. To produce 2048-bin power spectrum, a hamming window and Fourier
transformation was applied, followed by the logarithmic compression to produce
23-D log-FB signal. The visual features were extracted from the Grid Corpus
videos recorded at 25 fps. The video files were processed by extracting a
sequence of individual frames and applying a Viola-Jones lip detector [24] and
object tracker [25]. Furthermore, to ensure appropriate lip tracking,
processed utterances were manually validated. Finally, the 2D-Discrete Cosine
Transformation (2D-DCT) was applied to produce vectors of pixel intensities,
followed by interpolation. More details are presented in our previous work
[21].

##  6 Performance Evaluation

(a) Original Image

(b) Encrypted Image

Figure 7: Encryption Results. It is to be noted that the proposed encryption
scheme completely concealed the plaintext information.

(a) Original Image Histogram.

(b) Encrypted Image Histogram.

Figure 8: Histogram Results. It is to be noted that the obtained histogram of
the encrypted image is flat which is ideally required.

###  6.1 Lip-Images Encryption and Security Analyses

In order to show the effectiveness of the encryption scheme, one test image is
selected with a specific lips position. The proposed lightweight encryption
scheme is applied to the plaintext image and the obtained results are shown in
Figure 7. It can be seen that the proposed scheme completely concealed the
plaintext information. Moreover, the histogram results shown in Figure 8
acquired the flat histogram which is required ideally. However, the histogram
results are not sufficient to prove a complete security of the cryptosystem.
Therefore, a large number of security metrics such as correlation coefficient,
entropy, contrast, energy, NPCR, and UACI defined in our previous work [26,
27, 28] are used. The degree of similarity between adjacent pixels are
generally analysed via correlation coefficient metric. Ideally, correlation in
all directions (horizontal $(H_{CC})$, vertical$(V_{CC})$ and diagonal
$(D_{CC})$) should be close to zero. Entropy is an other important metric
which can evaluate the resistance capability against statistical attacks. For
a good cryptosystem, entropy of a gray-scale encrypted image should be 8 bits
ideally. A contrast of an encrypted image is defined as intensity between a
pixel and its neighbour pixels. In image encryption, higher values of contrast
indicates higher quality of encrypted image. A sum of squared elements (SSE)
in gray level co-occurrence matrix returns energy of an image. The energy
value of a secure encrypted image is desired to be low. In case of complete
constant pixels, energy value is 1. NPCR and UACI show resistant against
differential attacks. Higher values of NPCR and UACI reflects higher
encryption quality. More detail on how these parameters prove the security of
our encryption scheme is defined in detail in our previous work [26, 27, 28].
The correlation plots in vertical directions are shown in Figure 9. From
correlation plots, it is evident that the distribution of adjacent pixels in
vertical direction is uncorrelated as compared to plaintext correlation.
Mathematical values of correlations in vertical, horizontal and diagonal
directions are outlined in Table 1. In the table, lower correlation values
prove the robustness of the proposed scheme. In addition, it can be seen that
the histogram plot for the proposed scheme is uniformly distributed; hence,
assures resistance to statistical attacks as compared to existing algorithms
[29][30]. NPCR greater than 99% also reveals higher security. Consequently,
the security metrics demonstrate the effectiveness and higher security of the
proposed scheme. Lastly, the required encryption/decryption time is less than
25 msec on a 60 GB RAM using MATLAB software. Such a low processing time
confirms that the proposed scheme is lightweight and could be an effective
solution for practical real-time applications. However, further real-time
optimization is an going work.

(a) Original Image correlation coefficient.

(b) Encrypted Image correlation coefficient.

Figure 9: Correlation coefficients in vertical direction. Table 1: Image security assessment. It is to be noted that the correlation values in all directions (horizontal $(H_{CC})$, vertical$(V_{CC})$ and diagonal $(D_{CC})$) are low which proves the robustness of the proposed scheme. In addition, note the resistance capability against statistical attacks and high security, evaluated using Contrast, Energy, NPCR and UACI tests. Security Parameter | Original Frame | Encrypted Frame
---|---|---
$V_{CC}$ | 0.9523 | 0.0044
$H_{CC}$ | 0.9711 | -0.0056
$D_{CC}$ | 0.9677 | 0.0089
$Entropy$ | 7.0025 | 7.9983
$Contrast$ | 0.1049 | 10.4608
$Energy$ | 0.2161 | 0.0156
$NPCR$ | NA | 99.4566
$UACI$ | NA | 33.1561
Table 2: LSTM Training and Testing Accuracy - Comparison For Different Visual Frames. The table presents an overall behaviour of the LSTM model when contextual information (i.e. previous frames) is added. It is to be noted that the LSTM model exploited the temporal correlation effectively but saturated at 18 prior visual frames. | LSTM
---|---
Visual Frames | $MSE_{train}$ | $MSE_{test}$
1 | 0.092 | 0.104
2 | 0.087 | 0.097
4 | 0.073 | 0.085
8 | 0.066 | 0.082
14 | 0.061 | 0.080
18 | 0.058 | 0.078

###  6.2 Speech Enhancement Results

###  6.3 Lip-Reading Results

For lip-reading, multiple prior visual frames (lip images) are used (ranging
from 1 visual frame to 27 prior visual frames). The simulation results are
shown in Table 2. It can be seen that by moving from 1 visual frame to 18
visual frames, a significant performance improvement could be achieved. The
LSTM model with 1 visual frame achieved the MSE of 0.092, whereas with 18
visual frames the model achieved the least MSE of 0.058. LSTM based learning
model exploited the temporal information (i.e. prior visual frames)
efficiently and showed consistent reduction in MSE while going from 1 to 18
visual frames. This is mainly because of its inherent recurrent architectural
property and the ability of retaining state over long time spans by using
gates. More details, critical analysis, and comparisons are comprehensively
presented in our previous work [21].

Table 3: Speech Enhancement Results. It can be seen that at low SNR levels,
EVWF significantly outperformed benchmark SS and LMMSE based speech
enhancement methods.

SNR | PESQ | MOS
---|---|---
| SS | LMMSE | EVWF | SS | LMMSE | EVWF
-12dB | 0.9 | 0.95 | 1.52 | 0.31 | 0.315 | 1.68
-6dB | 1.01 | 1.03 | 1.58 | 0.51 | 0.48 | 1.965
-3dB | 1.17 | 1.18 | 1.60 | 1.17 | 1.165 | 2
0dB | 1.21 | 1.20 | 1.63 | 1.925 | 1.97 | 2.11
3dB | 1.25 | 1.34 | 1.72 | 2.025 | 2.085 | 2.22
6dB | 1.26 | 1.39 | 1.69 | 2.295 | 2.345 | 2.33
12dB | 1.54 | 1.60 | 1.74 | 2.58 | 2.61 | 2.54

####  6.3.1 Objective Test

For objective testing and comparison with the state-of-the-art audio only
speech enhancement methods (spectral subtraction (SS) and Log-Minimum Mean
Square Error (LMMSE)), perceptual evaluation of speech quality (PESQ) is used
to evaluate the quality of restored speech. PESQ is one of the reliable
methods to evaluate speech quality. The PESQ score is computed as a linear
combination of the average disturbance value and the average asymmetrical
disturbance values. The PESQ score ranges from -0.5 to 4.5 corresponding to
low to high speech quality. The PESQ scores for proposed EVWF and state-of-
the-art benchmark audio only speech enchantment approaches are depicted in
Table 3. It can be seen that at low SNR levels, EVWF significantly
outperformed both SS and LMMSE based speech enhancement methods.

####  6.3.2 Subjective Listening Tests

The subjective listening test was conducted in terms of MOS with self-reported
normal-hearing listeners. The listeners were presented with both clean
(target) and enhanced speech, and were asked to rate the re-constructed speech
on a scale of 1 to 5. The five rating choices were: (5) Excellent (when the
listener feels unnoticeable difference compared to the target clean speech)
(4) Good (perceptible but not annoying) (3) Fair (slightly annoying) (2) Poor
(annoying), and (1) Bad (very annoying). The EVWF performance was compared
with two state-of-the-art speech enhancement methods (SS and LMMSE). A total
of 10 listeners took part in the evaluation session. In Table 3, it can be
seen that at low SNRs (-12dB, -6dB, and -3dB), the proposed EVWF outperformed
audio-only speech enhancement methods. On the other hand, for high SNRs, our
AV approach performed comparably to the Audio only approach.

(a) Original Audio

(b) Encrypted Audio

(c) Decrypted Audio

Figure 10: Audio encryption results.

###  6.4 Audio Encryption and Security Analyses

The proposed lightweight chaotic encryption was applied to the enhanced audio
signal. A single bit in the original speech signal was modified and encrypted
via same key [31]. The two ciphered speech signals, i.e., $S_{c}1$ and
$S_{c}2$ were generated and the obtained results are shown in Figure 10. It
can be seen that the proposed scheme completely concealed the audio
information and decrypted speech signal accurately. Furthermore, to check the
robustness and security strength of the proposed audio encryption scheme,
security analysis tests such as number of sample change rate (NSCR), UACI,
correlation coefficient, and key length are conducted, where low correlation
value, 99.95% NSCR, and 33.3% UACI demonstrate the robustness of the proposed
encryption scheme. In addition, the key length of the proposed approach is
much greater than the minimum required length ($2^{100}$) that shows
resistance against brute force attack. Table 4 presents the security analysis
and effectiveness of the proposed scheme.

Table 4: Audio security assessment. Security Parameter | Encrypted Signal
---|---
NSCR | 99.95%
UACI | 33.39%
Correlation coefficient | 0.00022
Key length | $10^{45}$

##  7 Conclusion

The next-generation multimodal (lip-reading driven) hearing-aids stand as a
major enabler for modern digital hearing aids, capable of restoring the
intelligibility and reducing cognitive load in environments where overwhelming
noise is present. However, the real-time implementation of such AV hearing-
aids demand high data rate, low latency, low computational complexity, and
high security. In this paper, we envisioned a 5G IoT enabled secure AV
hearing-aid that leverages the complementary strengths of 5G and IoT
technology to address the aforementioned challenges. As part of the envisioned
technology, two main contributions (AV speech enhancement in the cloud and
lightweight AV encryption) are reported here. The critical analysis in terms
of both speech enhancement and AV encryption demonstrate the potential of the
envisioned technology in acquiring high quality speech reconstruction in
extreme noisy situations and secure mobile AV hearing aid communication.
Specifically, the comparative performance evaluation of the proposed speech
enhancement method under real noisy environments revealed that the proposed
approach significantly outperformed benchmark audio-only approaches at low SNR
with comparable performance at high SNRs. The audio and video encryption
results revealed the effectiveness of the proposed real-time lightweight
encryption scheme in terms of both high security and processing time. The
ongoing and future work includes the software integration of the proposed AV
mobile hearing aid with 5G-CRAN and its hardware prototype implementation for
real-time testing.

## Acknowledgments

This work was supported by the UK Engineering and Physical Sciences Research
Council (EPSRC) Grant No. EP/M026981/1. The authors would like to gratefully
acknowledge Mandar Gogate from the University of Stirling for his contribution
in implementing LSTM driven AV mapping, which was published in our previous
work and cited here for reference.

## References

  * [1] AJ Ruggles and IW Ekoto.  Ignitability and mixing of underexpanded hydrogen jets.  international journal of hydrogen energy, 37(22):17549–17560, 2012\.
  * [2] S Kortlang, S Ewert, H Meister, S Rählmann, J Kießling, et al.  Combination of controlled laboratory tests and structured field trials for a comprehensive evaluation of a model-based hearing aid.  Int J Audiol, 2016.
  * [3] Amir Hussain, Jon Barker, Ricard Marxer, Ahsan Adeel, William Whitmer, Roger Watt, and Peter Derleth.  Towards multi-modal hearing aid design and evaluation in realistic audio-visual settings: Challenges and opportunities.
  * [4] William H Sumby and Irwin Pollack.  Visual contribution to speech intelligibility in noise.  The journal of the acoustical society of america, 26(2):212–215, 1954.
  * [5] Quentin Summerfield.  Use of visual information for phonetic perception.  Phonetica, 36(4-5):314–331, 1979.
  * [6] Harry McGurk and John MacDonald.  Hearing lips and seeing voices.  Nature, 264(5588):746, 1976.
  * [7] Michelle L Patterson and Janet F Werker.  Two-month-old infants match phonetic information in lips and voice.  Developmental Science, 6(2):191–196, 2003.
  * [8] Ben Almajai, Milner.  Visually derived wiener filters for speech enhancement.  IEEE Transactions on Audio, Speech, and Language Processing, 19(6):1642–1651, 2011.
  * [9] Ahsan Adeel, Hadi Larijani, and Ali Ahmadinia.  Random neural network based novel decision making framework for optimized and autonomous power control in lte uplink system.  Physical Communication, 19:106–117, 2016.
  * [10] Richard Einhorn.  Hearing aid technology for the 21st century: A proposal for universal wireless connectivity and improved sound quality.  IEEE pulse, 8(2):25–28, 2017.
  * [11] Navrati Saxena, Abhishek Roy, Bharat JR Sahu, and HanSeok Kim.  Efficient iot gateway over 5g wireless: A new design with prototype and implementation results.  IEEE Communications Magazine, 55(2):97–105, 2017.
  * [12] Mamta Agiwal, Abhishek Roy, and Navrati Saxena.  Next generation 5g wireless networks: A comprehensive survey.  IEEE Communications Surveys & Tutorials, 18(3):1617–1655.
  * [13] Jeffrey G Andrews, Stefano Buzzi, Wan Choi, Stephen V Hanly, Angel Lozano, Anthony CK Soong, and Jianzhong Charlie Zhang.  What will 5g be?  IEEE Journal on selected areas in communications, 32(6):1065–1082, 2014.
  * [14] Naga Bhushan, Junyi Li, Durga Malladi, Rob Gilmore, Dean Brenner, Aleksandar Damnjanovic, Ravi Sukhavasi, Chirag Patel, and Stefan Geirhofer.  Network densification: the dominant theme for wireless evolution into 5g.  IEEE Communications Magazine, 52(2):82–89, 2014.
  * [15] Min Chen, Jun Yang, Yixue Hao, Shiwen Mao, and Kai Hwang.  A 5g cognitive system for healthcare.  Big Data and Cognitive Computing, 1(1):2, 2017.
  * [16] William J Buchanan, Shancang Li, and Rameez Asif.  Lightweight cryptography methods.  Journal of Cyber Security Technology, 1(3-4):187–201, 2017.
  * [17] Claude E Shannon.  Communication theory of secrecy systems.  Bell Labs Technical Journal, 28(4):656–715, 1949.
  * [18] Xiaoling Huang.  Image encryption algorithm using chaotic chebyshev generator.  Nonlinear Dynamics, 67(4):2411–2417, 2012.
  * [19] Xingyuan Wang, Dapeng Luan, and Xuemei Bao.  Cryptanalysis of an image encryption algorithm using chebyshev generator.  Digital Signal Processing, 25:244–247, 2014.
  * [20] Yicong Zhou, Long Bao, and CL Philip Chen.  A new 1d chaotic system for image encryption.  Signal processing, 97:172–182, 2014.
  * [21] Hussain Amir Whitmer William M Adeel Ahsan, Gogate Mandar.  Lip-reading driven deep learning approach for speech enhancement.  arXiv preprint arXiv:1808.00046, 2018.
  * [22] Martin Cooke, Jon Barker, Stuart Cunningham, and Xu Shao.  An audio-visual corpus for speech perception and automatic speech recognition.  The Journal of the Acoustical Society of America, 120(5):2421–2424, 2006.
  * [23] Jon Barker, Ricard Marxer, Emmanuel Vincent, and Shinji Watanabe.  The third ‘chime’speech separation and recognition challenge: Dataset, task and baselines.  In Automatic Speech Recognition and Understanding (ASRU), 2015 IEEE Workshop on, pages 504–511. IEEE, 2015.
  * [24] Paul Viola and Michael Jones.  Rapid object detection using a boosted cascade of simple features.  In Computer Vision and Pattern Recognition, 2001. CVPR 2001. Proceedings of the 2001 IEEE Computer Society Conference on, volume 1, pages I–I. IEEE, 2001.
  * [25] David A Ross, Jongwoo Lim, Ruei-Sung Lin, and Ming-Hsuan Yang.  Incremental learning for robust visual tracking.  International Journal of Computer Vision, 77(1-3):125–141, 2008\.
  * [26] Jawad Ahmad, Muazzam A Khan, Seong Oun Hwang, and Jan Sher Khan.  A compression sensing and noise-tolerant image encryption scheme based on chaotic maps and orthogonal matrices.  Neural Computing and Applications, 28(1):953–967, 2017.
  * [27] Fadia Ali Khan, Jameel Ahmed, Jan Sher Khan, Jawad Ahmad, and Muazzam A Khan.  A novel substitution box for encryption based on lorenz equations.  In Circuits, System and Simulation (ICCSS), 2017 International Conference on, pages 32–36. IEEE, 2017.
  * [28] Jan Sher Khan, Jawad Ahmad, and Muazzam A Khan.  Td-ercs map-based confusion and diffusion of autocorrelated data.  Nonlinear Dynamics, 87(1):93–107, 2017.
  * [29] Jawad Ahmad and Seong Oun Hwang.  Chaos-based diffusion for highly autocorrelated data in encryption algorithms.  Nonlinear Dynamics, 82(4):1839–1850, 2015.
  * [30] Amir Anees, Adil Masood Siddiqui, and Fawad Ahmed.  Chaotic substitution for highly autocorrelated data in encryption algorithm.  Communications in Nonlinear Science and Numerical Simulation, 19(9):3106–3118, 2014.
  * [31] P Sathiyamurthi and S Ramakrishnan.  Speech encryption using chaotic shift keying for secured speech communication.  EURASIP Journal on Audio, Speech, and Music Processing, 2017(1):20, 2017.


# Numerical Approximation of Kramers-Kronig
Relations to Transform Discretized Absorption Data

Patrick D Fitzgerald

(December 20, 2017)

###### Abstract

The Kramers-Kronig relations describe a pair of integral transforms relating
the real and imaginary components of an analytic function in the complex
plane. These relations are particularly useful in extracting refractive index
characteristics of a given physical test species, wherein more direct
measurements are not terribly reliable. In this report, a method of performing
this integral transformation is derived and discussed. Special attention is
given to the precision of these methodologies, and a specific implementation
of this in Matlab is presented.

##  1 Introduction

In the field of spectroscopy, ascertaining material properties such as
absorption coefficient and refractive index is frequently of interest in order
to characterize them, and one way to measure these properties is through
directly passing an electromagnetic pulse through a sample, and observing how
the waveform is distorted and attenuated. Due to the nature of this type of
measurement scheme, called Time Domain Spectroscopy, the resulting absorption
profile is generated with higher certainty than that of the refractive index.
This presents an opportunity to apply some interesting underlying physics and
mathematics to transfer some of that certainty over to the less certain
quantities.

Resulting from the type of measurement, it is most straightforward to
calculate the absorption coefficient directly for a uniform range of
frequencies, and to computationally generate the refractive index from the
absorption data.

###  1.1 Absorption and Refractive Index

Generally speaking, the absorption coefficient and the refractive index are
two electromagnetic properties which help describe how certain materials
behave under the effect of excitation by light. Further implications regarding
the microscopic behaviors and phenomena of the material can be derived from
these quantities. The absorption coefficient, $\alpha$, represents the
proportion of the light passing through a material which it absorbs per unit
length. The refractive index, $n$, determines the amount that light slows down
within the medium. Both of these quantities depend on the frequency of the
light in question, and as such they can be seen as functions of the angular
frequency $\omega$ (rad/s) or (ordinary) frequency $\nu$ (Hz), where
$\omega=2\pi\nu$. Interestingly, these quantities are strongly related, and
have a more obvious representation as the dielectric function

$\varepsilon(\omega)=n(\omega)+i\frac{c}{2\omega}\alpha(\omega),$ (1)

where $c$ is the speed of light in vacuum. Due to physical and temporal
symmetry considerations, this function is assumed to have the property that

$\varepsilon(-\omega)=\varepsilon(\omega)^{*}$ (2)

where the asterisk designates the complex conjugate. Accordingly, the
refractive index and the absorption coefficient become even functions about
the origin. $\varepsilon$ is an analytic function in the upper half of the
complex plane, which means that the results of Cauchy’s complex contour
integral theorems become applicable and exploitable.

###  1.2 The Kramers-Kronig Relations and Derivatives

One result of those contour integral theorems is the following relations,
which may be used to link together the real and imaginary components of an
analytic function $f$:

$\text{Re}(f(z))=\frac{1}{\pi}\mathcal{P}\int_{-\infty}^{\infty}\frac{\text{Im}\big{(}f(z^{\prime})\big{)}}{(z^{\prime}-z)}dz^{\prime}$
(3a) and
$\text{Im}(f(z))=-\frac{1}{\pi}\mathcal{P}\int_{-\infty}^{\infty}\frac{\text{Re}\big{(}f(z^{\prime})\big{)}}{(z^{\prime}-z)}dz^{\prime}$
(3b)

where $z$ is a free parameter (fixed for any given integration), and
$z^{\prime}$ is a dummy integration variable. These relations are called the
Kramers-Kronig relations, or the dispersion relations. It is clear to see that
under most scenarios for $f$, there is a pole present at $z^{\prime}=z$. The
$\mathcal{P}$ represents the Cauchy Principal Value, which indicates that the
integral must be computed carefully around the pole so as to have it converge
on a nontrivial, nonarbitrary value (if it does converge at all). For example,
one way of representing this careful integration for some general function in
the numerator $g$ is

$\mathcal{P}\int_{-\infty}^{\infty}\frac{g(z^{\prime})}{z^{\prime}-z}dz^{\prime}=\lim_{\epsilon\rightarrow
0}\bigg{(}\int_{-\infty}^{z-\epsilon}\frac{g(z^{\prime})}{z^{\prime}-z}dz^{\prime}+\int_{z+\epsilon}^{\infty}\frac{g(z^{\prime})}{z^{\prime}-z}dz^{\prime}\bigg{)},$
(4)

which just means that the integration approaches the pole from either side at
equal rates.111If this choice were not taken, and the pole was integrated
closer to the poles at different rates from the left and the right, then the
whole integral could result in any arbitrary value, so long as each side was
divergent. This choice of equal rates corresponds directly to the area-under-
the-curve perspective of integrals, whereby an equal amount of infinite area
above and below the axis logically allows for cancellation, yielding a finite
answer. This care must be taken for each pole present in the integrand.,222In
addition to specifying equal rates when approaching a singularity, this
prescription also similarly describes the rates at which one must approach
positive and negative infinity. This is not of concern, since the functions
being integrated in this context are assumed to decay very quickly at higher
frequencies, and thus converge to the same value regardless of the rates of
approach.

By applying this to the present context, $f$ may be replaced with
$\varepsilon$. From the resulting relation corresponding to (3a), we can
generate two integral relations, one where the free parameter is relabeled to
$\omega$, the other where the free parameter is relabeled to $\omega_{a}$, and
in both the dummy integration variable is relabeled to $\omega^{\prime}$.
These two relations may be combined by taking their difference and combining
the integrals into one, whose bounds may be halved due to the symmetry present
in $\alpha$, with slight modification to the integrand of each. This
combination and simplification leaves the Singly Subtractive Kramers-Kronig
Relations (SSKKR)[1]:

$n(\omega)=n(\omega_{a})+\frac{c}{\pi}\mathcal{P}\int_{0}^{\infty}\frac{\alpha(\omega^{\prime})(\omega^{2}-\omega_{a}^{2})}{(\omega^{\prime
2}-\omega^{2})(\omega^{\prime 2}-\omega_{a}^{2})}d\omega^{\prime}.$ (5)

This result allows for the determination of the refractive index at an
arbitrary $\omega$ from information derived from the absorption coefficient
over a large region, and yielded relative to the refractive index at the
frequency $\omega_{a}$, appropriately named the anchor frequency. The intent
of this anchor frequency is to select a frequency at which the refractive
index is known with decent certainty, so that the rest of the desired
refractive index values may be defined relative to this reference. Another
benefit of this form is the relative size of the factor which scales $\alpha$
in the integrand: because it decays rapidly away from respective poles, it
reduces the total contribution of the part of the integral which covers the
arbitrarily high $\omega^{\prime}$ values, or equivalently values of
$\omega^{\prime}$ close to zero. This means that truncating the integral’s
bounds produces less total error. This is especially useful in mitigating the
effects of a finite bandwidth of data for $\alpha$.

##  2 Numerically Approximating the SSKKR

Due to the fact that values for the absorption coefficient have only been
measured as a discrete set of data at equispaced frequency values, the
integral in (5) must be approximated by a discrete summation of weighted
absorption values of the form

$n_{j}=n(\omega_{a})+\sum_{i=1}^{N}{}w_{i}^{(j)}\alpha_{i},$ (6)

where $N$ is the total number of data points available for $\alpha$, and the
subscripts on $n_{j}$ and $\alpha_{i}$ are the refractive index and absorption
coefficient corresponding to frequencies $\omega_{j}$ and $\omega_{i}$,
respectively (or equivalently $\nu_{j}$ and $\nu_{i}$, respectively).
$w_{i}^{(j)}$ serves as a weighting coefficient for each of the measured
values of $\alpha$, and will need to be chosen to best approximate the
integral. The subscript on $w$ defines which $\alpha$ value it should scale
(corresponding to $\omega^{\prime}$ in (5)), and the superscript defines where
the value for $n$ is being calculated (corresponding to $\omega$ in (5)).

Since the $\alpha$ data is equispaced, the available methods for standard
numerical quadrature are restricted to schemes which support this structure,
and the fact that no extra intermediate data points are available333 While it
is technically possible to generate more intermediate data points from the
discrete Fourier transform from which all of this data was generated, this
would be an overly complicated way to tackle this problem. Furthermore,
choosing finer meshes is not a proper substitute for an intelligent management
of the precision of the computation. similarly restricts the available
choices. Accordingly, Newton-Cotes (NC) quadrature schemes will primarily be
used herein. Where these standard NC formulas are not used, new formulas will
be derived, but the construction of these will involve the same logic from
which these standard NC formulas were originally derived.

Before delving into the details of the evaluation, let us first make an
alteration of our perspective of this integral: it is more convenient to work
with respect to $\nu$, since the corresponding numerical values for frequency
are easier to interpret as ordinary frequency, instead of angular
frequency444This integral transform will be used primarily on data generated
from terahertz spectroscopy, and the frequencies listed in THz ($10^{12}$Hz)
is often more convenient in this respect. To avoid issues of integrating over
very large domains, the units for $\nu$ are just left in THz within the Matlab
code, and other physical constants like $c$ are converted to this same set of
units.:

$n(\nu)=n(\nu_{a})+\frac{c}{2\pi^{2}}\mathcal{P}\int_{0}^{\infty}\frac{\alpha(\nu^{\prime})(\nu^{2}-\nu_{a}^{2})}{(\nu^{\prime
2}-\nu^{2})(\nu^{\prime 2}-\nu_{a}^{2})}d\nu^{\prime}.$ (7)

No factor of $2\pi$ is generated from this change of variables inside $n$ or
$\alpha$, since as far as these functions are concerned, their values are just
labeled by the corresponding $\omega$ and $\nu$ values, but are not explicitly
dependent on them, in a numerical sense.

Figure 1: A sample of $S$ functions, where the $\nu$ and $\nu_{a}$ are placed
at a variety of locations within the domain. The dashed vertical lines
correspond to $\nu^{\prime}=\nu$, and the dotted vertical lines correspond to
$\nu^{\prime}=\nu_{a}$.

Figure 2: The same as in Figure 2, except that $\nu_{a}$ is fixed at
$\nu^{\prime}=1.3$, while $\nu$ alone is varied. This fixed nature of
$\nu_{a}$ is properly representative of the different scenarios which are seen
in this discussion.

###  2.1 Consideration of Different Cases

In attempting to approximate anything, it is often a good first step to
appreciate what details are most important, so that they may be the primary
focus of the approximation. After those are understood and accounted for, it
is also important to understand any sort of cases which require specific
behavior to maintain a high quality of approximation.

Regarding the former, in order to gain an appreciation of the important
details of this integrand, we may simply focus on just the function which
scales $\alpha$, which we may name $S(\nu^{\prime})$, where the value for
$\nu$ and $\nu_{a}$ are simply interpreted as parameters:

$S(\nu^{\prime})=\frac{(\nu^{2}-\nu_{a}^{2})}{(\nu^{\prime
2}-\nu^{2})(\nu^{\prime 2}-\nu_{a}^{2})}$ (8)

A variety of scenarios for $S(\nu^{\prime})$ are presented in Figures 2 and 2,
which each have different values for $\nu$ and $\nu_{a}$. From this graphical
representation, it is clear to see that the region just next to the poles will
contribute strongly to the integral and therefore also to the weights
$w_{i}^{(j)}$. Outside these neighborhoods, the functions are relatively slow
to change, and in many cases, rather small in magnitude. From these
observations, it is reasonable to assume that more care needs to be given to
the regions around the poles555 In addition to needing care due to the Cauchy
Principal Value, this region needs extra care due to its asymptotic behavior.
Since all of the NC schemes assume some sort of terminating polynomial form in
the domain of integration, they are especially poorly equipped to handle
functions which are not well approximated by their lower degree terms, as is
characteristic of asymptotic forms. , and a less nuanced attack will suffice
elsewhere.

Since for a given set of computations for determining $n(\nu)$, the value for
$\nu_{a}$ will remain fixed and the value for $\nu$ will iterate over a range
of values, there are necessarily more cases than just ’near pole’ and ’not
near pole’. The ’near pole’ case must be broken into three categories: only
one pole nearby, two distinct poles nearby, and one second-degree pole nearby.
In order to properly discuss these cases, let us label them as shown in Table
1:

Table 1: Node proximity-to-poles labels Case | Description
---|---
I | Not near any poles
II | Near one pole
III | Near two distinct poles
IV | Near one second degree pole

###  2.2 Case I: Regions Away from Poles

Outside any asymptotic regions, we will simply apply standard composite NC
schemes. Since this case represents all the locations where nothing especially
drastic is happening, the regions to which this composite scheme will need to
be applied might not always be of convenient or consistent size666This region
essentially acts as a filler, in order to gather up the contribution of the
$\alpha$ values not described elsewhere by more nuanced means.. Accordingly, a
plan needs to be in place on how to handle nonstandard size. If a minimum NC
order has $N_{I,min}+1$ points, then order could be increased to represent
more points ($N_{I}+1$) without an increase in truncation error. This would
eventually create significant round-off error, so trying to increase the order
of many segments by a small amount is better than increasing one drastically.
If, however, this Case I region is too small due to being pinched between a
region of Case II and that of another Case II or the edge of the available
data, then the Case II region should be extended to fully represent that
otherwise-Case-I region. A depiction of this Case II extension can be most
easily seen at the top and bottom corners of Figure 3.

In these segments of Case I points, the weight values will contain the
inherent NC weights, and the contribution of the integrand weighting function,
$S(\nu^{\prime})$, will manifest simply as scaling the NC weights by the $S$
value at the respective $\nu^{\prime}$ values. Since these segments are far
from poles, $S$ will be well behaved.

###  2.3 Cases II-IV: Regions Next to Poles

Case II regions will require a method of approximating the respective
integrals with a scheme with the same underlying structure as the NC scheme.
If inside this region of $N_{II,min}+1$ points, and the corresponding integral
over them, it becomes beneficial to find an interpolating function for
$\alpha$ in this range. In order to do this, let us employ a Lagrange
interpolating polynomial, in the same fashion as the NC scheme does. Thus,

$\alpha(\nu^{\prime})\approx\sum_{i=0}^{N_{II}}\bigg{(}\alpha_{i}\prod_{\begin{subarray}{c}k=0\\\
k\neq{}i\end{subarray}}^{N_{II}}\frac{\nu^{\prime}-\nu_{k}}{\nu_{i}-\nu_{k}}\bigg{)}$
(9)

in the region from $\nu_{0}$ to $\nu_{N_{II}}$ (which may easily be relabeled
to represent any such integration without loss of generality). Since there is
only one pole to be careful of, it would be most convenient to evaluate the
rest in the previous manner. To this end, we may expand the integrand
weighting function $S$ to be

$S(\nu^{\prime})=\frac{-1}{2\nu}\bigg{(}\frac{1}{\nu^{\prime}+\nu}\bigg{)}+\frac{1}{2\nu}\bigg{(}\frac{1}{\nu^{\prime}-\nu}\bigg{)}+\frac{-1}{2\nu_{a}}\bigg{(}\frac{1}{\nu^{\prime}-\nu_{a}}\bigg{)}+\frac{1}{2\nu_{a}}\bigg{(}\frac{1}{\nu^{\prime}+\nu_{a}}\bigg{)},$
(10)

by partial fraction decomposition. In this form, it becomes very easy to
extract one of the pole’s respective components and leave the rest to be
evaluated as in the Case I scenario. For the remaining pole, which would be at
either $\nu^{\prime}=\nu$ or $\nu^{\prime}=\nu_{a}$, labeled here as $\nu_{p}$
for generality, the remaining integrand becomes

$(\pm)_{p}\frac{1}{2\nu_{p}}\frac{1}{(\nu^{\prime}-\nu_{p})}\sum_{i=0}^{N_{II}}\bigg{(}\alpha_{i}\prod_{\begin{subarray}{c}k=0\\\
k\neq{}i\end{subarray}}^{N_{II}}\frac{\nu^{\prime}-\nu_{k}}{\nu_{i}-\nu_{k}}\bigg{)},$
(11)

where $(\pm)_{p}$ just defines the sign of the term depending on whether the
pole is at $\nu$ or $\nu_{a}$. Upon commuting sums and factors around, and
comparing the integral this represents to the desired weighted sum, we see
that the part of the weights $w_{i}$ created from the pole at $\nu_{p}$ is

$(w_{i})_{p}=(\pm)_{p}\frac{1}{2\nu_{p}}\bigg{(}\prod_{\begin{subarray}{c}k=0\\\
k\neq{}i\end{subarray}}^{N_{II}}\frac{1}{\nu_{i}-\nu_{k}}\bigg{)}\int_{\nu_{0}}^{\nu_{N_{II}}}\bigg{(}\prod_{\begin{subarray}{c}k=0\\\
k\neq{}i\end{subarray}}^{N_{II}}(\nu^{\prime}-\nu_{k})\bigg{)}\frac{1}{\nu^{\prime}-\nu_{p}}d\nu^{\prime}.$
(12)

Since the current integrand is factored, we can guarantee that the pole will
cancel out creating only a removable discontinuity, so long as $i\neq{}p$. In
this majority of cases, the integral can be expanded and evaluated directly,
since it is just a polynomial. However, if $i=p$ then the polynomial will need
to be long divided to separate out the quotient of the division as a
polynomial (again, easily integrated), and the portion which did not fully
divide: the remainder. The remainder’s contribution can be solved
analytically, since it will take the form of a constant divided by the pole,
where applying the Cauchy Principal Value, we find that

$\int_{\nu_{0}}^{\nu_{N_{II}}}\frac{A}{\nu^{\prime}-\nu_{p}}d\nu^{\prime}=\lim_{b\rightarrow\nu_{p}}\log|\nu^{\prime}-\nu_{p}|\Big{|}_{\nu_{0}}^{b}+\lim_{b\rightarrow\nu_{p}}\log|\nu^{\prime}-\nu_{p}|\Big{|}_{b}^{\nu_{N_{II}}}=\log\bigg{|}\frac{\nu_{N_{II}}-\nu_{p}}{\nu_{0}-\nu_{p}}\bigg{|}.$
(13)

So long as the pole does not occur at either of the ends of the region, these
two scenarios will always work. This will hold true regardless of stretching,
which will further remove the points from the ends of the integration domain.
To ensure this is the case, no fewer than three points may be present within
any Case II region.

Figure 3: A graphical representation of the weights matrix $w$, where the
value for $w_{i}^{(j)}$ is the (j,i)th element in the matrix. The different
cases are labeled below the matrix. The black boxes designate the locations of
the poles for a given iteration (row of matrix). The vertical line of poles is
the fixed pole at $\nu^{\prime}=\nu_{a}$, and the diagonal line of poles is
the pole which scans across the data at $\nu^{\prime}=\nu$. Where the two
lines get close exemplifies when the more nuanced integration methods are
needed. Regions which are still classified by Case I but are too close to a
region of Case II or an end of the available data, the neighboring Case II
domains are simply increased so as to avoid needing to use a lower order NC
approximation due to low number of data points.

A very similar approach may be taken for the Case III regions. Since two
distinct poles are present in the region, if $\alpha$ is again represented
with a Lagrange interpolating polynomial and $S$ is expanded again, there will
only be one pole in each integral and the results reduce down to those of Case
II. The two poles at $\nu^{\prime}=-\nu$ or $\nu^{\prime}=-\nu_{a}$ can again
be treated in the Case I way. The only reason for separating this out into a
new case is mostly to account for the more dynamic number of nodes included,
and so it is more straightforward to place neighboring regions consistently.
This variability in included nodes can be seen in the blue section of Figure
3, where each of these regions contains anywhere from $N_{II,min}+2$ to
$2N_{II,min}+1$ points.

Finally, Case IV is even simpler. If the two poles are at the same location,
that means that $\nu=\nu_{a}$. Looking at the original definition for the
integral transform, this causes $S$ to identically go to zero, making the
entire integral go to zero similarly. Accordingly, there is no need to worry
about different regions and arranging them for this row of the $w_{i}^{(j)}$
matrix. This can be seen in Figure 3 as the single red row.

##  3 Error Analysis

As described by the derivation shown by Isaacson and Keller[3], the truncation
error, $R_{M}$ created by replacing $\alpha$ with an $M^{th}$ degree Lagrange
interpolating polynomial (with $M+1$ points) is [3]

$R_{M}(\nu^{\prime})=\bigg{(}\prod_{j=0}^{j=M}(\nu^{\prime}-\nu_{j})\bigg{)}\alpha[\nu_{0},\nu_{1},\nu_{2},...,\nu_{M},\nu^{\prime}],$
(14)

where the latter term is the divided difference of $\alpha$. Since $\alpha$ is
a well behaved function, we may assume the continuity of this divided
difference regardless of $M$. Correspondingly, the error $E_{M}$ created by
this in the standard Newton Cotes scheme is the integral of this truncation
error,

$E_{M}=\int_{\nu_{0}}^{\nu_{M}}\bigg{(}\prod_{j=0}^{j=M}(\nu^{\prime}-\nu_{j})\bigg{)}\alpha[\nu_{0},\nu_{1},\nu_{2},...,\nu_{M},\nu^{\prime}]d\nu^{\prime}.$
(15)

Upon making a substitution into the integral of $\nu^{\prime}=\nu_{0}+ht$,
where $h$ is the spacing between the nodes, Isaacson and Keller work out a
clean closed-form solution for the error for the standard NC scheme as

$E_{M}=\frac{K_{M}}{(M+2)!}h^{M+3}\alpha^{(M+2)}(\xi)=\mathcal{O}(h^{M+3})$
(16)

for even $M$, and

$E_{M}=\frac{K^{\prime}_{M}}{(M+1)!}h^{M+2}\alpha^{(M+1)}(\xi)=\mathcal{O}(h^{M+2})$
(17)

for odd $M$. In both of these, $\xi$ is some value between $\nu_{0}$ and
$\nu_{M}$. This statement is so concise due to the mean value theorem. The
fact that the even and odd $M$ values have different relative accuracy, with
respect to order of $h$, is from the convenient structure available inside the
integrand, which those authors were able to exploit.

This covers Case I, but we now need to characterize the error associated with
the pole integrals investigated in this report. For Case II, there is only one
pole inside the domain, so the truncation error of this pole’s integral will
take the form

$R_{M}(\nu^{\prime})=\bigg{(}\prod_{\begin{subarray}{c}j=0\\\
j\neq{}p\end{subarray}}^{j=M}(\nu^{\prime}-\nu_{j})\bigg{)}\alpha[\nu_{0},\nu_{1},\nu_{2},...,\nu_{M},\nu^{\prime}],$
(18)

and therefore the error associated with this

$E_{M}=\int_{\nu_{0}}^{\nu_{M}}\bigg{(}\prod_{\begin{subarray}{c}j=0\\\
j\neq{}p\end{subarray}}^{j=M}(\nu^{\prime}-\nu_{j})\bigg{)}\alpha[\nu_{0},\nu_{1},\nu_{2},...,\nu_{M},\nu^{\prime}]d\nu^{\prime},$
(19)

and upon using that same substitution again, $\nu^{\prime}=\nu_{0}+ht$, we
have

$E_{M}=h^{M+1}\int_{0}^{M}\bigg{(}\prod_{\begin{subarray}{c}j=0\\\
j\neq{}p\end{subarray}}^{j=M}(t-j)\bigg{)}\alpha[...,\nu_{0}+ht]dt=\mathcal{O}(h^{M+1}).$
(20)

Doing the same for the two distinct roots of Case III through the same
procedures, we have

$E_{M}=h^{M}\int_{0}^{M}\bigg{(}\prod_{\begin{subarray}{c}j=0\\\
j\neq{}p_{1}\\\
j\neq{}p_{2}\end{subarray}}^{j=M}(t-j)\bigg{)}\alpha[...,\nu_{0}+ht]dt=\mathcal{O}(h^{M}).$
(21)

The convenient pattern present in the default NC scheme allowed for the
improvement of the odd $M$ case, but for the rest of these, that pattern is
lost due to the division by the respective poles, and therefore no
simplification would be straightforward or easy to find. For this same reason,
it would also be difficult to apply the mean value theorem to get a concise
bound for the error in terms of $\xi$. However, this is not an issue, since
the scale factors (as functions of $M$) present in these error terms are
necessarily finite, since they are integration of polynomials of finite degree
scaled by a well behaved function ($\alpha$), evaluated over finite bounds.
Accordingly, these error bounds are indeed the listed orders in $h$.

If we were to consider Case IV in this manner, we would arrive at a potential
issue where the integrand is not finite everywhere, since only degree of the
double pole would cancel. Without further thought, this would present a
serious issue with respect to bounding the error with the mean value theorem.
However, this approach is not needed, since this case identically goes to
zero, so the error is necessarily zero.

Since the minimum number of points present in each case has been described as
having a minimum, we can list out a worst case bound for the truncation error
with respect to $h$ for each of the regions, as shown in Table 2. The error
mentioned for each case only really describes the part of the integral which
was not evaluated with a simpler case’s integration scheme.

Table 2: The orders of the errors inherent in the different types of integration used. The worst case error is described for the minimum size of each region, so that when a region is stretched, the truncation error only shrinks, and therefore is still within the listed order. Case | $M_{min}$ | $E_{worst}$
---|---|---
I |  $N_{I}$ even | $\mathcal{O}(h^{N_{I}+3})$
I |  $N_{I}$ odd | $\mathcal{O}(h^{N_{I}+2})$
II | $N_{II}$ | $\mathcal{O}(h^{N_{II}+1})$
III | $N_{II}+1$ | $\mathcal{O}(h^{N_{II}+1})$
IV | N/A | 0

##  4 Implementation

These procedures were implemented inside of Matlab, using $N_{I}+1=3$ and
$N_{II}+1=5$, which ensures that the maximum error for this implementation is
$\mathcal{O}(h^{5})$. Anything much higher would begin to be unreasonable, and
anything too much lower would start having issues. An outline for the
different functions made for this, and how they interact and depend on one
another is depicted in a flow chart in Figure 4.

Figure 4: An outline of the Matlab code implementing the procedures discussed.
The arrows shown indicate their dependencies on the other functions created
for this purpose.

The source code for all of these files may be found in Appendices A-I, with an
additional Appendix J which acts as a script/wrapper to call the functions
properly, assuming the right data is available to process. Table 3 details
which files are in which appendix.

Table 3: My caption Appendix | Function/Script
---|---
A | GenerateFullWeights()
B | OrganizeRegions()
C | OrganizeRegions_I()
D | GenerateNewtonCotesWeights()
E | GenerateSinglePoleWeights()
F | PolynomialIntegration()
G | PolynomialDivision()
H | ReExpandCoefficient()
I | GetExpandedCoefficients()
J | Test_Script.m

The description of what each file does can be found in Table 4.

Table 4: My caption Function/Script | Purpose
---|---
GenerateFullWeights() |  Generates the weight matrix to multiply $\alpha$
OrganizeRegions() |  Handles breaking the elements of $\alpha$ into regions and cases
OrganizeRegions_I() | Handles the placement and stretching of Case I regions
GenerateNewtonCotesWeights() | Generates the Case I NC weights for any number of points
GenerateSinglePoleWeights() | Generates the Case II weights for any number of points
PolynomialIntegration() | Performs integration on arbitrary degree polynomials
PolynomialDivision() | Long divides a polynomial by a specified root
ReExpandCoefficient() |  Converts polynomials of $x^{n}$ into a different basis $(x-a)^{n}$
GetExpandedCoefficients() | Generates a vector representation of a polynomial
Test_Script.m | A script to demonstrate how the functions should be used

###  4.1 Notes on Implementation

In addition to the lower level approach needed to actually implement the
theoretical formulas into code, there were a couple minor details which might
not be self-explanatory.

Firstly, there is the major use of a certain scheme of polynomial
manipulation, including expanding, long dividing, and integrating. These are
simply writing a given polynomial in a certain basis, often $(x)^{n}$ where
$n$ may take any integer value from zero upwards. This vector representation
made for a relatively fast means for computing things like large product of
first degree terms. This methodology is extended straightforwardly to
integrating, since the place in the vector simply indicates what the factor
becomes, and where it gets shifted to. The general scheme for writing a
polynomial is

$a_{0}+a_{1}x+a_{2}x^{2}+a_{3}x^{3}+...=[a_{0},a_{1},a_{2},a_{3},...].$ (22)

A basis centered at $x=a$ would take the form

$b_{0}+b_{1}(x-a)+b_{2}(x-a)^{2}+b_{3}(x-a)^{3}+...=[b_{0},b_{1},b_{2},b_{3},...].$
(23)

From this notation it follows that multiplication by some power of $(x-a)$ in
a basis centered at $a$ acts as essentially a bit shift to the right of
corresponding amount. Therefore, it follows that multiplication of one vector
by another could be rewritten as a sum of simple vector-times-power-of-x
terms, which is straightforward from that perspective. The rest of the
functionality should be similarly straightforward.

One more interesting feature is the methods by which the long division is
done. A given polynomial is rewritten in the same basis as the divisor term,
so that the division becomes trivial and is just a bit shift to the left,
keeping the constant term off as the remainder. After the shift, the remaining
polynomial can be returned to its original basis and yielded as the quotient
of the division.

On a different note, the two functions which are frequently called upon to
generate the Case I and Case II weights are built with persistent variables.
In Matlab this allows the data stored in a variable to persist to future calls
of that function. This functionality is used to not regenerate data if it is
already stored, which just requires that the same spatial parameters are held
constant between a given call. This drastically reduces the run-time of the
algorithm, because it doesn’t need to recalculate these smaller weights as
frequently.

One last detail regards how the Case II and Case III regions were assumed to
be oriented within their minimum-sized sizes. Although this doesn’t matter
given the minimum sizes used here are both odd, if they were changed to
different values then it is useful to note that a left-sided centering is
used. This means that if a minimum-sized region spans an even number of
points, then it will be placed on the left of the two spots in the center,
instead of the right one. Accordingly, preference is given to stretch to the
left before stretching to the right in order to compensate slightly. This
could have been done more generally, by choosing which side offered the fewer
number of stretched regions, but that would have required a much uglier case
analysis within the code.

##  5 Results and Discussion

Two examples of the weight matrix generated by the appended Matlab code are
shown in Figures 5 and 6. The strongest features are, as expected, highly
concentrated on the locations of the poles.

Using some actual data ($\alpha$ in Figure 7) with the weight matrix shown in
Figure 6, we use the script file (Appendix J) to calculate and display the
corresponding refractive index $n$, which is displayed in Figure 8.777It is
not recommended to use this data as a reference for the absorption and
refractive index of water. While this data is accurate to itself, it us not
necessarily absolutely physically accurate.

Figure 5: A visualization of the weight matrix for a small set of alpha
values. This smaller set was chosen to make the patterns more visible. Figure
6: A visualization of the weight matrix for a larger set of alpha values. This
perspective allows for easier appreciation of how localized the integration
really is, since the poles are really the only visible contributors. Figure 7:
A sample of data for the absorption coefficient $\alpha$. Figure 8: The
refractive index corresponding to the $\alpha$, compared to the original
refractive index and an externally corrected refractive index. Note that
before any correction, anchor frequency would lock the resultant refractive
index to the original refractive index at that point–however, here the
corrected value was used as the reference instead for more direct appreciation
of how the integration did. Recall that the reference only contributes an
overall vertical shift, so no significant change in the form of $n$ was made
under this modified display.

From these figures, it is apparent that the integral transform was successful
in pulling out the refractive index, but there are significant features around
the anchor frequency which are not satisfactory. It would appear the points
which would be heavily influenced by Case III scenarios (when the two poles
are very close) have a great deal of noise in that neighborhood. This could be
a simple issue of a typo in the code, or it could be more significant,
pointing to a systematic problem with the approach used in evaluating those
regions.

Upon further investigation, it becomes clear that the single pole integral
results in quantities which do not depend on $h$. This makes sense from the
formulas derived, but as far as an integration technique, seems quite bizarre.
The next steps for this code need to be to re-examine the implementation for
mistakes, and upon not finding any, re-examine how the result can be good when
the points contributing the most to the integral are not affected by their
spacing (and therefore their quantity).

##  6 Conclusions

In this report, the singly-subtractive Kramers-Kronig relations were prepared
for general use in extracting information from the absorption measurements to
ascertain a better understanding of the refractive index, which otherwise has
a significant degree of indeterminacy. A theoretical groundwork was laid for
the different scenarios that would be encountered when approximating the
absorption coefficient with Lagrange interpolating polynomials, and schemes
were developed for evaluating the integrals in these scenarios with heightened
care taken around the iteratively moving poles of the integral transform. An
implementation was discussed briefly, with the source code accompanying the
report, and qualitative results confirming its functionality were presented.
Features were discovered which seem counter intuitive for results of
integration, but they should really be a lot worse if they are actually a
problem. Further investigation will follow, digging into how significant of a
problem this observed noise could be.

## References

  * [1] Bernier, Maxime, et al. “Accurate Characterization of Resonant Samples in the Terahertz Regime Through a Technique Combining Time-Domain Spectroscopy and Kramers-Kronig Analysis.” IEEE Transactions on Terahertz Science and Technology, vol. 6, no. 3, 2016, pp. 442-450., doi:10.1109/tthz.2016.2535244.
  * [2] Burden, Richard L., et al. Numerical analysis. 10th ed., Cengage Learning, 2016.
  * [3] Isaacson, Eugene, and Herbert Bishop Keller. Analysis of Numerical Methods. Wiley, 1966. Chapter 7.

See pages - of appendices.pdf