# Activity Monitoring of Islamic Prayer (Salat) Postures using Deep Learning

Anis Koubaa1,5,3, Adel Ammar1, Bilel Benjdira1, 6, Abdullatif Al-Hadid1, Belal
Kawaf1,
Saleh Ali Al-Yahri1, Abdelrahman Babiker, Koutaiba Assaf, Mohannad Ba Ras

Email: (akoubaa, aammar<EMAIL_ADDRESS>1Robotics and Internet-of-
Things Lab (RIOTU), Prince Sultan University, Riyadh, Saudi Arabia.
2 CISTER Research Centre, ISEP, Polytechnic Institute of Porto, Porto,
Portugal
5Gaitech Robotics, China.  6Research Laboratory SEICT, LR18ES44. National
Engineering School of Carthage, Tunisia.

###### Abstract

In the Muslim community, the prayer (i.e. Salat) is the second pillar of
Islam, and it is the most essential and fundamental worshiping activity that
believers have to perform five times a day. From a gestures’ perspective,
there are predefined human postures that must be performed in a precise
manner. However, for several people, these postures are not correctly
performed, due to being new to Salat or even having learned prayers in an
incorrect manner. Furthermore, the time spent in each posture has to be
balanced. To address these issues, we propose to develop an artificial
intelligence assistive framework that guides worshippers to evaluate the
correctness of the postures of their prayers. This paper represents the first
step to achieve this objective and addresses the problem of the recognition of
the basic gestures of Islamic prayer using Convolutional Neural Networks
(CNN). The contribution of this paper lies in building a dataset for the basic
Salat positions, and train a YOLOv3 neural network for the recognition of the
gestures. Experimental results demonstrate that the mean average precision
attains 85% for a training dataset of 764 images of the different postures. To
the best of our knowledge, this is the first work that addresses human
activity recognition of Salat using deep learning.

##  I INTRODUCTION

Human Activity Recognition has been extensively investigated using different
techniques including sensing technologies [1, 2, 3], computer vision [4], and
more recently using deep learning [5, 6]. It consists of the classification of
the activity of a person from the data collected from sensors (e.g.
accelerometer, camera, and laser scanner ). The advances in human activity
recognition have enabled different applications in different areas such as
healthcare, sports activities, violence detection, older people monitoring,
postures’ recognition, to name a few.

In this paper, we consider a particular human activity application with
special interest to the Muslim community around the world, namely the
recognition of postures of the Islamic prayer, also known as Salat. Salat is
the second pillar in Islam and is the most important worshipping activity that
is repeated five times a day by any Muslim. Besides its spiritual value, it
consists of a series of postures that must be executed in a predefined
sequence as instructed by the Prophet Mohamed, peace be upon him. The
motivation of this work is that several people (e.g., kids, beginners) are
likely not to perform the postures correctly. There are different reasons for
this to happen like, for example, the lack of proper learning of prayer
movements, or being careless. Besides, another requirement in Salat is the
tranquility; that is, every posture should be executed in a proper and
sufficient amount of time while making invocations and reading the Quran. Only
a few research works in the literature have attempted to address the activity
recognition of Salat. In [1, 2], the authors used smartphone technology to
recognize Salat activities. In [3, 7], the authors assessed the activities of
Salat using Electromyographic (EMG) signals. On the other hand, several
research works used deep learning approaches to address the problem of human
activity recognition in different applications such as sports, healthcare,
well-being using wearable sensors [5, 6, 8]. To the best of our knowledge,
there is no previous work that addressed the problem of activity monitoring of
Islamic prayer using deep learning approaches. In this paper, we propose a
novel solution to the posture recognition of Salat using convolutional neural
networks (CNN), which aims at identifying the four basic postures of Salat
namely, Standing (Qiam), Bowing (Ruku), Prostration (Sujud), Sitting (Julus),
using state-of-the-art CNN algorithms. CNN has been used in a wide variety of
applications such as vehicle detection [9, 10], semantic segmentation of urban
environments [11], and self-driving vehicles [12]. The contribution presented
in this paper represents the first step towards the main objective of this
project that consists in developing an AI-based tool for the assessment of the
Islamic prayer and an assistive system that helps beginners and kids to
correct the postures during Salat. It has to be noted that we only consider
the recognition of correct postures in our neural network model, and we ignore
wrong positions, which is left as a future extension of this work in which we
will address the anomalies during prayers. The principal added value of this
paper is the development of a specific dataset for the human activities in the
context of Islamic prayer that provides four classes related to each of the
basic postures of Salat. Furthermore, we trained the state-of-the-art YOLOv3
neural network on the constructed dataset to accurately detect and recognize
the four postures of Salat. We have tested the resulting trained network on
videos showing people praying, and all the postures were most of the time
correctly recognized. The remainder of the paper is organized as follows.
Section II discusses the literature on human activity recognition and also
related works on Islamic prayer monitoring using sensing activities. Section
III gives a brief overview of the YOLOv3 algorithm. Section IV presents the
different postures of Salat and the corresponding dataset. Section V discusses
the experimental results. Finally, Section VI concludes the paper and presents
future extensions.

##  II RELATED WORKS

The earliest work that targets the task of automatic recognition of prayer
movements is the paper introduced in 2009 by El-Hoseiny et al.[13]. They used
a camera to capture the side view of the prayer. Then, they used morphological
operations to extract the polygon corresponding to the contour of the prayer
body. The polygon information is used to calculate the backbone axis angle and
to extract four key-points of the human body. The backbone axis angle is the
angle between the y-axis and the back-axis of the detected body. The four key
points determined from the polygon are the center point of the polygon, the
ankle point, the head point, and the back point. Based on the coordinates of
the key-points and the value of the backbone axis angle, the prayer postures
and movements can be determined using a set of given inequations. The method
is totally based on hand-crafted features without any use of machine learning
classifiers. This work is the only work that used ordinary camera sensor for
this task. The other works leveraged the use of accelerometer and Kinect
sensors. Al Ghannam et al.[2] used a set of classifiers for prayer activity
monitoring and recognition based on accelerometer data collected from a mobile
phone. They used three machine learning classifiers: J48 Decision Trees, IB1
(Instance-Based Learning ) Algorithm, and Naive Bayes. The accuracy of all
these algorithms exceeded 90%. Eskaf et al.[14] presented a framework for
daily life activities (sitting, standing, walking, …) from accelerometer data.
Then, they aggregated these simple activities for recognition of prayer using
supervised machine learning classifiers. Ali et al.[15] presented similarly a
system for automatic recognition and monitoring of prayer postures based on
triaxial accelerometer data collected from a smartphone. They added an
analysis of group prayer activities using dynamic time warping.

In [16], the authors proposed a human activity recognition technique based on
a deep learning model designed for low power devices. The recognition gives
significant contextual data for well-being. Alobaid et al.[1] also studied the
use of mobile accelerometer data for prayer activity recognition. They made a
performance comparison of three feature extraction approaches and eight
machine learning classifiers. They concluded that Random Forest is the most
appropriate for this task, with an accuracy of 90%. They proposed a two-level
classifier to solve the confusion between two similar prayer stages and to
improve the accuracy to 93%. They added a brief study on the effect of
personal characteristics like height and age on the performance of the
classifier. Jaafer et al.[17] used another sensor, the Kinect RGB-Depth
camera. Two Kinect sensors are put in a fixed position on the body, and the
skeleton information is extracted using Kinect Software Development Kit. They
used the Hidden Markov Model as the machine learning classifier for learning
the prayer movements from the skeleton information.

We can see from the above works that only one work had used the ordinary
camera sensor, although its practicability and easiness of use compared to
other sensors. Moreover, this work used only geometric properties without
profiting from the potential of machine learning algorithms. The significant
advance noted in machine learning and especially the emergence of deep
learning algorithms will surely make a significant improvement in the
efficiency of prayer postures recognition from a video camera. This is why we
targeted this limitation in the current state of the art approaches. We
selected one of the most performing one-stage deep learning algorithms in
object detection task (YOLOv3[18]), and we applied it in the context of prayer
posture recognition. We considered in our study that the camera could be
placed in different views from the prayer and not only in a fixed position
from him.

##  III ALGORITHMS BACKGROUND

Among the deep learning algorithms used in computer vision, YOLOv3[18] is the
most attractive for the task of prayer posture recognition. Two facts are
defending this choice. Firstly, YOLOv3[18] had proven its efficiency compared
to other object detection algorithms[9]. Secondly, it has an efficient
inference time (up to 45 frames per second). This allows a real-time
recognition of the prayer postures. In the following subsections, we present
the architecture of the first version of YOLO (YOLOv1[19]) and the different
improvements implemented in YOLOv2[20] and YOLOv3[18].

###  III-A YOLOv1

YOLO was first introduced in 2016 as a different approach to treat generic
object detection problems. The approach is based on training only one CNN for
the task of localization and the task of classification at the same time.
YOLOv1 is composed of 24 convolutional layers for feature extraction and 2
fully connected layers for generating the final output. The total architecture
is detailed in $Figure$ 1.

Figure 1: YOLO v1 Architecture

The input image is divided by YOLO into a grid of $S\times S$. Every grid cell
can be associated with only one object. The grid cell has also a fixed number
$B$ of boundary boxes to be predicted for this object. A confidence score is
calculated for every bounding box. As a result, YOLO calculates for every grid
cell a vector of class probabilities among the $C$ classes we are targeting.
YOLO calculates also for every bounding box of the cell a vector containing 5
parameters: $(x,y,w,h,box\\_confidence\\_score)$. Hence, the YOLO network
generates for every input image a tensor of the form:

$S\times S\times(B*5+C)$ (1)

where:

  * •

$S\times S$: corresponds to the number of grid cells.

  * •

$B$: corresponds to the number of bounding boxes.

  * •

$C$: corresponds to the number of targeted classes.

To train the YOLO network, we used a combination of three loss functions. The
first is the classification loss (loss of the conditional probabilities for
every class). The second is the localization loss (the position of the
estimated bounding box compared to the ground truth). The third is the
confidence loss (the box confidence score compared to the ground truth). To
measure the error between the predicted value and the ground truth, YOLO uses
sum-squared error as a metric. The expression of the loss is detailed in
$Equation$ 2.

$\begin{gathered}Loss=\lambda_{coord}\sum_{i=0}^{S^{2}}\sum_{j=0}^{B}{\rm
1\\!l}_{ij}^{obj}[(x_{i}-\hat{x}_{i})^{2}+(y_{i}-\hat{y}_{i})^{2}]\\\
+\lambda_{coord}\sum_{i=0}^{S^{2}}\sum_{j=0}^{B}{\rm
1\\!l}_{ij}^{obj}[(\sqrt{\omega_{i}}-\sqrt{\hat{\omega}_{i}})^{2}+(\sqrt{h_{i}}-\sqrt{\hat{h}_{i}})^{2}]\\\
+\sum_{i=0}^{S^{2}}\sum_{j=0}^{B}{\rm
1\\!l}_{ij}^{obj}(C_{i}-\hat{C}_{i})^{2}\\\
+\lambda_{noobj}\sum_{i=0}^{S^{2}}\sum_{j=0}^{B}{\rm
1\\!l}_{ij}^{noobj}(C_{i}-\hat{C}_{i})^{2}\\\ +\sum_{i=0}^{S^{2}}{\rm
1\\!l}_{i}^{obj}\sum_{c\in
classes}(p_{i}(c)-\hat{p}_{i}(c))^{2}\end{gathered}$ (2)

where:

  * •

$\lambda_{coord}$ is the weight for the bounding box coordinate prediction
loss. During the training, it is set to 5.

  * •

$\lambda_{noobj}$ is the weight for the confidence score prediction loss for
boxes that don’t include objects. During the training, It is set to 0.5.

  * •

${\rm 1\\!l}_{i}^{obj}$ indicates if the object appears in cell $i$

  * •

${\rm 1\\!l}_{i}^{obj}$ indicates that the bounding box of index $j$ is the
responsible for the prediction.

  * •

$x_{i}$ is the ground truth value of $x$, $\hat{x}_{i}$ is the predicted value
for $x$, The same for $y$, ( $\omega$ (width) and $h$ (height).

  * •

$C$ corresponds to the confidence score of the bounding box

  * •

$p_{i}(c)$ is the probability that the cell $i$ belong to the class $c$.
$\hat{p}$ denotes the predicted probability.

When it was introduced, YOLO surpassed other object detection algorithms in
terms of speed. Its mAP (mean Average Precision) was comparable or exceeded
the mAP of other state-of-the-art algorithms.

###  III-B YOLOv2

In YOLOv2[20], many improvements had been introduced to increase the accuracy
and decrease the processing time. Among them, we can note:

  * •

Use of Batch Normalization (BN). This technique was introduced in 2015 to
improve the convergence of the loss during the training. BN was added to all
convolutional layers in YOLO, which improved the mAP by 2%.

  * •

Replacing the input image size $224\times 224$ by $448\times 448$, which
improved the mAP by 4%.

  * •

The adoption of convolution with anchor boxes. The class prediction is moved
from the level of the grid cell to the level of the boundary box. This made a
small increase in the mAP by a margin of 0.3% but increased the recall from
81% to 88%. This enhanced the ability to detect all the objects existing in
the image and reduced the False Negative Rate.

  * •

The use of K-means clustering algorithm to select the best anchor box from the
training set of bounding boxes. Euclidian distance is replaced by the IoU
(Intersection Over Union) for the clustering.

  * •

The predictions are made on the offsets to the anchors. Instead of predicting
$(x,y,\omega,h,C)$, YOLOv2 predicts
$(t\textsubscript{x},t\textsubscript{y},t\textsubscript{w},t\textsubscript{h},t\textsubscript{C})$.
This makes the network converge better.

  * •

The use of fine-grained features. Similarly to the identity mapping in ResNet
[21], YOLOv2 concatenates low resolution features to high resolution features
to improve the ability to detect small objects. This improves the mAP by 1%.

  * •

The adoption of multi-scale training. Instead of using a fixed size for the
input image, YOLOv2 randomly selects an image size every 10 batches. This
improves the ability to predict well over different sizes of input images.

###  III-C YOLOv3

In April 2018, the YOLOv3[18] was introduced as an incremental improvement to
the previous versions. Among the improvements made, we can note:

  * •

The use of the multi-label classification. Instead of the mutual exclusive
labeling in the previous versions, YOLOv3 uses a logistic classifier to
estimate the likeliness of the object being of a specific label. The
classification loss is changed to use for every label the binary cross-entropy
loss instead of the general mean square loss used in the previous versions.

  * •

The use of another bounding box prediction. During the training, the
objectness score 1 in YOLOv3 is associated with the bounding box anchor that
best overlaps the ground truth object. Besides, if the IoU (Intersection Over
Union) between the bounding box anchor and the ground truth is less than a
threshold (0.7 in the implementation), it is ignored. In the end, YOLOV3
associates for every ground truth object one bounding box anchor.

  * •

The change of the output 3d tensor. In YOLOv3, the prediction is done for one
grid cell at 3 different scales and then concludes the final bounding box from
those scales. This was inspired by the feature pyramid networks. Hence, the
new dimension of the output 3d sensor is then:

$S\times S\times(3*(5+C))$ (3)

where:

    * –

$S\times S$: corresponds to the number of grid cells.

    * –

$B$: is omitted because only one bounding box anchor is kept at the end.

    * –

$C$: corresponds to the number of targeted classes.

  * •

The adoption of a new feature extractor (the darknet-53). It has 53 layers and
uses skip connection similarly to ResNet [21]. It uses both $3\times 3$ and
$1\times 1$ convolutions. It gave the state of the art accuracy but with
better speed and fewer computations.

##  IV DATASET

Figure 2: Sample images of our dataset, showing the different prayer
positions. The images on the top were captured specifically for this study,
while the images on the bottom were collected from the Internet.

We have built the Salat Postures Dataset from a set of images of students and
laboratory members performing different Salat positions using mobile phones,
from various view angles, especially captured for the purpose, complemented
with images (and video frames) of people in prayer collected from the
Internet. In fact, our ultimate scenario is that a person records himself
using a mobile phone, and then the video stream will be processed in real-time
to identify the different postures and finally provides postures flaws and
recommendation. For this reason, most of the collected images are based on
mobile phone data collection to serve the final scenario, added to this some
more general images. We have also considered images that show the whole body
of the person making prayer, as we impose that the person records the whole
body from any angle to perform the detection and classification.

The total number of collected images is 849 that we manually labeled into four
classes using rectangular bounding boxes. The four classes correspond to the
four main positions in Islamic prayer:

  * •

Qiyam: it is the act of standing, and this posture is related to when reading
the Quran, or after rising from bowing (Ruku). It typically spans over one or
more minutes. Sub-figure 5 of Figure 2 shows an example of a man in Qiyam
posture.

  * •

Ruku: it means bowing, and this posture happens after Qiam to make
invocations. It typically spans over a few seconds. Sub-figures 1 and 3 of
Figure 2 shows an example of a man in Ruku posture.

  * •

Sujud: It is the prostration and happens after raising from bowing for also
making invocation. It spans over several seconds. Sub-figure 4 of Figure 2
shows an example of a man in Sujud posture.

  * •

Julus: it is the position of sitting after finishing Sujud for making
invocation and closing Salat at its end. Sub-figures 2 and 6 of Figure 2 shows
an example of a man in Julus posture.

Figure 2 shows some sample images of the dataset, and Table I shows the number
of images and instances in the training and testing datasets. We have
subdivided the training data and test data to be 90% and 10%, respectively.

Table I: Number of images and instances in the training and testing datasets. | Training dataset | Testing dataset
---|---|---
Number of images | 764 | 85
Percentage | 90.0% | 10.0%
Instances of class "Standing" | 303 | 37
Instances of class "Bowing" | 210 | 27
Instances of class "Prostrating" | 174 | 11
Instances of class "Sitting" | 179 | 22

##  V EXPERIMENTAL EVALUATION

Figure 3: Average Precision (AP) at an IoU threshold of 0.5, for each input
size and each class.

For the experimental setup, we used a workstation powered by an Intel Core
i9-9900K (Octa-core) processor, with 64GB RAM, and an NVIDIA GeForce RTx2080T
(11 GB) GPU, running on Ubuntu 16.04 LTS.

We used the following metrics to assess the performance of the object
detection networks:

  * •

IoU: Intersection over Union, which measures the overlap between the predicted
bounding box and the ground-truth one.

  * •

mAP: mean average precision, which corresponds to the area under the precision
vs. recall curve, averaged over all classes. The mAP was measured for
different values of IoU (0.5, 0.6, 0.7, 0.8 and 0.9).

  * •

Inference time (in millisecond per image), which measures the inference
processing speed on test images.

  * •

TP (True Positives): number of objects correctly detected and classified.

  * •

FP (False Positives): number of objects detected but wrongly classified.

  * •

FN (False Negatives): number of non-detected objects.

Figure 4: Comparison of the trade-off between mAP and inference time for the
three different input sizes of YOLO v3, and for different values of IoU
threshold ranging from 0.5 to 0.9. Figure 5: Average intersection over Union
(IoU) for the three input sizes.  Figure 6: Number of true positives (TP),
false positives (FP), and false negatives (FN) for an IoU threshold of 0.5,
and an input size of 320x320.

Figure 3 shows the value of AP50 (Average Precision for an IoU threshold of
0.5), for each input size and each class. For all three networks, the
performance obtained for the detection of Ruku and Sujud classes are markedly
lower than for qiyam and sitting classes. This may be explained by the visual
resemblance between these two classes, which causes them to be conflated.

Figure 4 displays the trade-off between mAP and inference time for the 3 input
sizes, and for different values of IoU threshold ranging from 0.5 to 0.9. For
an IoU threshold of 0.5, passing from an input size of 320*320 to 416*416, the
mAP50 increases by 2%, while the inference time increases by 8%. Then, when
passing from an input size of 416*416 to 608*608, the mAP50 increases by 7%,
while the inference time jumps by 30%. For higher values of IoU threshold,
increasing the input size has little or negative effect on the mAP. This
indicates that higher input sizes do not enhance the precision of the bounding
boxes, which can also be seen in figure 5 that depicts the average IoU for
each input size. The value of mAP90 (at IoU=0.9) is close to zero, which shows
that YOLOv3 has difficulties in getting the bounding boxes to be perfectly
aligned with the object, as already noticed in the original YOLOv3 paper[18].

On another hand, the increase in inference time for higher input sizes is due
to the fact that a larger input size in the network generates a higher number
of parameters, which raises the number of operations necessary to obtain the
output.

Figure 6 depicts the number of true positives (TP), false positives (FP), and
false negatives (FN) for an IoU threshold of 0.5, and an input size of
320x320. We notice that the main issue with the ruku and sujud classes is the
relatively high number of false negatives compared to the two other classes.
This indicates that we should enrich the training dataset with a higher number
of diverse images, especially for these two classes.

##  VI CONCLUSIONS

In this paper, we conducted a performance evaluation study for detecting
Islamic prayer postures. The detection was performed by training the YOLOv3
network (with different configurations) on a dataset that was especially
collected for the purpose. The performance of the different network
configurations was assessed using several metrics (mAP, IoU, Inference time,
TP, FP, and FN). The mAP at IoU=0.5 ranges from 78% (for a network input size
of 320x320) to 85% (for an input size of 608x608), which is very encouraging
for a training dataset of only 764 images.

This study is a first step that will be further developed in future works, by
enlarging the training dataset, investigating other network architectures,
conducting thorough hyperparameter optimization, and integrating the object
detection model into a broader artificial intelligence assistive framework
that aims to guide Muslim worshippers to evaluate the correctness of their
postures during prayer.

## ACKNOWLEDGMENT

This work is supported by the Robotics and Internet-of-Things Lab of Prince
Sultan University.

## References

  * [1] O. Alobaid and K. Rasheed, “Prayer activity recognition using an accelerometer sensor,” in Proceedings on the International Conference on Artificial Intelligence (ICAI), pp. 271–277, The Steering Committee of The World Congress in Computer Science, Computer …, 2018.
  * [2] R. Al-Ghannam and H. Al-Dossari, “Prayer activity monitoring and recognition using acceleration features with mobile phone,” Arabian Journal for Science and Engineering, vol. 41, no. 12, pp. 4967–4979, 2016.
  * [3] M. F. Rabbi, N. Wahidah Arshad, K. H. Ghazali, R. Abdul Karim, M. Z. Ibrahim, and T. Sikandar, “Emg activity of leg muscles with knee pain during islamic prayer (salat),” in 2019 IEEE 15th International Colloquium on Signal Processing Its Applications (CSPA), pp. 213–216, March 2019\.
  * [4] M. Ramzan, A. Abid, H. U. Khan, S. M. Awan, A. Ismail, M. Ahmed, M. Ilyas, and A. Mahmood, “A review on state-of-the-art violence detection techniques,” IEEE Access, vol. 7, pp. 107560–107575, 2019.
  * [5] D. Ravì, C. Wong, B. Lo, and G. Yang, “A deep learning approach to on-node sensor data analytics for mobile or wearable devices,” IEEE Journal of Biomedical and Health Informatics, vol. 21, pp. 56–64, Jan 2017.
  * [6] W. Xu, Z. Miao, J. Yu, and Q. Ji, “Deep reinforcement learning for weak human activity localization,” IEEE Transactions on Image Processing, pp. 1–1, 2019.
  * [7] F. Ibrahim and S. A. Ahmad, “Assessment of upper body muscle activity during salat and stretching exercise: A pilot study,” in Proceedings of 2012 IEEE-EMBS International Conference on Biomedical and Health Informatics, pp. 412–415, Jan 2012.
  * [8] A. Gumaei, M. M. Hassan, A. Alelaiwi, and H. Alsalman, “A hybrid deep learning model for human activity recognition using multimodal body sensing data,” IEEE Access, vol. 7, pp. 99152–99160, 2019.
  * [9] B. Benjdira, T. Khursheed, A. Koubaa, A. Ammar, and K. Ouni, “Car Detection using Unmanned Aerial Vehicles: Comparison between Faster R-CNN and YOLOv3,” in 2019 1st International Conference on Unmanned Vehicle Systems-Oman (UVS), pp. 1–6, IEEE, 2019.
  * [10] A. Ammar, A. Koubaa, M. Ahmed, and A. Saad, “Aerial Images Processing for Car Detection using Convolutional Neural Networks: Comparison between Faster R-CNN and YoloV3,” arXiv pre-print 1910.07234, October 2019.
  * [11] B. Benjdira, Y. Bazi, A. Koubaa, and K. Ouni, “Unsupervised Domain Adaptation Using Generative Adversarial Networks for Semantic Segmentation of Aerial Images,” Remote Sensing, vol. 11, no. 11, 2019.
  * [12] B. Schoettle and M. Sivak, “A survey of public opinion about autonomous and self-driving vehicles in the us, the uk, and australia,” tech. rep., University of Michigan, Ann Arbor, Transportation Research Institute, 2014.
  * [13] M. H. El-Hoseiny and E. Shaban, “Muslim prayer actions recognition,” 2009 Second International Conference on Computer and Electrical Engineering, vol. 1, pp. 460–465, 2009.
  * [14] K. Eskaf, W. M. Aly, and A. Aly, “Aggregated activity recognition using smart devices,” in 2016 3rd International Conference on Soft Computing & Machine Intelligence (ISCMI), pp. 214–218, IEEE, 2016.
  * [15] M. Ali, M. Shafi, U. Farooq, et al., “Salat activity recognition using smartphone triaxial accelerometer,” in 2018 5th International Multi-Topic ICT Conference (IMTIC), pp. 1–7, IEEE, 2018.
  * [16] D. Ravi, C. Wong, B. Lo, and G. Yang, “Deep learning for human activity recognition: A resource efficient implementation on low-power devices,” in 2016 IEEE 13th International Conference on Wearable and Implantable Body Sensor Networks (BSN), pp. 71–76, June 2016.
  * [17] N. A. Jaafar, N. A. Ismail, and Y. A. Yusoff, “An investigation of motion tracking for solat movement with dual sensor approach,” 2015.
  * [18] J. Redmon and A. Farhadi, “Yolov3: An incremental improvement,” CoRR, vol. abs/1804.02767, 2018.
  * [19] J. Redmon, S. K. Divvala, R. B. Girshick, and A. Farhadi, “You Only Look Once: Unified, Real-Time Object Detection,” in 2016 IEEE Conference on Computer Vision and Pattern Recognition, CVPR 2016, Las Vegas, NV, USA, June 27-30, 2016, pp. 779–788, 2016.
  * [20] J. Redmon and A. Farhadi, “YOLO9000: better, faster, stronger,” in 2017 IEEE Conference on Computer Vision and Pattern Recognition, CVPR 2017, Honolulu, HI, USA, July 21-26, 2017, pp. 6517–6525, 2017.
  * [21] K. He, X. Zhang, S. Ren, and J. Sun, “Deep Residual Learning for Image Recognition,” Arxiv.Org, 2015.


* Hariharan and Girshick [2017] B. Hariharan, R. B. Girshick,  Low-shot visual recognition by shrinking and hallucinating features,  in: IEEE International Conference on Computer Vision, ICCV 2017, Venice, Italy, October 22-29, 2017, IEEE Computer Society, 2017, pp. 3037–3046. URL: https://doi.org/10.1109/ICCV.2017.328. doi:10.1109/ICCV.2017.328.
  * Antoniou et al. [2017] A. Antoniou, A. J. Storkey, H. Edwards,  Data augmentation generative adversarial networks,  CoRR abs/1711.04340 (2017). URL: http://arxiv.org/abs/1711.04340. arXiv:1711.04340.
  * Arjovsky et al. [2017] M. Arjovsky, S. Chintala, L. Bottou,  Wasserstein generative adversarial networks,  in: D. Precup, Y. W. Teh (Eds.), Proceedings of the 34th International Conference on Machine Learning, ICML 2017, Sydney, NSW, Australia, 6-11 August 2017, volume 70 of Proceedings of Machine Learning Research, PMLR, 2017, pp. 214–223. URL: http://proceedings.mlr.press/v70/arjovsky17a.html.
  * Zhang et al. [2018] R. Zhang, T. Che, Z. Ghahramani, Y. Bengio, Y. Song,  Metagan: An adversarial approach to few-shot learning,  in: S. Bengio, H. M. Wallach, H. Larochelle, K. Grauman, N. Cesa-Bianchi, R. Garnett (Eds.), Advances in Neural Information Processing Systems 31: Annual Conference on Neural Information Processing Systems 2018, NeurIPS 2018, December 3-8, 2018, Montréal, Canada, 2018, pp. 2371–2380. URL: https://proceedings.neurips.cc/paper/2018/hash/4e4e53aa080247bc31d0eb4e7aeb07a0-Abstract.html.


# Transparent gradient index lens for underwater sound based on phase advance

Theodore P. Martin<EMAIL_ADDRESS>U.S. Naval Research
Laboratory, Code 7160, Washington, DC 20375, USA  Christina J. Naify
Elizabeth A. Skerritt  Christopher N. Layman  Michael Nicholas  David C. Calvo
Gregory J. Orris  U.S. Naval Research Laboratory, Code 7160, Washington, DC
20375, USA  Daniel Torrent  José Sánchez-Dehesa  Department of Electronic
Engineering, Wave Phenomena Group, Universidad Politecnica de Valencia, C/
Camino de Vera s7n, E-46002 Valencia, Spain

###### Abstract

Spatial gradients in refractive index have been used extensively in acoustic
metamaterial applications to control wave propagation through phase delay.
This study reports the design and experimental realization of an acoustic
gradient index lens using a sonic crystal lattice that is impedance matched to
water over a broad bandwidth. In contrast to previous designs, the underlying
lattice features refractive indices that are lower than the water background,
which facilitates propagation control based on a phase advance as opposed to a
delay. The index gradient is achieved by varying the filling fraction of
hollow, air-filled aluminum tubes that individually exhibit a higher sound
speed than water and matched impedance. Acoustic focusing is observed over a
broad bandwidth of frequencies in the homogenization limit of the lattice,
with intensity magnifications in excess of 7 dB. An anisotropic lattice design
facilitates a flat-faceted geometry with low backscattering at 18 dB below the
incident sound pressure level. Three dimensional Rayleigh-Sommerfeld
integration that accounts for the anisotropic refraction is used to accurately
predict the experimentally measured focal patterns.

###### pacs:

43.20.+g,43.30+m,43.40.+s

Metamaterial lattices composed of sub-wavelength scattering components have
been used increasingly in recent years to control the propagation of both
electromagnetic and acoustic waves through two- and three-dimensions. One
primary goal of acoustic metamaterial design has been to achieve effective
fluid, or “metafluid,” material properties that minimize shear coupling and
propagation. The coupling between a fluid-born acoustic wave and a composite
elastic material at off-normal incidence results in mode-mixing of dilatation
and shear modes that complicates the interaction. For underwater applications
requiring strong acoustic coupling, compliant materials such as rubbers have
been a standard material used for coupling and encapsulation due to their
relatively close impedance match with water and low shear moduli. Traditional
compliant materials often have larger densities and/or lower dilatation moduli
than water, making them useful for applications requiring low relative sound
speed. For example, a sonic crystal lattice of rubber cylinders with a
gradient filling fraction was recently used to achieve an underwater
omnidirectional absorption coating Naify  _et al._ (2014), which requires a
radially decreasing sound speed. Li  _et al._ (2011); Climente  _et al._
(2012)

Metafluids with complementary material properties compared to rubbers, e.g.
high relative sound speeds with matched impedance, have been more difficult to
achieve in underwater applications because they require high stiffness but low
relative density. Sonic crystal lattices feature metafluid functionality in
the homogenization limit of the lattice, including shear decoupling and
broadband performance, and hence they have the potential to expand the range
of realizable metafluid material properties to include high sound speeds.
Thin-walled, hollow elastic shells have been recently proposed as high-sound
speed sonic crystal components in water. Martin  _et al._ (2012) By carefully
tuning the wall thickness of the shell, an impedance-matched condition can be
obtained assuming the shell material has higher impedance than the background
fluid. Additionally, by using shell elements which are individually impedance
matched, it has been shown that any phononic crystal configuration of the
elements will also be impedance matched in the homogenization limit of the
crystal lattice. Martin  _et al._ (2012) Given that local variations in the
element filling fraction create a variation in refractive index, acoustically
transparent devices with tunable wave-guiding capability should be achievable
using hollow shell elements. The concept of hollow shell elements has also
been extended to include structural components inside the shells to improve
control over the material properties.Titovich and Norris (2014)

Here we report the design and experimental demonstration of an acoustically
transparent gradient index (GRIN) lens that focuses sound by _advancing_ the
phase of waves propagating through an aqueous background. GRIN geometries have
been utilized for a broad range of metamaterial applications including
scattering reduction Zhang  _et al._ (2010); Yin  _et al._ (2012), wave
focusing Torrent and Sánchez-Dehesa (2007); Climente  _et al._ (2010); Martin
_et al._ (2010); Peng  _et al._ (2010); Sanchis  _et al._ (2010); Zigoneanu
_et al._ (2011); Lin  _et al._ (2012); Chang  _et al._ (2012); Hladky-Hennion
_et al._ (2013) and bending Ren  _et al._ (2010); Layman  _et al._ (2011); Li
_et al._ (2011); Climente  _et al._ (2012); Naify  _et al._ (2014). In
contrast to previous lens designs, we present a lens composed of impedance-
matched, hollow shell elements with sound speeds higher than water. The higher
sound speed enables the phase of propagating waves to be advanced beyond what
is possible in the background propagation medium, which is important to a
number of metamaterial applications (for a review, see Refs. 19; 20). The lens
is intended as a proof of concept to demonstrate that both transparency and
broadband wave-guiding functionality can be achieved simultaneously using a
lattice of non-compliant, low refractive index scattering elements.

The GRIN lens is constructed using cylindrical, air-filled aluminum tubes
arranged in a lattice that function as a broadband, uniform effective medium
at wavelengths $\lambda$ larger than the lattice spacing $a$ ($\lambda>4a$).
With a wall thickness of $1/20$th the diameter, each individual cylindrical
tube scatters sound as an impedance-matched effective fluid cylinder Martin
_et al._ (2012) resulting in negligible backscattering over the lattice as a
whole. An anisotropic lattice spacing is used to vary the filling fraction
similar to what was proposed by Lin, Tittmann, and Huang Lin  _et al._ (2012).
The anisotropic design simplifies the lens construction by using a single
cylindrical tube geometry and features a flat lens facet. While the objective
of many metafluid applications is to achieve anisotropy in the wave speed, the
transparency of the lattice results in an approximately isotropic scattering
configuration despite the underlying anisotropy of the design. Our design goal
is to tune the anisotropy to produce a convex refraction gradient that focuses
sound while simultaneously minimizing lens aberration.

Figure 1(a) shows the rectangular unit cell of the cylinder lattice. A
variable refractive index can be achieved by lengthening one of the lattice
constants $a_{y}$ compared to the lattice constant $a_{x}$ in the orthogonal
direction. Phase speeds $\bar{c}_{xx}$ and $\bar{c}_{yy}$ in the primary
Cartesian directions are shown in Fig. 1(b) as a function of the anisotropy
ratio $a_{y}/a_{x}$ for a fixed aluminum tube outer radius $R=0.4a_{x}$. The
double index on the phase speeds delineates components in an anisotropic
tensor. Phase speeds $\bar{c}$ with an over-bar are normalized to the water
background, which is assumed to have density $\rho_{b}=1000$ kg$/$m3 and sound
speed $c_{b}=1480$ m$/$s. The phase speeds are derived from the longitudinal
dispersion bands of the acoustic band structure calculated for each anisotropy
ratio. The anisotropy in phase speed can also be derived by considering
multiple scattering effects in the lattice. Torrent and Sánchez-Dehesa (2008)

Examples of the acoustic band structure are shown in Fig. 1(c,d). Blue and red
lines, with slopes corresponding to the phase speeds in Fig. 1(b), demonstrate
that the longitudinal bands have linear dispersion up to an upper cutoff
frequency $\omega a_{x}/2\pi c_{b}\simeq 0.3$. The cutoff frequency is
constrained by the wrapping of the longitudinal band at the Brillioun zone
boundary in the $\Gamma Y$ direction for the highest anisotropy ratio
considered. The longitudinal band in the $\Gamma M$ direction in Fig. 1(c) has
a slowness that deviates slightly from circularity due to the underlying cubic
symmetry. In addition to the longitudinal bands, narrow resonance bands
arising from the core-shell architecture of the unit cell are predicted at
various frequencies.Titovich and Norris (2014) While we do not observe
evidence of coupling to these bands in our measurements, the resonant bands
should not, in general, be overlooked depending on the lattice geometry.

Our lens design consists of a constant lattice spacing in the $x$-direction,
which produces a flat lens facet, while a prescribed lattice anisotropy in the
$y$-direction produces the convex GRIN profile [see Fig. 2(b)]. The primary
impact of increasing the anisotropy ratio is to decrease the effective sound
speed irrespective of propagation direction. As is evident in Fig. 1(b), the
difference between $\bar{c}_{xx}$ and $\bar{c}_{yy}$ is small even at
$a_{y}/a_{x}=2$, resulting in approximately isotropic acoustic transport over
the entire Brillioun zone. Therefore, an approximate, isotropic sound speed
$\bar{c}_{avg}=\left(\bar{c}_{xx}+\bar{c}_{yy}\right)/2$ was used to initially
guide the design of the GRIN geometry. A number of GRIN profiles that result
in convex lensing are commonly used in the literature; our design combines two
common refraction profiles in order to minimize lens aberration in the
“optical” (ray) transport limit,

$\displaystyle\bar{n}^{2}(y)$ $\displaystyle=$
$\displaystyle\frac{1}{\left[\bar{c}_{avg}(y)\right]^{2}}$ (1)
$\displaystyle=$
$\displaystyle\left(\beta\frac{\bar{c}_{0}^{2}}{1-\left(\alpha_{1}y\right)^{2}}+\left(1-\beta\right)\bar{c}_{0}^{2}\left(\cosh{\alpha_{2}y}\right)^{2}\right)^{-1}$
$\displaystyle\alpha_{1}$ $\displaystyle=$
$\displaystyle(2/h)\sqrt{1-\left(\bar{c}_{0}/\bar{c}_{h}\right)^{2}}$ (2)
$\displaystyle\alpha_{2}$ $\displaystyle=$
$\displaystyle(2/h)\operatorname{arccosh}\,(\bar{c}_{h}/\bar{c}_{0})\,$ (3)

where $h$ is the lens height in the gradient direction,
$\bar{c}_{0}=\bar{c}_{avg}(0)$ and $\bar{c}_{h}=\bar{c}_{avg}(h/2)$ are the
extremal sound speeds obtained from the chosen range of $a_{y}/a_{x}$, and
$\beta$ is a parameter that mixes the two GRIN profiles.

An iterative ray tracing routine based on the eikonal approximation Parazzoli
_et al._ (2006) was used to determine the optimal value of $\beta$ that
minimizes the lens aberration. The isotropic refraction gradient in Eq. (1)
prescribes $\bar{c}_{avg}(y)$ and hence the lattice geometry $a_{y}/a_{x}$ for
a given value of $\beta$; however, to optimize $\beta$ an eikonal
approximation that accounts for the anisotropic sound speeds in the lattice
was employed to improve the accuracy of the physical geometry. The
differential equation governing the eikonal function $\xi(\mathbf{x})$ in the
presence of anisotropy can be derived for acoustic metafluids in a similar
manner to the electromagnetic case, Parazzoli  _et al._ (2006)

$\displaystyle\bar{c}_{xx}^{2}\xi_{x}^{2}+\bar{c}_{yy}^{2}\xi_{y}^{2}-2\bar{c}_{xy}^{2}\xi_{x}\xi_{y}=1$
(4)

where $\xi_{[x,y]}=d\xi/d[x,y]$. The off-diagonal sound speed $\bar{c}_{xy}$
can be calculated from the band structure or by considering multiple
scattering effects. Torrent and Sánchez-Dehesa (2008) The off-diagonal term
was found to be negligible compared to the on-diagonal terms for each
anisotropy ratio in our lattice.

Figure 2(a) shows the ray paths derived using Eq. (4) for a GRIN lens
insonified by a plane wave after optimizing the parameter $\beta$. The lens
has thickness $10a_{x}$, height $h=27.1a_{x}$, and lattice spacing
$a_{x}=15.6$ mm; these parameters were also used to manufacture the physical
lens. The choice of $\beta=0.65$ produces a focal point at $x\simeq 1.2$ m
from the center of the lens with minimal aberration of the ray trajectories.
The $\beta$-optimization routine starts with the average sound speed profile
$\bar{c}_{avg}(y)$ as input and converts to the anisotropic phase speeds from
Fig. 1(b). The optimization routine also ensures that the gradient in the
lattice spacing $a_{y}(y)$ is properly discretized to be consistent with the
total height $h$ of the optimized GRIN profile.

A schematic of the upper half of the $\beta$-optimized lens design is shown in
Fig. 2(b). The full lens is symmetric about the $y$-axis. A photograph of the
assembled lens is shown in Fig. 2(c). The lens is constructed of identical 1 m
long hollow aluminum cylinders with air inside. Each cylindrical tube has a
diameter of 12.5 mm and nominal wall thickness of 0.625 mm. Rubber end caps
are adhered to both ends of each cylinder to prevent water infiltration.
Additionally, $\sim 2.5$ cm of foam was inserted into the tube ends to damp
the excitation of axial modes. The tubes are mounted on plexiglass end-plates
into which holes with the designed lattice pattern are drilled. Two additional
rows of cylinders with spacing $a_{y}=a_{x}$ are placed at the top and bottom
of the lens to help mitigate the sharp transition where the edge of the lens
($\bar{c}_{avg}=1.36$) meets the open water ($\bar{c}_{avg}=1$).

The lens was submerged at the center of a $6\times 6\times 4$ m3 water tank.
Acoustic waves were produced by a 0.1 m diameter spherical source located in
the plane that bisects the lens midpoint in the axial direction ($z$-axis). In
order to demonstrate the lens directionality, two in-plane source locations
were considered at $(x,y)=(-1.6,0)$ m and $(x,y)=(-1.6,0.43)$ m corresponding
to incident angles of $0^{\circ}$ and $15^{\circ}$ respectively from the lens
central axis ($x$-axis). Wave propagation was measured using hydrophones at a
sampling rate of 0.8 MHz. Hydrophones were mounted on a three-axis translation
positioning system to record measurements of the sound pressure levels,
$P(x,y)$ and $P_{0}(x,y)$, in the presence and absence of the lens
respectively. Measurements were taken at 10 mm increments. The sound pressure
level was measured at individual frequencies by averaging over a 10-cycle
pulse in the transmission region (as indicated in Figs. 3 and 4) and a
20-cycle pulse in the reflection region. The longer pulse cycle used in the
reflection region ensured overlap between incident and reflected waves so that
the measured total field could be properly compared with simulations. The
length of the incident pulses were short enough to isolate reflections from
the tank walls.

Measurements were performed in the $xy$-plane perpendicular to the cylindrical
axis of the lens over a broad range of frequencies $\omega a_{x}/2\pi
c_{b}<0.3$. Examples of the measured and numerically modeled pressure
intensity $(|P|/|P_{0}|)^{2}$ in the vicinity of the lens are shown in Figs. 3
and 4 for frequencies $\omega a_{x}/2\pi c_{b}=0.077$, 0.138, and 0.256. In
both the measured and predicted results, the pressure intensity is normalized
(at each position) to that of the source intensity in the absence of the lens.
Given the symmetry of the lens about the $y$-axis, measurements were obtained
over the bottom half of the scattering plane with a small overlap into the
upper half-plane to detail the on-axis forward scattering pattern. The regions
scanned in the measurement are also outlined in the upper panels of Figs. 3
and 4 for ease of comparison with the numerical predictions. Note that
although the measurements at $\omega a_{x}/2\pi c_{b}=0.077$ and 0.256 lie
within two of the resonance bands identified in Fig. 1, no obvious additional
resonant features are observed in the intensity maps that can be attributed to
these resonance bands.

The upper panels (a-c) of Figs. 3 and 4 show a prediction of the transmitted
and reflected pressure intensity calculated using two-dimensional multiple
scattering theory (2D-MST). Martin  _et al._ (2012); Torrent and Sánchez-
Dehesa (2008) The MST accounts for elastic scattering in the cylindrical tubes
and assumes insonification by a cylindrical monopole source located at
positions that match the experiment. The experimentally measured acoustic
intensities in panels (d-f) show good qualitative agreement with the MST-
predicted results. The reflected signal, $R=20\log_{10}[(P-P_{0})/P_{0})]$,
was estimated relative to the incident amplitude $P_{0}$ using the sound
pressure level $P$ measured in front of the incident face of the lens (regions
to the left of the lens in Figs. 3 and 4). The measured reflection was at or
below $13\%$ of the incident amplitude ($R\leq-18$ dB) over the range of
operational frequencies, which demonstrates significant transparency and is
commensurate with recent Fresnel lens designs in air. Molerón  _et al._
(2014); Li  _et al._ (2014) The measured reflection includes additional
constructive interference compared to the MST due to diffraction from the
finite aperture size in the axial direction.

In the forward direction a focusing peak is observed that strengthens in
intensity and moves out away from the lens with increasing frequency. The
frequency-dependence of the focusing peak is expected. At very low frequencies
the aperture is diffraction-limited; the wavelength is too long to adequately
resolve the index gradient and there is also significant diffraction around
the aperture. As the frequency increases, the transmission approaches the
geometric limit of ray-acoustics where a focusing peak would be located beyond
the focal point predicted by the ray tracing in Fig. 2(a). At the lowest
example frequency, $\omega a_{x}/2\pi c_{b}=0.077$ [panels (a,d)], the lens is
close to the diffraction limit with height-to-wavelength ratio $h/\lambda<2$.
As the frequency is further reduced the focusing peak intensity becomes
significantly suppressed. Therefore, an operational bandwidth can be
identified ranging between $\omega a_{x}/2\pi c_{b}\approx 0.05$–0.30, where
the lower cutoff is constrained by the diffraction limit and the upper cutoff
by the homogenization limit of the lattice.

It is apparent in Figs. 3 and 4 that the 2D-MST does not accurately predict
the precise location of the intensity maxima nor the small-scale, near-field
interference features in the forward direction. The discrepancy is due to the
assumption of infinite cylinders in the 2D theory compared to the finite size
of the lens aperture in the experiment. The diffractive corrections arising
from the finite size in the axial direction can be predicted using a three-
dimensional (3D) Rayleigh-Sommerfeld integral. The pressure at a position
$(x,y,z)$ in front of an acoustic aperture located at $x=0$ can be calculated
by integrating over the aperture, Born and Wolf (1999)

$\displaystyle P(x,y,z)$ $\displaystyle=$
$\displaystyle\frac{x}{i\lambda}\iint\limits_{A}\tilde{P}(y^{\prime},z^{\prime})\frac{e^{ikr}}{r^{2}}dy^{\prime}dz^{\prime}$
(5) $\displaystyle\tilde{P}(y^{\prime},z^{\prime})$ $\displaystyle=$
$\displaystyle\psi(y^{\prime},z^{\prime})e^{i\phi(y^{\prime},z^{\prime})}$ (6)
$\displaystyle r^{2}$ $\displaystyle=$ $\displaystyle
x^{2}+\left(y-y^{\prime}\right)^{2}+\left(z-z^{\prime}\right)^{2}\,$ (7)

where $\psi(y^{\prime},z^{\prime})$ and $\phi(y^{\prime},z^{\prime})$ are the
real-valued acoustic amplitude and phase of the complex pressure
$\tilde{P}(y^{\prime},z^{\prime})$ on the aperture facet, and $\lambda$ and
$k$ are the wavelength and wave number respectively. Here the integration
aperture is simplified by the flat facet of our lens design, with the
integration carried out over a simple rectangular area at the _forward_ face
of the lens.

The pressure on the aperture facet $\tilde{P}(y^{\prime},z^{\prime})$ can be
calculated using anisotropic ray-tracing under the eikonal approximation of
Eq. (4). A schematic in Figure 5(a) depicts the ray trajectories of wave-
fronts that emerge from the spherical source and traverse the forward facet of
the lens. The calculation iterates over a discrete set of rays that enter the
_rear_ facet (i.e. the incident face) at equally spaced intervals. The phase
function $\phi(y^{\prime},z^{\prime})$ on the forward facet was determined by
advancing the phase of each ray along the path lengths defined by
$\xi(\mathbf{x})$ assuming a constant phase at the source. The effective path
length is altered by the spatially-dependent local sound speed in the lens.
The phase difference $\Delta\phi=\phi(y^{\prime},z^{\prime})-\phi(0,0)$
relative to the phase at the lens central axis is shown in Fig. 5(b), where
positive or negative values represent a phase advance or delay respectively.
The calculation is performed for a source position on the lens central axis.
Given the symmetry of the lens only one quadrant is shown. Curve fits in Fig.
5(d) indicate that $\Delta\phi$ is approximately quadratic in both the $y$\-
and $z$-directions.

The amplitude function $\psi(y^{\prime},z^{\prime})$ is inversely related to
the square root of the local areal density of the rays. This was approximated
by calculating the average nearest-neighbor separation of rays intersecting
the forward facet; a relative amplitude
$\bar{\psi}(y^{\prime},z^{\prime})=\psi(y^{\prime},z^{\prime})/\psi(0,0)$ can
then be estimated from the inverse of this nearest-neighbor separation. Figure
5(c) shows the variation in $\bar{\psi}(y^{\prime},z^{\prime})$ after
normalizing to the predicted amplitude in the absence of the lens. It is clear
that the lens refraction has minimal impact on the amplitude ($<3\%$) and that
the diffractive correction is primarily influenced by the change in phase over
the aperture.

Following a similar procedure as was used by Gao _et. al._ , Gao  _et al._
(2012) Eq. (5) is evaluated numerically within the 2D plane of the measurement
using quadratic approximations to $\phi(y^{\prime},z^{\prime})$ and
$\psi(y^{\prime},z^{\prime})$. Intensity maps produced by the Rayleigh-
Sommerfeld integration are shown in Fig. 3(g-i) and demonstrate significantly
improved quantitative agreement with the experiment. As was the case for the
measurement and MST modeling, the Rayleigh-Sommerfeld-based pressure
intensities are normalized to the pressure intensity of the spherical source
in the absence of the lens. The 3D Rayleigh-Sommerfeld integration more
accurately predicts both the location of the focal positions and their
magnitude. The measured magnification in sound pressure level ranges between
4.5–7.3 dB over the operational bandwidth. The 2D-MST predicts magnifications
that are as much as $40\%$ lower than the measured values, whereas the
magnifications predicted by the Rayleigh-Sommerfeld integration agree to
within $10\%$. We emphasize that the 2D lattice only focuses along one axis
($y$-axis) resulting in lower magnifications than would be produced by
axisymmetric 3D designs. Molerón  _et al._ (2014); Li  _et al._ (2014) The
magnification can be increased by extending our 2D design to three dimensions
using tubes bent into a toroid configuration. Sanchis  _et al._ (2010); Chang
_et al._ (2012)

We now place our results in context with other metafluid design concepts.
Traditional underwater and ultrasound applications have achieved an impedance-
match by utilizing rubbers that have sound speeds less than water. While
rubbers have been used as components in transparent metafluid lattices, Naify
_et al._ (2014) there have been other recent advances in acoustic impedance-
matched designs. For example, negative-index complementary metamaterials have
been proposed to significantly enhance the transparency through aberrating
materials. Shen  _et al._ (2014) A class of metafluids based on the concept of
space-coiling have also featured prominently in the literature. Kock and
Harvey (1949); Liang and Li (2012); Xie  _et al._ (2013) Space-coiled Fresnel
lens designs have been reported that show significant focusing while utilizing
a thin aperture compared to the wavelength. Molerón  _et al._ (2014); Li  _et
al._ (2014) Although the strength of the space-coiling design is the ability
to significantly modify phase within a confined space, there is a drawback
that it can only delay the phase over a lengthened path; similar to rubbers in
water, the effective sound speed of these devices is less than the background
fluid. The space-coiling design has not yet been demonstrated in water, where
a higher viscosity and lack of rigid boundary conditions may impose additional
challenges.

Our design extends the reach of aqueous transparent metafluids to include the
option of _phase advance_ and decreased effective path length. While not
necessarily required for traditional lensing, there are metamaterial
applications that require higher sound speeds compared to the propagation
medium, Maldovan (2013); Kadic  _et al._ (2013); Ren  _et al._ (2010); Layman
_et al._ (2011); García-Meca  _et al._ (2014) the most prominent of which is
scattering reduction. The requirement of higher relative sound speeds features
prominently in scattering reduction designs based on both coordinate
transformation Cummer and Schurig (2007); Chen and Chan (2007); Cummer  _et
al._ (2008) and scattering cancellation. Guild  _et al._ (2012); Martin and
Orris (2012) The realization of a high-sound speed, transparent GRIN lens
represents an important proof-of-concept: it demonstrates that phase-advance
metafluids can be constructed with significant tunability in both sound speed
and impedance, including the option of matched impedance. Given that the
concept of core-shell lattice components has now been extended to also include
core modification, Titovich and Norris (2014) a wide range of GRIN and
anisotropic designs should be achievable in water based on metafluid lattices.

[Work is supported by the Office of Naval Research.]

Figure 1: (a) Schematic unit cell of the anisotropic lattice with
$a_{y}=2a_{x}$. (b) Phase speeds in the homogenization limit of the lattice
plotted as a function of anisotropy ratio $a_{y}/a_{x}$. Lines are fits to the
data based on Ref. [21]. (c) Acoustic band structure calculated using the
Finite Element Method (FEM) for an isotropic lattice with tube outer radius
$R=0.4a_{x}$. Blue lines plot the acoustic sound speed in the $\Gamma$X
direction. (d) Acoustic band structure calculated for an anisotropic lattice
with tube outer radius $R=0.4a_{x}$. Blue and red lines plot the acoustic
sound speed in the $\Gamma$X and $\Gamma$Y directions respectively. Insets
show the directions of the bands in $k$-space. Figure 2: (a) Ray paths through
the lens (blue lines), which originate from a plane wave at normal incidence,
are calculated using an eikonal approximation in the $xy$-plane; the red box
indicates the location of the lens. Only rays incident on the upper half-plane
of the lens are shown. (b) Schematic showing the locations of the lattice
sites in the upper half-plane of the lens; the resulting sound speed profile
is plotted to the right. The lattice sites are symmetric about $y=0$. (c)
Photograph of the assembled lens composed of 1 m long, hollow aluminum
cylinders arranged in the desired lattice pattern. Figure 3: Sound pressure
intensity in the vicinity of the lens at $\omega a_{x}/2\pi c_{b}=0.077$ (7.27
kHz), 0.138 (13.05 kHz), and 0.256 (24.24 kHz) for a source location on the
$x$-axis. (a-c) Pressure intensity calculated using 2D-MST. (d-f) Measured
pressure intensity. (g-i) Pressure intensity calculated using a 3D Rayleigh-
Sommerfeld approximation. Grey boxes show the position of the lens. Red
outlines indicate the experimentally mapped areas. The reduced observational
range of the reflected signal at 24.24 kHz is due to a limited overlap between
the incident and reflected signals in the time domain at high frequency.
Figure 4: Sound pressure intensity in the vicinity of the lens at $\omega
a_{x}/2\pi c_{b}=0.077$ (7.27 kHz), 0.138 (13.05 kHz), and 0.256 (24.24 kHz)
for a source positioned at a $15^{\circ}$ angle with respect to the $x$-axis.
(a-c) Pressure intensity calculated using 2D-MST. (d-f) Measured pressure
intensity. Grey boxes show the position of the lens. Red outlines indicate the
experimentally mapped areas. Figure 5: (a) Schematic showing the experimental
setup; blue lines indicate ray paths emerging from the spherical source and
Rayleigh-Sommerfeld integration vectors from the forward facet of the lens to
hydrophone locations in the measurement plane. (b) Change in phase
$\Delta\phi(y^{\prime},z^{\prime})$ (degrees) over the forward facet (due to
symmetry only the 1st quadrant is shown). Blue regions represent a phase
delay, red regions represent a phase advance. (c) Normalized amplitude
$\bar{\psi}(y^{\prime},z^{\prime})$ over the forward facet, where
normalization is with respect to the amplitude of the spherical wave in the
absence of the lens. (d) Red and blue lines show the phase change $\Delta\phi$
plotted as a function of position along the $y^{\prime}$\- and
$z^{\prime}$-axes respectively. Data points are calculated by advancing the
phase along ray paths based on an eikonal approximation, lines are quadratic
fits to the data.

## References

  * Naify  _et al._ (2014) C. J. Naify, T. P. Martin, C. N. Layman, M. Nicholas, A. L. Thangawng, D. C. Calvo, and G. J. Orris, Applied Physics Letters 104, 073505 (2014).
  * Li  _et al._ (2011) R.-Q. Li, X.-F. Zhu, B. Liang, Y. Li, X.-Y. Zou, and J.-C. Cheng, Applied Physics Letters 99, 193507 (2011).
  * Climente  _et al._ (2012) A. Climente, D. Torrent, and J. Sánchez-Dehesa, Applied Physics Letters 100, 144103 (2012).
  * Martin  _et al._ (2012) T. P. Martin, C. N. Layman, K. M. Moore, and G. J. Orris, Phys. Rev. B 85, 161103 (2012).
  * Titovich and Norris (2014) A. S. Titovich and A. N. Norris, The Journal of the Acoustical Society of America 136, 1601 (2014).
  * Zhang  _et al._ (2010) B. Zhang, T. Chan, and B.-I. Wu, Phys. Rev. Lett. 104, 233903 (2010).
  * Yin  _et al._ (2012) M. Yin, X. Yong Tian, H. Xue Han, and D. Chen Li, Applied Physics Letters 100, 124101 (2012).
  * Torrent and Sánchez-Dehesa (2007) D. Torrent and J. Sánchez-Dehesa, New Journal of Physics 9, 323 (2007).
  * Climente  _et al._ (2010) A. Climente, D. Torrent, and J. Sánchez-Dehesa, Applied Physics Letters 97, 104103 (2010).
  * Martin  _et al._ (2010) T. P. Martin, M. Nicholas, G. J. Orris, L.-W. Cai, D. Torrent, and J. Sánchez-Dehesa, Applied Physics Letters 97, 113503 (2010).
  * Peng  _et al._ (2010) S. Peng, Z. He, H. Jia, A. Zhang, C. Qiu, M. Ke, and Z. Liu, Applied Physics Letters 96, 263502 (2010).
  * Sanchis  _et al._ (2010) L. Sanchis, A. Yánez, P. L. Galindo, J. Pizarro, and J. M. Pastor, Applied Physics Letters 97, 054103 (2010).
  * Zigoneanu  _et al._ (2011) L. Zigoneanu, B.-I. Popa, and S. A. Cummer, Phys. Rev. B 84, 024305 (2011).
  * Lin  _et al._ (2012) S.-C. S. Lin, B. R. Tittmann, and T. J. Huang, Journal of Applied Physics 111, 123510 (2012).
  * Chang  _et al._ (2012) T. M. Chang, G. Dupont, S. Enoch, and S. Guenneau, New Journal of Physics 14, 035011 (2012).
  * Hladky-Hennion  _et al._ (2013) A.-C. Hladky-Hennion, J. O. Vasseur, G. Haw, C. Croënne, L. Haumesser, and A. N. Norris, Applied Physics Letters 102, 144103 (2013).
  * Ren  _et al._ (2010) C. Ren, Z. Xiang, and Z. Cen, Applied Physics Letters 97, 044101 (2010).
  * Layman  _et al._ (2011) C. N. Layman, T. P. Martin, K. M. Moore, D. C. Calvo, and G. J. Orris, Applied Physics Letters 99, 163503 (2011).
  * Maldovan (2013) M. Maldovan, Nature 503, 209 (2013).
  * Kadic  _et al._ (2013) M. Kadic, T. Bückmann, R. Schittny, and M. Wegener, Reports on Progress in Physics 76, 126501 (2013).
  * Torrent and Sánchez-Dehesa (2008) D. Torrent and J. Sánchez-Dehesa, New J. Phys. 10, 023004 (2008).
  * Parazzoli  _et al._ (2006) C. G. Parazzoli, B. E. C. Koltenbah, R. B. Greegor, T. A. Lam, and M. H. Tanielian, J. Opt. Soc. Am. B 23, 439 (2006).
  * Molerón  _et al._ (2014) M. Molerón, M. Serra-Garcia, and C. Daraio, Applied Physics Letters 105, 114109 (2014).
  * Li  _et al._ (2014) Y. Li, G. Yu, B. Liang, X. Zou, G. Li, S. Cheng, and J. Cheng, Sci. Rep. 4, 6830 (2014).
  * Born and Wolf (1999) M. Born and E. Wolf,  _Electromagnetic theory of propagation, interference and diffraction of light_ , seventh ed. (Cambridge University Press, Cambridge, 1999).
  * Gao  _et al._ (2012) Y. Gao, J. Liu, X. Zhang, Y. Wang, Y. Song, S. Liu, and Y. Zhang, Opt. Express 20, 1320 (2012).
  * Shen  _et al._ (2014) C. Shen, J. Xu, N. X. Fang, and Y. Jing, Phys. Rev. X 4, 041033 (2014).
  * Kock and Harvey (1949) W. E. Kock and F. K. Harvey, J. Acoust. Soc. Am. 21, 471 (1949).
  * Liang and Li (2012) Z. Liang and J. Li, Phys. Rev. Lett. 108, 114301 (2012).
  * Xie  _et al._ (2013) Y. Xie, A. Konneker, B.-I. Popa, and S. A. Cummer, Applied Physics Letters 103, 201906 (2013).
  * García-Meca  _et al._ (2014) C. García-Meca, S. Carloni, C. Barceló, G. Jannes, J. Sánchez-Dehesa, and A. Martínez, Phys. Rev. B 90, 024310 (2014).
  * Cummer and Schurig (2007) S. A. Cummer and D. Schurig, New Journal of Physics 9, 45 (2007).
  * Chen and Chan (2007) H. Chen and C. T. Chan, Applied Physics Letters 91, 183518 (2007).
  * Cummer  _et al._ (2008) S. A. Cummer, B.-I. Popa, D. Schurig, D. R. Smith, J. Pendry, M. Rahm, and A. Starr, Phys. Rev. Lett. 100, 024301 (2008).
  * Guild  _et al._ (2012) M. D. Guild, M. R. Haberman, and A. Alù, Phys. Rev. B 86, 104302 (2012).
  * Martin and Orris (2012) T. P. Martin and G. J. Orris, Applied Physics Letters 100, 033506 (2012).