\section{Introduction}
Cardiovascular diseases (CVDs) have been a major concern for human health for a long time and touted as the leading cause of death \cite{CDC_01}. World Health Organization (WHO) declared that CVDs were responsible for the death of almost 17.9 million people in 2016 - 31\% of all global deaths. 85\% of these deaths are due to heart attack and stroke \cite{Aa01}. CVDs can also cause permanent or temporary disabilities and reduce life-quality \cite{Ab02}. As a result, numerous research in this field with a view to preventing people from cardiac diseases have been carried out for the past decades. For in-time CVD diagnosis, it is vital to contentiously monitor HRV parameters that requires accurate estimation of IBI \cite{Harju2018MonitoringOH, Ahmed2010HeartRA}.

IBI is one of the most important parameters that can be extracted from ECG and photoplethysmography (PPG) signals. IBI in an ECG signal is the time interval between two consecutive beats as shown in Figure~\ref{fig:ibi}. It is generally measured in units of milliseconds and in normal heart function, each IBI value varies from beat to beat. This natural variation is known as HRV. However, certain cardiac conditions may cause the individual IBI values to become nearly constant, resulting in lower HRV values. This can happen, for example, during periods of exercise as the heart rate (HR) increases and the beats become more regular. Certain illnesses can cause the HR to increase and become uniform as well, such as when a subject is afflicted by an infection. In fact, IBI and HRV can be used as indicators of degraded cardiac system and can be early indicator of certain cardiac diseases \cite{Chessa2002RoleOH}.

Wearable sensors provide opportunity for continuous and convenient measurement of health parameters such as IBI. For the last decade, there has been immense development in the domain of wearable devices for physiological parameter monitoring, disease diagnosis and early prevention. Sleep apnea monitoring \cite{Ac03,Ad04}, cardiac anomalies or arrhythmia detection/classification \cite{Ae05,Af06,Qayyum2019ECGHC}, HRV estimation \cite{Ag07}, ECG monitoring \cite{Ah08}, and ophthalmic disease diagnosis \cite{Islam2019SourceAC} are some of the tasks worth mentioning.

\begin{figure}[!htbp]
  \centering
  \includegraphics[scale = 0.425]{figs/ibi.png}
  \caption{R-R interval or IBI in an ECG.}
  \label{fig:ibi}
\end{figure}

Although PPG sensors embedded in various wearables such as smart watches and smart rings enable convenient and continuous measurement of IBI, extracting IBI values from PPG is challenging as PPG is highly prone to noise. On the other hand, ECG signal that is encrypted with Gaussian out-of-band noise can be denoised easily using different types of finite impulse response (FIR) filters and empirical mode decomposition (EMD) techniques \cite{Chang2011GaussianNF}. However, an ECG signal that is encrypted with motion artifacts cannot be dealt easily as the frequency spectrum of some of these noises (1-10 Hz) \cite{Kher15} overlaps with that of the PQRST contents (0-50 Hz) \cite{Tereshchenko2015FrequencyCA,Kirst2011UsingDF}. Such situations emerge specially when the subject is in constant and/or periodic motion or performs exercise, where extracting IBI values from these motion encrypted ECG signals is critical. Even if IBI is directly estimated under this heavy motion situation, it will be nowhere near perfect. So, there is an unmet need for noise rejection techniques that can deal with high level of noise and motion artifact.

To overcome aforementioned unmet need, we propose a novel framework to determine IBI values from ECG signals totally buried in intense motion artifact noise. We have modified \textit{The One Hundred Layers Tiramisu} network \cite{Jgou29}, made it suitable for time-series denoising process, and used it for our case. The intention of this study is to eliminate motion artifact and make the R-peaks of an ECG signal more prominent in a way such that IBI values can be estimated from it as accurately as possible. Our modified network works as an autoencoder. An autoencoder is a well-known framework for making compact representation of a signal and taking it back to its original dimensions if required. Our methodology intends to squeeze the given ECG signal using autoencoder, suppress the noise, get back to the original dimensions, and represent the ECG signal with its peaks sufficiently visible to calculate IBI values. The autoencoder we have used is comprised of fully convolutional dense networks (FC DenseNet) - a Tiramisu model - which is familiar for image segmentation, image classification and more recently time domain analysis.

In the process of estimating IBI values from noisy ECGs, the contribution of this paper can be highlighted as follows.
\begin{itemize}
    \item[1.] An FC DenseNet based deep learning approach - derived from \textit{The One Hundred Layers Tiramisu model} \cite{Jgou29} - which performs as a robust autoencoder to suppress the noise and makes the R-peaks of the noisy ECG signal more notable. To the best of our knowledge, such a stacked tiramisu model is yet to be applied in beat detection from ECG with high intensity noise and IBI estimation as well.

    \item[2.] Our proposed methodology is capable of reckoning IBI values from ECG with noise level up to -30dB with some considerable amount of error. So far, we haven't noticed other works going beyond this range without the inclusion of high frequency or Gaussian white noise.

    \item[3.] In contrast to many concurrent works, our proposed approach does not require any sort of pre-processing or post-processing algorithms or tools. However, a simple peak picking algorithm has been used to facilitate the IBI calculation procedure.

\end{itemize}

The remaining of this paper is organized as follows. Related works are discussed in Section \ref{rltdwrks}. The dataset used in this work, noise signals, their characteristics, and noise addition protocols are discussed in Section \ref{dtst}, development of the proposed method is explained in Section \ref{prpsdmthd}, results, detailed comparison with state of the art work in Section \ref{rslts}, observations, drawbacks, and discussions are presented in Section \ref{disc}. In the end, we conclude with a brief conclusion in Section \ref{conc}.

\section{Related work}
\label{rltdwrks}

Research in the field of IBI estimation and ECG R-peak detection, which is the requirement of IBI estimation, from noisy ECG is high in number. Researchers have used a plethora of signal processing tools and techniques for accurate estimation of IBI in presence of motion artifacts. Mostly, they have focused on some publicly available ECG datasets like MIT-BIH Arrhythmia database or IEEE SP Cup Dataset (2015).

Removing high frequency noise or Gaussian noise from ECG and then estimating IBI from it can be achieved by applying  straightforward moving average \cite{Pandey09}, low-pass, high-pass or band-pass FIR filters \cite{Almal10,Rakshit11}, Empirical Mode Decomposition (EMD) technique \cite{Karagiannis12,Talbi13}, Pan-Tompkins peak detection \cite{Pan1985ARQ}, and wavelet transformation \cite{Sundararaj14}. However, the challenge starts when we tend to do the same in presence of heavy motion artifacts. Again, baseline wander and power line noise can be removed by using band-pass filter and notch filter. But, electrode motion artifact noise - which is caused by skin stretching and alters the impedance of the skin around the electrode - takes rigorous process to be encountered as its frequency spectrum (1 to 10 Hz) overlaps with that of the PQRST complex of ECG waveform \cite{Kher15}.

\textit{Rahel et al.} designed a study to evaluate the IBI signal qualities of a Holter device and a heart rate chest belt monitor while the subjects were at rest and were performing 5 different levels of activities like sitting and reading, doing household chores, walking, jogging, and training \cite{GilgenAmmann16}. But there is no indication of the specific amount of noise they dealt with. \textit{Shalini et al.} used Slope Sum Function and Teager-Kaiser Energy operator method for R-peak artifacts detection in cardiovascular and non-cardiovascular signals like Electroencephalogram (EEG), Electrooculogram (EOG), and Electromyogram (EMG) \cite{Rankawat22}. However, they too didn't mention the SNR level up to which their algorithm works perfectly.

\textit{Sonia et al.} discussed a new time delay estimation technique which is essentially derived from operational calculus, differential algebra, and non commutative algebra and helps to estimate RR interval from noisy ECG with SNR level up to 6dB with considerable number of false peak detection \cite{Rezk19}. \textit{Aygun et al.} proposed a technique to select the fiducial points from noisy ECG and PPG signal using shortest path algorithm which takes the advantage of time-continuity of heartbeats \cite{Aygun17}. With this, they accurately measured IBI up to -2dB SNR level which they later used for determining HRV \cite{Aygun18}. Here, the derivatives of the ECG or PPG signal have calculated using Savitzky-Golay method to detect the possible fiducial points like R-peaks, systolic peaks, points with maximum slope, and onset points. Since the signals are noisy, they ended up getting too many candidate points and most of them were false fiducial points if not all. With all these points, a graph was generated for each signal where the edges refer to IBIs, both true and false IBIs, and the vertices refer to morphological points, also both true and false. From this graph, it was clear that the starting point of one heartbeat is the end point of the last one and there is no disruption between them. Any candidate path that satisfies certain conditions derived from average HR can be considered as a true IBI and the two associated nodes/vertices as true morphological point for that specific signal. Next, each of the vertices were connected to some previous vertices which fell within a specific time window. If the time difference between the reference vertex and any neighboring vertex deviated the average IBI (which was calculated using average HR), the corresponding edge i.e. the 'time interval' was assigned zero weight. Otherwise, a numerical weight was assigned to it where the weight varied according to its difference with the average IBI. The intuition behind this weighting is the fact that the true IBIs should remain close to the average IBI calculated from the average HR. The weights of the vertices were also derived from the weights of corresponding edges. After all these weight distribution, the vertex with the minimum accumulated weight, is the chosen true vertex inside the time window. This is indeed leveraging the shortest path algorithm and enables detecting the IBIs that are closest to the average IBI. Finally, they combined all the shortest paths and formed an array for each of the morphological features of each ECG and PPG signal. An overview of this process is depicted in Figure~\ref{fig:aycya_flow}. Although both \cite{Aygun18} and \cite{Rezk19} have employed novel mathematical tools for IBI estimation without necessarily denoising the entire ECG signal, they haven't been able to go far in terms of SNR level.

\begin{figure*}[t]
\centering
\includegraphics[width=0.7\textwidth]{figs/Aycya_flow.png}
\caption{Overview of the work performed by \textit{Aygun et al.} \cite{Aygun18}.}
\label{fig:aycya_flow}
\end{figure*}

Many researchers have focused on R-peak detection from noisy ECG signals. As mentioned earlier, a robust peak detection algorithm can also lead to accurate IBI estimation. In addition to applying a band-pass filter for removing power line and baseline shifting noise, \textit{Pulavskyi et al.} performed a two-stage smoothing using 'boxcar' and 'parzen' kernel. This methodology allowed them to detect peaks with good precision and recall values up to -15dB, however, there was presence of white and pink noise which decreased the relative weight of electrode motion and muscle artifacts  \cite{Pulavskyi20}.

With advanced deep learning mechanisms, denoising tasks have gone further as they have been able to obtain higher accuracy and hardly require human supervision. To talk about a few deep learning based approaches, \textit{Ansari et al.} exploited a simple convolutional neural network (CNN) to differentiate between usable and unusable ECG segments. The usable segments have higher probability of QRS detection. They carried out their research up to 0dB with some error \cite{Ansari21}. \textit{Juho et al.}, on the other hand, utilized a bidirectional Long-Short Term Memory (LSTM) network to suppress the noise up to 0.1dB and detect peaks from it \cite{Laitala23}. A CNN encoder-decoder is used by \textit{Natasa et al.} to denoise the QRS complexes of the long term ECG signals acquired with their wearable armband device and used them later to calculate HR \cite{Reljin24}. They too were able to use the denoising process up to -17dB with inclusion of white and brown noise, but no muscle artifact. \textit{Sricharan et al.} modified the conventional U-Net \cite{Ronneberger25} algorithm for time-series data and used it along with distance transformation to determine the position of the ECG R-peaks for noise level up to 0dB \cite{Vijayarangan26}. The peak detection task was framed as a regression task in \cite{Vijayarangan26}. They obtained the Distance Transformation (DT) of all the ECG signals which disclose the distance of each point on the ECG from its nearest peak. It results in a zig-zag wave where the lowest points refer to the R-peaks of the corresponding ECG signal. The DT is of the same size as its input ECG. They have exploited the U-Net architecture with conventional encoder-decoder where the encoder performs downsampling with 8 layers of strided-convolutions, the decoder does upsampling in a similar way but in opposite direction. Both the encoder and decoder perform 1D convolution and there is a bottleneck between them to hold the minimum representation. Additionally, a residual inception block was placed at each layer which exploits skip connections to discard vanishing gradient problem and ensure fast convergence. The architecture for this work is illustrated in Figure~\ref{fig:Sri}.

\begin{figure*}[h]
\centering
\includegraphics[width=0.8\textwidth]{figs/Sricharan.png}
\caption{The U-Net architecture used by \textit{Sricharan et al.} \cite{Vijayarangan26}.}
\label{fig:Sri}
\end{figure*}

Likely, \textit{Lishen et al.} took the benefit of U-Net and DR-Net to perform a two-stage denoising process successfully \cite{Qiu27}. A generative adversarial network (GAN) based ECG synthesizer has been designed by \textit{Karol} to generate his own synthesized ECG dataset. Later, noise was added to the artificial ECG and a CNN autoencoder was trained  utilizing this dataset to achieve an MSE of 0.017 (direct from optimizer during training) \cite{Antczak2020AGA}. Finally, \textit{Brosnan et al.} \cite{Yuen28} demonstrated peak detection with a machine learning pipeline that consists of a Butterworth filter, two wavelet convolutional neural networks (WaveletCNNs) autoencoders, an optional QRS complex inverter, a Monte Carlo k-nearest neighbours (k-NN), and a convolutional long short-term memory (ConvLSTM) as shown on Figure~\ref{fig:pipeline}.

\begin{figure*}[h]
\centering
\includegraphics[width=0.9\textwidth]{figs/MLPipeline.png}
\caption{The machine learning pipeline developed by \textit{Brosnan et al.} \cite{Yuen28}.}
\label{fig:pipeline}
\end{figure*}

The noisy one-channel ECG signal is fed into the Butterworth filter with cut-off frequency set at 5Hz to perform the task of removing baseline shifting noise by attenuating the low frequency components. Next, the output is decomposed using symlet-4 wavelet to get the wavelet coefficients. The coefficients are fed into a CNN autoencoder where the noisy unnecessary coefficients are rejected but the significant parts retain. Thus it produces clean wavelet coefficients which are significant for further analysis. The CNN autoencoder that takes on the wavelet coefficients has a 2 layered encoder, a bottleneck and a 2 layered decoder and Leaky-ReLU activation layer at each hidden layer. A wavelet reconstruction is performed to get the coefficients back to original form i.e. single channel ECG. Now, a difference filter is deployed prior to the next wavelet autoencoder to refine the QRS complexes by suppressing the low frequency components. Next, the signal is moved to the second wavelet autoencoder which has same functionality as the previous one. Both of the two wavelet autoencoder work as bandpass filters to remove the unnecessary components. To deal with the inverted peaks, an inverter was placed right after the second autoencoder and it just inverted the already inverted peaks. Then Monte Carlo k-NN method was utilized to normalize the ECGs by scaling the beats to unit mili-volt. Finally, the ConvLSTM network was applied to detect the peak points as it has pattern matching ability to detect any specific characteristics from a time-series signal.They used the MIT-BIH, the European ST-T and the Long Term ST database Noise Stress Test databases, went up to -6dB and recorded some good result. All these deep learning based methods, however, failed to leverage different deep learning algorithms' power to improve the accuracy enough for denoising or IBI estimation and explore intense noise scenarios.

\section{Datasets}
\label{dtst}
We have opted to use clean ECG recordings from publicly available database like MIT-BIH Arrhythmia Database v1.0.0 \cite{Moody2001TheIO,Goldberger2000PhysioBankPA}, European ST-T Database from PhysioNet \cite{Taddei1992TheES} and IEEE Signal Processing Cup 2015 \cite{Zhang2015TROIKAAG}.

The sampling frequency of MIT-BIH Arrhythmia Database is 360 Hz i.e. 1 sample in every 2.78 ms. For ease of calculation and noise addition, we upsampled both European ST-T Database and IEEE SP Cup 2015 dataset from 250 Hz and 125 Hz to 360 Hz respectively.

We then added electrode motion (EM) noise, motion artifact (MA) noise and baseline wandering (BW) noise from  MIT-BIH Noise Stress Test (NST) Database v1.0.0 \cite{Moody1992TheMN} and generated noisy data for different SNR levels starting from 36dB to -36dB with a decrement of 6dB. So, for each clean ECG signal, we generated 13 variants with different SNR levels. Needless to mention, noise have been added to the signals thoroughly instead of adding them between consecutive rest periods. For this noise addition task we followed the rules given in NST Database-
\begin{equation*}
\textit{ECG}\textsubscript{noisy} = \textit{ECG}\textsubscript{clean} + [a_1]*\textit{EM} + [a_2]*\textit{MA} + [a_3]*\textit{BW}
\end{equation*}

where, $\textit{ECG}\textsubscript{clean}$ refers to clean ECG signal and $a_1, a_2, a_3$ gain values are adjusted according to the expected SNR value. Since we are more interested in estimating IBI values from a signal totally buried under muscle artifact and electrode motion, these two noises have been assigned higher weights meaning $a_2 > a_1 > a_3$.

\begin{figure}[!htbp]
\centering
\includegraphics[scale = 0.28]{figs/noise_new.png}
\caption{A clean ECG segment and its noisy variants. EM, MA, BW noises have been added to clean signals and the numbers on the left of each box refer to corresponding SNR values. IBIs estimated from raw and noisy ECG signal are highly inaccurate compared to the true IBIs}
\label{fig:noise}
\end{figure}

\begin{table*}[b]
\begin{center}
\caption{List of used ECG recordings from different database}
\setlength{\arrayrulewidth}{0.4mm}
\setlength{\tabcolsep}{18pt}
\resizebox{0.7\textwidth}{!}{%
\begin{tabular}{|cc|c|c|}
\hline
\multicolumn{2}{|c|}{MIT-BIH Arrhythmia Database} &
European ST-T Database &
IEEE SP Cup 2015 \\
\hline
\multicolumn{2}{|c|}{Train data} & Train data & Train data \\ \hline
100 & 102 & e0103 & 1 \\
 & & e0104 & 5 \\ \hline
\multicolumn{2}{|c|}{Test data} & Test data & Test data \\ \hline
101 & 209 & e0105 & 2 \\
103 & 213 & e0106 & 3 \\
109 & 219 & e0107 & 4 \\
112 & 220 & e0113 & 6 \\
113 & 223 & e0114 & 7 \\
115 & 228 & e0115 & 8 \\
116 & 230 & e0118 & 9 \\
122 & 231 & e0122 & 10 \\
123 & 234 & e0123 & 11 \\
&  & e0126 & 12 \\
&  & e0127 &  \\
&  & e0129 &  \\

\hline
\end{tabular}}
\label{tab1}
\end{center}
\end{table*}

In this work, we used first 409600 data points of each ECG signal since the noise signals have limited data points. Figure~\ref{fig:noise} shows a clean segment of ECG data from MIT-BIH Arrhythmia database (recording \# 101) and its corresponding noisy variants starting from 6db to -36dB. As the noise increases, specifically after 0 dB, the R-peaks start to heavily get buried in noise and detecting them with conventional signal processing techniques would be very inaccurate and challenging if even possible. These artifacts enforce redundant tests, additional costs, and specialists' intervention sometimes. Table~\ref{tab1} lists the recordings from different dataset used for training and testing. Usually, autoencoder type models require very small training data \cite{Vijayarangan26} and our case is no different.

\section{Proposed Methodology}
\label{prpsdmthd}

Our working principle starts with a fully convolutional-dense deep neural network (tiramisu model) \cite{Jgou29}. Noisy signals are fed to this model which has multiple convolution and pooling layers intuitively responsible for compressing the signal to a smaller representation, eliminating the noise, retaining the most patterned and prominent features, and making the peaks more visible. Later, we identify the peaks with a simple peak picking algorithm and then calculate IBI values from there. The overall procedure of this work is illustrated in Figure~\ref{fig:flow}.

\begin{figure*}[!htbp]
\centering
\includegraphics[width=0.9\textwidth]{figs/flow.png}
\caption{An overview of the proposed IBI estimation framework.}
\label{fig:flow}
\end{figure*}

\subsection{Tiramisu Model}
It is challenging to calculate IBI values from noisy ECG signals. As a solution, we designed a robust tiramisu model which consists of stacked dense blocks (DB) and numerous feedback loops. High number of dense blocks and feedback loops eliminate the requirement of pre-processing as the model itself is capable of removing high frequency or Gaussian white noise. While the convolution and pooling layers perform the job of squeezing out the noise, the added skip connections between layers allow the output of each dense block and up convolution to encode finer details from actual signal and later layers, solve the gradient vanishing problem, and provide improved results than without them. In contrast to the standard U-Net architecture, the dense blocks in the tiramisu model have skip connections within themselves which is a fact that eliminates the possibility of gradient vanishing within the contraction path.

The overall architecture of the proposed model is illustrated in Figure~\ref{fig:tiramisu}. As depicted, the model incorporates a contraction path, a bottleneck, an expanding path, and certain skip connections. The skip connections here help the expanding path (or upsampling path) to get back information from the contraction path, reuse them and thus abolish gradient vanishing. The main theme of this entire model is to take advantage of feature reuse by taking the already complicated DenseNet a step ahead.

\begin{figure}[t]
\centering
\includegraphics[scale=0.18]{figs/tiramisu.png}
\caption{Big picture view of the proposed model. It consists of a contraction and an expanding path. Both paths have several convolution blocks, dense blocks, transition downs, transition ups, concatenations and skip connections. The dense block between these two paths is the bottleneck.}
\label{fig:tiramisu}
\end{figure}

The contraction path (or downsampling path) starts with a convolution with kernel size = 3. Then comes a series of consecutive dense blocks and transition down (TD) blocks. Prior to the bottleneck, 3 dense blocks have been placed on the contraction path with each having 4, 5, and 7 layers respectively. Right after each dense block, a TD block has been attached. Note that, the TD blocks perform the pooling operation and cause some information loss and resolution reduction along the downsampling path. In between, there are concatenation blocks to create feedback loops from previous layers.

The expanding path is almost a mirror symmetry of the contraction path except it has lesser feedback loops and the TD blocks are replaced by transition up (TU) blocks along the path. These TU blocks does transposed convolution operation and upsample the previous mappings. The upsampled mappings are then concatenated to the mappings coming from the downsampling path via skip connections and form a new input for the next dense blocks along expanding path. The last dense block in this trajectory sums up the information from all the previous dense blocks. Unlike the contraction path, we do not concatenate the input of a dense block with its output since it does not make any significant difference. The final convolution in this path is done with widow size = 3. The dense block that lies between these two paths is the bottleneck.

The dense blocks in the model have different number of layers which perform one dimensional convolution with kernel size = 3, \textit{ReLU} activation function, Batch Normalization, and drop out value of 0.2. The batch normalization function is given by-

\begin{equation}
    BN(x) = \frac{x-\mu_x}{\sigma_x}
\end{equation}

where, $\mu_x$ is the mean and $\sigma_x$ is the standard deviation of input matrix $x$.

The TD layers also use one dimensional convolution with kernel size = 1, \textit{ReLU} activation, Batch Normalization, drop out value of 0.2 and Max Pooling with pool size = 2. However, the TU layers perform transposed convolution with kernel size = 3 and stride = 2.

On the final output layer we have employed \textit{tanh} activation since we desire to have our output between [-1 1] and gain non-linearity for getting advantage during derivative calculation. The \textit{tanh} function can be written as-

\begin{equation}
    tanh(x) = \frac{e^{x}-e^{-x}}{e^{x}+e^{-x}}
\end{equation}

As mentioned earlier, we have 3 dense blocks on either path with each block having 4, 5, and 7 layers respectively. Also, we have 3 TUs or TDs and 1 input or output block on either side. Since each layer has one convolution layer, each TU or TD also perform one convolution and each input/output block has one convolution layer, we have 2*(4+5+7+3+1) = 40 convolution layers on these two paths combined. Additionally, the bottle neck dense block has 10 layers i.e. 10 convolution layers. So, a total of 50 convolutions are done in this entire model. All the convolutions (except last one) are equipped with \textit{ReLU} activation, L2 regularization with factor of 0.01 and zero padding with a view to redeeming the initial shape. Finally, we have employed the Adam optimizer with default learning rate i.e. 0.001 and mean square error (MSE) based loss function-
\begin{equation}
    MSE = \frac{1}{N}\sum_{i=1}^{N}(y_i - y_i^p)^2
\end{equation}

where, $y_i$ and $y_i^p$ refer to actual target value and predicted value respectively and $N$ is the number of data points.

Table \ref{tab2} shows the full architecture of our tiramisu model and the constituents of layer, TD, and TU blocks.

\begin{table}[!htbp]
    \caption{Full tiramisu model architecture and contents of layer, TD and TU blocks.}
    \begin{minipage}{.32\linewidth}
      \centering
        \begin{tabular}{||p{4cm}||}
        \hline
           \begin{center}\textbf{Full tiramisu model architecture}\end{center}\\
         \hline\hline
         Input\\
         \hline
         Convolution: kernel size = 3\\
         \hline
         DB (4 layers) + TD \\
         \hline
         DB (5 layers) + TD\\
         \hline
         DB (7 layers) + TD\\
         \hline
         Bottleneck: DB (10 layers)\\
         \hline
         DB (7 layers) + TU\\
         \hline
         DB (5 layers) + TU\\
         \hline
         DB (4 layers) + TU\\
         \hline
         Convolution: kernel size = 3\\
         \hline
         \textit{tanh} activation\\
         \hline
         Output\\
         \hline
        \end{tabular}
    \end{minipage}%
    \begin{minipage}{.4\linewidth}
      \centering
        \begin{tabular}{||p{3.2cm}||p{3.2cm}||p{3cm}||}
        \hline
            \begin{center}\textbf{Layer block}\end{center} & \begin{center}\textbf{Transition Down (TD) block}\end{center}&\begin{center}\textbf{Transition Up (TU) block}\end{center}\\
        \hline\hline
        Batch Normalization & Batch Normalization & Trans. convolution: kernel size = 3\\
        \hline
        \textit{ReLU} activation & \textit{ReLU} activation & l2 regularizer (0.01)\\
        \hline
        l2 regularizer (0.01) & l2 regularizer (0.01) & \textit{stride} = 2\\
        \hline
        Convolution: kernel size = 3 & Convolution: kernel size = 1 &\\
        \hline
        Dropout $p$ = 0.2 & Dropout $p$ = 0.2 &\\
        \hline
        &Max Pooling: size = 2&\\
        \hline
        \end{tabular}
        \label{tab2}
    \end{minipage}
\end{table}

\subsection{Peak Detection}
For peak detection, we have employed a general algorithm \cite{peak}. Usually, resting HR varies from 60-100 bpm \cite{American_Health}. But for athletes, it could be as low as 40 bpm during rest periods and as high as 200 bpm during work out sessions \cite{Athletes}. Our sampling frequency is 360 Hz. It ensures that there must be at least one beat in every 108 to 540 data points. So we have used an algorithm that takes on a given signal and finds out all the local maxima which must be at least 108 data points apart from each other. In this way, we may end up getting some false peaks because, when the noise is intense and a motion artifact-induced peak is higher than its neighboring R-peak, the tiramisu model fails to detect the R-peak and removes it while the motion artifact remains. Subsequently, the peak detection algorithm detects this motion artifact as an R-peak causing some error.

\subsection{IBI Calculation}
Under any circumstance, to estimate IBI precisely, the position of the beats/peaks need to be detected accurately. After denoising the signal using the proposed AE-based model and detecting the R-peaks, in this step, the IBI value using the identified R-peaks can be extracted with Equation \ref{eq1}.
\begin{equation} \label{eq1}
    IBI_{i} = t_i-t_{i-1}
\end{equation}

\noindent where, $t_i$ is the occurrence time of the $i^{th}$ R-peak and IBI is extracted for all peaks in an ECG signal.

\section{Results}
\label{rslts}

In this study, we aim to make the R-peaks of noisy ECG more prominent to estimate IBI from it as accurately as possible. We evaluate the performance of our approach by testing it on different levels of noisy ECG, which have been generated by adding certain level of noise, explained in Section \ref{dtst}. In this section, we define our evaluation metrics and measure the performance of our model in terms of estimated IBI values extracted from ECGs with different SNR levels. Lastly, we will compare our approach against a few state-of-the-art methods designed for IBI calculation and/or peak detection from noisy ECG.

For further explanation we refer to Figure~\ref{fig:demo}. The clean ECG in this figure is the first 40960 data points of recording 234 from MIT-BIH Arrhythmia Database. Three different noises have been added to the clean ECG to get a noisy signal of -30dB SNR. -30dB refers to a very high amount of noise and implies that the power of noise is $10^3$ times that of actual signal.  Without any pre-processing, this noisy ECG is then fed to the tiramisu model to make the R-peaks more visible. Next, our peak detection algorithm has been applied to the output of the model. We estimated IBI values from the beats and in the final box, we have compared our IBI values with that of the ground truth and recorded a root mean squared error of 3.44 ms.

\begin{figure*}[t]
\centering
\includegraphics[width=\textwidth, height=9cm]{figs/Demo.png}
\caption{\textbf{\textit{Clean ECG:}} An ECG segment from MIT-BIH Arrhythmia Database (recording \# 234). \textbf{\textit{Noisy ECG:}} Muscle artifact, baseline wander and electrode motion noise have been added to the clean signal. Resulting SNR of the ECG is -30dB. \textbf{\textit{Model output:}} The noisy ECG is fed to the tiramisu model and the output has the R-peaks more prominent. \textbf{\textit{Peak Detection:}} Peak detection algorithm has been applied to the output. \textbf{\textit{IBI Comparison:}} IBI values estimated from the prominent peaks and compared against the ground truth IBI values. The RMSE of IBI values in this case is only 3.44 ms.}
\label{fig:demo}
\end{figure*}

\begin{figure}[t]
     \centering
     \begin{minipage}[!htbp]{0.24\textwidth}
         \centering
         \includegraphics[width=\textwidth]{figs/BA230.png}\\
         \subcaption{recording \# 230}
         \label{BA1}
     \end{minipage}%
     \hfill
     \begin{minipage}[!htbp]{0.24\textwidth}
         \centering
         \includegraphics[width=\textwidth]{figs/BA220.png}\\
         \subcaption{recording \# 220}
         \label{BA2}
     \end{minipage}%
     \hfill
     \begin{minipage}[!htbp]{0.24\textwidth}
         \centering
         \includegraphics[width=\textwidth]{figs/BA115.png}
         \subcaption{recording \# 115}
         \label{BA3}
     \end{minipage}%
     \hfill
     \begin{minipage}[!htbp]{0.24\textwidth}
         \centering
         \includegraphics[width=\textwidth]{figs/BA213.png}
         \subcaption{recording \# 213}
         \label{BA4}
     \end{minipage}%

     \begin{minipage}[!htbp]{0.24\textwidth}
         \centering
         \includegraphics[width=\textwidth]{figs/PCP230.png}
         \subcaption{recording \# 230\\\textit{r} = 0.9595}
         \label{C1}
     \end{minipage}%
     \hfill
     \begin{minipage}[!htbp]{0.24\textwidth}
         \centering
         \includegraphics[width=\textwidth]{figs/PCP220.png}
         \subcaption{recording \# 220\\\textit{r} = 0.9273}
         \label{C2}
     \end{minipage}%
     \hfill
     \begin{minipage}[!htbp]{0.24\textwidth}
         \centering
         \includegraphics[width=\textwidth]{figs/PCP115.png}
         \subcaption{recording \# 115\\\textit{r} = 0.9895}
         \label{C3}
     \end{minipage}%
     \hfill
     \begin{minipage}[!htbp]{0.24\textwidth}
         \centering
         \includegraphics[width=\textwidth]{figs/PCP213.png}
         \subcaption{recording \# 213\\\textit{r} = 0.9006}
         \label{C4}
     \end{minipage}%

        \caption{Bland-Altman plot with LOAs are depicted in Figs. \ref{BA1} to \ref{BA4} for different recordings. Also, correlation plot for same recordings are shown with Pearson correlation coefficient in Figs. \ref{C1} to \ref{C4}. For some ECGs, the proposed method achieves r values like 0.9895 though nominal values like 0.9 also exist.}
        \label{fig:seven graphs}
\end{figure}

\begin{figure}[!htbp]
     \centering
     \begin{minipage}[!htbp]{0.25\textwidth}
         \centering
         \includegraphics[width=\textwidth]{figs/BR_BA_230.png}\\
         \subcaption{\textit{Brosnan et al.}}
         \label{BA5}
     \end{minipage}%
     \hfill
     \begin{minipage}[!htbp]{0.25\textwidth}
         \centering
         \includegraphics[width=\textwidth]{figs/SR_BA_230.png}\\
         \subcaption{\textit{Sricharan et al.}}
         \label{BA6}
     \end{minipage}%
     \hfill
     \begin{minipage}[!htbp]{0.25\textwidth}
         \centering
         \includegraphics[width=\textwidth]{figs/AY_BA_230.png}
         \subcaption{\textit{Aygun et al.}}
         \label{BA7}
     \end{minipage}%

     \begin{minipage}[!htbp]{0.25\textwidth}
         \centering
         \includegraphics[width=\textwidth]{figs/BR_PCP_230.png}
         \subcaption{\textit{Brosnan et al.}\\\textit{r} = 0.8883}
         \label{C5}
     \end{minipage}%
     \hfill
     \begin{minipage}[!htbp]{0.25\textwidth}
         \centering
         \includegraphics[width=\textwidth]{figs/SR_PCP_230.png}
         \subcaption{\textit{Sricharan et al.}\\\textit{r} = 0.8663}
         \label{C6}
     \end{minipage}%
     \hfill
     \begin{minipage}[!htbp]{0.25\textwidth}
         \centering
         \includegraphics[width=\textwidth]{figs/AY_PCP_230.png}
         \subcaption{\textit{Aygun et al.}\\\textit{r} = 0.7897}
         \label{C7}
     \end{minipage}%

    \caption{Bland-Altman plot for IBI estimation on recording \# 230 using methods of \textit{Brosnan et al.}, \textit{Sricharan et al.}, and \textit{Aygun et al.} are depicted in Figs. \ref{BA5} to \ref{BA7} respectively. Corresponding correlation plots are also shown along with correlation coefficients in Figs. \ref{C5} to \ref{C7}.}
    \label{fig:six graphs}
\end{figure}

Prior studies in the field of ECG denoising and beat detection exercised \textit{SNR\textsubscript{imp}}, \textit{Sensitivity} or \textit{F\textsubscript{1} score} to evaluate their performance regarding detecting correct peaks. However, herein, we are interested in estimating IBI from noisy signals. Therefore, we selected our performance metrics based on that objective and leveraged \textit{Root Mean Squared Error (RMSE)} and \textit{\% error} of IBI for performance evaluation. We define those metrics as follows

\begin{align}
\operatorname { RMSE_{\textit{IBI}} } = \sqrt { \frac { \sum _ { i = 1 } ^ { n } ( P _ { i } - O _ { i } ) ^ { 2 } } { n } }
\end{align}

\begin{align}
\% error = \frac{1}{n}\sum _ { i = 1 } ^ { n }\frac { | O _ { i } - P _ { i } | } {O _ { i }  } \times 100\%
\end{align}
where, \textit{P\textsubscript{i}} refers to predicted IBI, \textit{O\textsubscript{i}} denotes observed IBI (groundtruth), and \textit{n} indicates number of IBIs in one segment of ECG.

The Bland-Altman plots in Figs. \ref{BA1} to \ref{BA4} show comparison of the IBIs estimated from noisy ECGs of SNR -30dB against their true IBIs. The ECGs considered in these plots are recordings \# 230, 220, 115 and 213. The limits of agreements (LOAs) - which work as boundary for 95\% of data and a good measurement of accuracy - are located at [49.92, -50.85] ms, [73.74, -72.44] ms, [20.14, -20.29] ms, and [33.89, -33.52] ms respectively. Also, the Pearson correlation plots are illustrated in \crefrange{C1}{C4} for the same set of recordings. The high correlation coefficients (\textit{r}) - 0.96, 0.93, 0.99, and 0.9 - emphasize that the estimated IBIs are very coherent to the true IBIs.

In comparison, similar plots are also shown in \crefrange{BA5}{BA7} and \crefrange{C5}{C7} highlighting the performance of \textit{Brosnan et al.}, \textit{Sricharan et al.}, and \textit{Aygun et al.} respectively on recording \# 230. The LOAs are at [90.25, -94.40] ms, [88.87, -88.38] ms, and [129.829, -134.76] ms respectively. So, our method is ensuring [40.33, -43.55] ms, [38.95, -37.53] ms, and [79.91, -83.91] ms of improvements respectively. As stated in corresponding papers, these methods also perform well under low noise situations. Correlation coefficients (\textit{r}) for these methods on this specific recording are 0.88, 0.87, and 0.79 respectively whereas our approach obtains 0.96 correlation for the same setting. For other recordings of same dataset, none of these 3 methods achieve correlation above 0.9 at an SNR value of -30dB. Performance of all methods (including ours) decline for recordings of IEEE SP Cup 2015 and European ST-T Database and the rationale has been clearly explained in Section \ref{disc}.

\begin{figure*}[!htbp]
\centering
\includegraphics[width = 0.80\textwidth]{figs/box.png}
\caption{Box and whisker plot for different SNR levels of test data of MIT-BIH Arrhythmia Database. On the x-axis we have SNR values starting from 36dB up to -36dB and on y-axis RMSE (ms) values are listed. The RMSEs are calculated by comparing the estimated IBIs with the true ones. Circular dots with low opacity are different RMSE values. The horizontal edges of colored boxes mark the interquartile range (IQR) (Q\textsubscript{1}-Q\textsubscript{3}) and contain 50\% data whereas the black horizontal lines inside each box mark the median value. The vertical whiskers hold other 50\% of the data and black horizontal lines mark the minimum and maximum. $\smblklozenge$ symbols are the outliers.}
\label{fig:box}
\end{figure*}

The accuracy of the proposed methodology can also be inferred from the box and whisker plot depicted in Figure~\ref{fig:box}. The IBI values have been estimated for all 13 noisy variants of the test ECGs in MIT-BIH Arrhythmia Database listed in Table \ref{tab1} and compared with that of the clean versions. Notice that, we hardly get any error for SNR values up to -12dB and the estimated IBI values are consistent and nearly as good as the true IBIs. However, starting from -18dB, we encounter some amount of error and above -30dB the noise becomes very intense as well as the RMSE. The median value at -30dB is 7.792 ms.

\begin{figure}[!htbp]
\centering
\includegraphics[scale = 0.2]{figs/IBI_comp.png}
\caption{Comparison of IBI plot with true IBIs and estimated IBIs obtained by \textit{Brosnan et al.}. For higher resolution, only first 3 minutes have been plotted.}
\label{fig:IBI_comp}
\end{figure}

\begin{table*}[t]

\begin{subtable}[h]{\textwidth}
    \centering
\resizebox{0.8\textwidth}{!}
\begin{tabular}{p{0.18\linewidth}p{0.1\linewidth}p{0.1\linewidth}p{0.1\linewidth}p{0.1\linewidth}p{0.1\linewidth}p{0.1\linewidth}}
\hline
 & 0dB & -6dB & -12dB & -18dB & -24dB & -30dB \\
\hline
\multirow{2}{10em}{Brosnan et al.\cite{Yuen28}} & 1.69 ms & 6.86 ms & 15.41 ms & 21.22 ms & 28.73 ms & 31.44 ms \\
& 2.38\% & 6.97\% & 9.06\% & 12.74\% & 14.34\% & 17.24\% \\ \hline
\multirow{2}{10em}{Sricharan et al.\cite{Vijayarangan26}} & 2.64 ms & 8.61 ms & 12.02 ms & 20.05 ms & 49.78 ms & 54.33 ms \\
& 1.54\% & 3.79\% & 9.87\% & 18.71\% & 20.97\% & 22.41\% \\ \hline
\multirow{2}{10em}{Aygun et al.\cite{Aygun18}} & 5.08 ms & 8.41 ms & 17.24 ms & 32.47 ms & 42.32 ms & 55.72 ms \\ & 3.54\% & 5.88\% & 8.77\% & 14.27\% & 17.31\% & 21.83\% \\ \hline
\multirow{2}{10em}{Proposed Method} & 0.08 ms & 0.17 ms & 0.68 ms & 1.09 ms & 4.72 ms & 8.39 ms \\
& 0.01\% & 0.08\% & 0.28\%  & 0.95\% & 4.23\% & 6.34\% \\ \hline

    \end{tabular}}
    \caption{MIT-BIH Arrhythmia Database}
    \label{tab3}
\end{subtable}

\hfill

\begin{subtable}[h]{\textwidth}
    \centering
\resizebox{0.8\textwidth}{!}
\begin{tabular}{p{0.18\linewidth}p{0.1\linewidth}p{0.1\linewidth}p{0.1\linewidth}p{0.1\linewidth}p{0.1\linewidth}p{0.1\linewidth}}
 &  &  &  &  &  &  \\
\hline
 & 0dB & -6dB & -12dB & -18dB & -24dB & -30dB \\
\hline
\multirow{2}{10em}{Brosnan et al.\cite{Yuen28}} & 1.02 ms & 4.27 ms & 12.57 ms & 23.29 ms & 25.13 ms & 30.53 ms \\
& 1.87\% & 1.61\% & 2.59\% & 9.61\% & 15.02\% & 18.72\% \\ \hline
\multirow{2}{10em}{Sricharan et al.\cite{Vijayarangan26}} & 2.21 ms & 7.07 ms & 13.74 ms & 19.69 ms & 27.92 ms & 37.81 ms \\
& 1.21\% & 3.17\% & 11.51\% & 17.94\% & 21.38\% & 27.33\% \\ \hline
\multirow{2}{10em}{Aygun et al.\cite{Aygun18}} & 5.32 ms & 7.95 ms & 15.33 ms & 33.41 ms & 37.97 ms & 47.86 ms \\ & 3.21\% & 5.12\% & 7.91\% & 12.02\% & 18.32\% & 19.27\% \\ \hline
\multirow{2}{10em}{Proposed Method} & 2.27 ms & 3.85 ms & 4.09 ms & 5.76 ms & 7.98 ms & 10.52 ms \\
& 0.85\% & 1.32\% & 1.77\% & 2.31\% & 3.87\% & 7.71\% \\ \hline

    \end{tabular}}
    \caption{European ST-T Database}
    \label{tab4}
\end{subtable}

\hfill
\linebreak[1]
\begin{subtable}[h]{\textwidth}
    \centering
\resizebox{0.8\textwidth}{!}
\begin{tabular}{p{0.18\linewidth}p{0.1\linewidth}p{0.1\linewidth}p{0.1\linewidth}p{0.1\linewidth}p{0.1\linewidth}p{0.1\linewidth}}
 &  &  &  &  &  &  \\
\hline
 & 0dB & -6dB & -12dB & -18dB & -24dB & -30dB \\
\hline
\multirow{2}{10em}{Brosnan et al.\cite{Yuen28}} & 2.35 ms & 8.41 ms & 15.14 ms & 18.73 ms & 21.61 ms & 32.64 ms \\
& 2.17\% & 4.47\% & 7.37\% & 9.22\% & 13.36\% & 16.04\% \\ \hline
\multirow{2}{10em}{Sricharan et al.\cite{Vijayarangan26}} & 4.97 ms & 8.82 ms & 15.26 ms & 25.41 ms & 32.63 ms & 38.59 ms \\
& 1.73\% & 3.59\% & 9.85\% & 15.81\% & 17.23\% & 21.75\% \\ \hline
\multirow{2}{10em}{Aygun et al.\cite{Aygun18}} & 6.21 ms & 9.91 ms & 17.13 ms & 31.04 ms & 35.31 ms & 41.23 ms \\ & 3.64\% & 5.89\% & 7.42\%  & 11.90\% & 18.47\% & 18.98\% \\ \hline
\multirow{2}{10em}{Proposed Method} & 5.14 ms & 7.11 ms & 9.27 ms & 11.91 ms & 17.84 ms & 26.33 ms \\
& 3.63\% & 4.22\% & 6.96\% & 8.05\% & 8.79\% & 11.23\% \\ \hline

    \end{tabular}}
    \caption{IEEE Signal Processing Cup 2015}
    \label{tab5}
\end{subtable}

\caption{Comparison of different methods on 3 dataset. For any method, the upper value refers to avg. RMSE and the lower value refers to avg. \% error. Performance have been measured on signals listed in Table~\ref{tab1} only.}
\label{tab_comp}
\end{table*}

Plots of true IBIs, IBIs estimated with our approach and that of \textit{Brosnan et al.} (both at -30dB) are pictured in Figure~\ref{fig:IBI_comp}. Some other methods have not been included with a view to keeping the plots clean to understand. IBIs estimated with our approach are not that much different from the true one following the trend in a cogent way.

Lastly, we preform comparative analysis of 3 state-of-ther-art and our proposed method on MIT-BIH Arrhythmia Database v1.0.0, European ST-T Database, and IEEE Signal Processing Cup 2015 in ~\crefrange{tab3}{tab5} (the details of these methods are described in Section 2). Once more, the performance metrics are average RMSE and average \% error and noise has been added to all the signals from each dataset. In these tables we have considered SNR levels starting from 0dB to -30dB since performance of all methods on higher quality signals are almost flawless. The weighted average of RMSE is 13.514 ms and \% error is 7.97\% for the proposed method whereas they are 31.467 ms (17.953 ms higher than us) and 17.384\% (9.414\% higher) for \textit{Brosnan et al.}, 45.439 ms (31.925 ms higher) and 23.72\% (15.75\% higher) for \textit{Sricharan et al.}, 49.739 ms (36.225 ms higher) and 20.349\% (12.379\% higher) for \textit{Aygun et al.}.

\section{Discussion}
\label{disc}

While analyzing the proposed method's performance, we have come across some of the observations, facts, and limitations which are worth mentioning-

The proposed algorithm can make the peaks more visible, detect them out of intense noisy signals and estimate the IBI values from it. However, it cannot redeem or rescue the typical shape of an ECG i.e. the PQRST wave. It can only recover the R-peaks.

Also, the IEEE SP Cup 2015 data is already encrypted with motion artifacts and we are adding additional noise to it. So, the unexpected performance on this dataset is quite reasonable.

It can be inferred from table \ref{tab_comp} that the performance of proposed algorithm on European ST-T Database and IEEE Signal Processing Cup 2015 is below than that of MIT-BIH Arrhythmia Database. Since the later is already sampled at 360Hz, we do not have to upsample it to add noise like the other two. So, upsampling of data degrades the performance on it.

\begin{figure}[h]
\centering
\includegraphics[scale = 0.225]{figs/dist.png}
\caption{A distorted ECG with missing peaks.}
\label{dist}
\end{figure}

Some of the signals are not clean enough, missing peaks, distorted and have some innate noise. They are beyond the scope of fixing with any sort of pre processing. Such a distorted signal segment is presented in Figure~\ref{dist}. Obtaining IBIs from these types of signals is nearly impossible even in their pure form. So we had to neglect such signals from all database.

\begin{figure}[h]
\centering
\includegraphics[scale = 0.225]{figs/inconsistent_peak_2.png}
\caption{An ECG with inconsistent peak size. There are too many such signals in each of the database. They have been omitted since they decrease the performance.}
\label{inc_p}
\end{figure}

Some of the signals have inconsistent peak size, like the one in Figure~\ref{inc_p}. This becomes a notable hindrance. The SNR levels are set according to the power of the ECG and power of noise. Now, in an ECG with inconsistent peak size, the consistent peaks play dominant role to set the power of ECG signal and the noise power is set accordingly. In such a case, the noise peaks adjacent to the inconsistent peaks become much higher and the tiramisu model fails to differentiate between the actual ECG peak and the noise. So, it suppresses the actual peak and retains the noisy peak leading to higher error in IBI estimation.

Same goes for signals with SNR beyond -30dB. The noise signal gets so intense that it cannot be suppressed anymore. The detected peaks are not perfect and the IBI values estimated are far from accurate. So, we claim, our method is good for estimation up to -30dB.

\section{Conclusion}
\label{conc}
A stacked tiramisu model based novel method for extracting IBI from an ECG that is submerged in motion artifacts is presented and analyzed in this paper. This method can be easily used for this purpose and can compete with traditional concurrent techniques. The model simply suppresses the noise and later, the most prominent feature of ECG i.e. position of R-peaks, has been leveraged to complete the task. Also, it doesn't require any pre-processing or post-processing method which saves some of the complexities. Results show the efficiency of the proposed model and the estimated IBIs are highly correlated to that of the true IBIs even at -30dB SNR value. The weighted average RMSE for IBI estimation is 13.514 ms. Although the size of the model is massive, with some lighter versions of deep learning tools, it can also be made suitable for embedded system applications. On the other hand, this approach is easily understandable, handy, highly accurate, and bears the proof of practical feasibility. Also, with some slight changes, this method can also be used for other physiological signals. Undoubtedly, it is a deep learning based approach implemented for the betterment of mankind and we hope it will have far reaching effect in the field of cardiology since it can work as an early indicator of many cardiac diseases.

\bibliographystyle{unsrt}


\section{Introduction} \label{sec:intro}

A relationship exists between galaxy properties and local environmental density \citep{DRESSLER}, such that early-type galaxies (ETGs) are preferably found in denser environments. Although the kinematic properties of galaxies do not correlate one-to-one with visual morphology, a tentative relationship between a proxy for the spin parameter, $\lambda_{R_e}$, and environmental density was presented by \cite{ATLAS3DVII}. This kinematic morphology-density relation (KMDR) suggests that the fraction  of slow rotating galaxies (low $\lambda_{R_e}$) increases towards denser environments. However, galaxy stellar mass also correlates with both environment and the intrinsic properties of galaxies. Thus, the question arises what the true physical driver of the KMDR is.

Further work supported the picture of a KMDR \citep{DEUGENIO2013, HOUGHTON2013,SCOTT2014,FOGARTY2014}. More recent results with larger galaxy samples find that the KMDR is driven mostly by stellar mass \citep{BROUGH2017,VEALE2017-2,GREENE2017}, but that the KMDR may also still exist at fixed stellar mass \citep{GRAHAM2019}. There is also evidence from simulations that points towards environmental dependence as a weaker secondary effect, but mass as the primary physical driver \citep[e.g.][]{LAGOS2017}. Some clear results have emerged, i.e. the trend with mass, but it is evident that any environmental dependence is likely to be a second order effect.

Multi-object integral field spectroscopy has revolutionised the number of galaxies with spatially resolved kinematic measurements. The Sydney-Australian-Astronomical-Observatory Multi-object Integral-Field Spectrograph (SAMI) galaxy survey has observed $\sim$3000 galaxies \citep{SAMII}, while the Mapping Nearby Galaxies at Apache Point Observatory (MaNGA) survey aims to observe $\sim$10,000 galaxies \citep{MANGAI}. Other ancillary surveys such as the Sloan Digital Sky Survey \citep[SDSS;][]{YORK2000} and the Galaxy and Mass Assembly Survey \citep[GAMA;][]{DRIVER2011} enable an accurate definition of environment, tracing the underlying large scale structure that exists in the Universe.

With the growing wealth of spatially resolved kinematic data, the statistical tool of correlation functions becomes more powerful. It allows us to connect large scale structure in galaxy clustering to internal galaxy properties, in our case galaxy spin, parameterised by \ifmmode{\lambda_{R_{\rm{e}}}}\else{$\lambda_{R_{\rm{e}}}$}\fi\. Correlation functions have already demonstrated a relation between environment and the star-forming and morphological properties of galaxies \citep[e.g.][]{Madgwick,HERMIT1996}. Marked correlation functions, where galaxies are marked by some physical parameter, are even more effective at detecting and quantifying weak correlations with environment \citep{SHETH2004,HARKER2006},  making this method ideal for detecting a possible relation between $\lambda_{R_e}$ and environment.

\update{Most papers look at the fraction of fast and slow rotators as a function of mass and environment \citep{ATLAS3DVII,BROUGH2017,VEALE2017,GREENE2017,GRAHAM2019}, whereas a broader analysis of the \ifmmode{\lambda_{R_{\rm{e}}}}\else{$\lambda_{R_{\rm{e}}}$}\fi\ distribution as a function of mass and environment shows that environment might have a small impact on \ifmmode{\lambda_{R_{\rm{e}}}}\else{$\lambda_{R_{\rm{e}}}$}\fi\ \citep{WANG2020}.}

In this paper we aim to investigate the correlation between $\lambda_{R_e}$ and environmental density. We present an analysis using marked cross correlation functions applied to SAMI \citep{SAMII} and GAMA \citep{DRIVER2011} data, as well as mock observations from the EAGLE Simulations \citep{EAGLEI}. We adopt a $\Lambda$CDM cosmology, with $H_0=70\text{ km} \text{ s}^{-1} \text{ Mpc}^{-1}$, $\Omega_m=0.3$, $\Omega_{\Lambda}=0.7$.

\section{Observations and Simulations}
\label{sec:data}
\subsection{Observations}

The SAMI instrument \citep{SAMII} is mounted on the Anglo-Australian Telescope and provides a 1 degree diameter field of view. SAMI employs 13 fused fibre bundles \citep[Hexabundles;][]{BLANDHAWTHORN2011,BRYANT2014} with a high (75\%) fill factor. Each bundle contains 61 fibres of 1.6\arcsec diameter resulting in each IFU having a diameter of 15 \arcsec. The IFUs, as well as 26 sky fibres, are fed to the AAOmega spectrograph \citep{SHARP2006}, using the 580V grating at 3570-5750A giving a resolution of R=1808 ($\sigma$=70.4 \ifmmode{\,\rm{km}\, \rm{s}^{-1}}\else{$\,$km$\,$s$^{-1}$}\fi), and the R1000 grating from 6300-7400A giving a resolution of R-4304 ($\sigma$=29.6 \ifmmode{\,\rm{km}\, \rm{s}^{-1}}\else{$\,$km$\,$s$^{-1}$}\fi) \citep{VANDESANDE2017}.

The SAMI Galaxy survey \citep{SAMII,BRYANT2015} selected galaxies from the GAMA \citep{DRIVER2011} survey, in addition to eight low-redshift clusters \citep{OWERS2017}. Reduced data cubes \citep{SHARP2015} and stellar kinematic maps are available with the SAMI Galaxy Survey data releases \citep{ALLEN2015,GREEN2018,SCOTT2018,CROOM2021}.

Our sample contains \update{1832 galaxies} with $\lambda_{R_e}$ measurements, derived from spatially resolved kinematic measurements as described in \cite{VANDESANDE2017}, and include an aperture correction \citep{VANDESANDE2017_MNRAS} and a seeing correction \citep{HARBORNE2020,VANDESANDE2020}. Furthermore, we define a volume-limited sample by selecting galaxies at: 1) $z < 0.06$ to avoid biases in the marked correlation function as the SAMI selection function results in different distributions of galaxy stellar masses in different volumes, and 2) $M/M_{\odot}>10^{10}$ to avoid low-completeness in the stellar kinematic sample. Alternatively, we could have treated each redshift section within the SAMI function selection individually, and taken a variance weighted mean of each resulting correlation function. While this achieves a stronger signal, the result becomes more difficult to interpret due to our ranking method for weights (Section \ref{sec:method}). The final sample contains \update{710 galaxies}, of all morphological types.\update{ This sample is unique amongst other similar surveys. For example, ATLAS$^{3\text{D}}$'s sample for kinematic analysis \citep{ATLAS3DVII} only contained 260 field ETGs, and MaNGA \citep[10,000 galaxies,][]{MANGAI} is only a narrow band in the mass redshift plane, with a complicated selection function that may make analyses that require volume limited samples more difficult to replicate.}

Galaxies from the GAMA survey \citep{DRIVER2011,LISKE2015} serve as a background galaxy distribution for our analysis. GAMA adopted an r-band magnitude limit of $r_{\text{pet}} < 19.8$ mag. For our analysis, we restrict the GAMA sample at $z<0.06$ ( $z_{\text{max}}$ of SAMI) and define a volume-limited sample in redshift and $r$-band apparent magnitude. This was done so that the background distribution was not biased by a higher density of fainter galaxies at lower redshifts in an apparent magnitude limited sample. The GAMA data used in this paper came from three equatorial regions centred at $9^h$, $12^h$ and $14.5^h$ in RA, each of $12\times4$ deg$^2$.

\subsection{Simulations}
\label{sec:sims}

We use galaxy mock-observations from the EAGLE hydrodynamical cosmological simulation suite \citep{EAGLEI} as presented by \cite{LAGOS2018}. A total of $7\times 10^8$ galaxies were extracted from the 100 Mpc$^3$ box, where each baryonic particle has an initial mass of $1.8\times 10^6$M$_{\odot}$, with a maximum gravitational softening length of 0.7 kpc. We adopt a stellar mass cut of M$_{\text{stars}} >5\times 10^9$ M$_{\odot}$, to ensure galaxies had angular momentum profiles that converged. $\lambda_{R_e}$ values were derived for mock-observations of this sample, leaving us with 5587 galaxies. A sample of 29737 galaxies not subjected to the stellar mass cut serves as the background galaxy distribution. These samples have an effective mass limit of $M/M_{\odot}>10^{8.5}$.

Additionally, as the mass distributions of SAMI and EAGLE galaxies are significantly different, we use a set of galaxies sampled from EAGLE in such a way to match the SAMI mass distribution. The initial distributions can be seen in Figure \ref{fig:mass_spin_dist}.
\begin{figure*}
 \includegraphics[width=\textwidth]{mass_spin_dist_final.pdf}
 \caption{The distribution of SAMI (red) and EAGLE (blue) galaxies in $\text{log}_{10}(M/M_{\odot})$-$\lambda_{R_e}$ space. The normalised histograms of both samples for both parameters is also included. In this paper we adopt a mass limit for SAMI of $\log_{10}(M/M_{\odot})> 10$.}
 \label{fig:mass_spin_dist}
\end{figure*}

\section{Marked Correlation Functions}
\label{sec:method}
Marked statistics are a powerful tool to determine whether correlations between galaxy parameters depend on environment or not. Galaxies are assigned a mark, corresponding to some physical parameter, and a marked correlation function \citep{SHETH2004,SHETH2005} is calculated. We begin by defining the 2-point real-space correlation function $\xi(r)$, given by \cite{PEEBLES}:
\begin{equation}\label{eq:corr_func_def}
    dP=\rho^2[1+\xi(r)]dV_1dV_2,
\end{equation}
where $dP$ is the probability for two galaxies to be located at a distance $r$ from each other, in volume elements $dV_1$ and $dV_2$. $\rho$ is the mean density of galaxies in the volume considered. The correlation function $\xi(r)$ measures the ``overdensity''  of galaxies. In our case, we assume that redshift gives the radial distance to a galaxy. This is known as a redshift-space ($s$) correlation function, $\xi(s)$.

In practice, estimators are used to calculate $\xi(s)$. We use an estimator introduced by \cite{PEEBLES}:
\begin{equation}
    1+\xi(s)=\frac{\langle SG(s)\rangle}{\langle SG_R(s)\rangle},
\end{equation}
where $\langle SG(s)\rangle$ and $\langle SG_R(s)\rangle$ are pair counts between a SAMI ($S$) and GAMA ($G$) galaxy, and a SAMI and random GAMA ($G_R$) galaxy respectively. These pair counts are calculated by first defining radial bins, equally spaced in $\log_{10}(s)$. A SAMI galaxy is selected, and the distance to each GAMA galaxy is calculated. A count is then added to each relevant radial bin. This is repeated for all SAMI galaxies. Random GAMA galaxies are a sample of GAMA-like galaxies, created in such a way to match GAMA's selection function.

We now assign SAMI galaxies marks ($m$). Galaxies are ranked in $\lambda_{R_e}$ and log stellar mass, with the mark taken as the rank, to ensure an equivalent dynamical range for both marks. Marked correlation functions can then be thought of as a ratio of galaxy marks to the mean mark, $\overline{m}$, as a function of galaxy separation \citep{SHETH2005}:
\begin{align}\label{eq:marked_def}
    M(s) & \equiv \frac{\sum m(\boldsymbol{x})m(\boldsymbol{y})\mathcal{I}(|\boldsymbol{x}-\boldsymbol{y}|-s)}{\overline{m_x}\overline{m_y}\sum\mathcal{I}(|\boldsymbol{x}-\boldsymbol{y}|-s)}\nonumber\\
    &=\frac{\sum m(\boldsymbol{x})\mathcal{I}(|\boldsymbol{x}-\boldsymbol{y}|-s)}{\overline{m}\rho^2[1+\xi(s)]},
\end{align}
where $m(\boldsymbol{x})$ is the mark of a SAMI galaxy. We have defined the mark of all GAMA galaxies [$m(\boldsymbol{y})$] to be $1$, $\mathcal{I}(x)=0$ unless $x=0$, and the sum is over all SAMI-GAMA galaxy pairs. As we have divided by $\overline{m}$, $M(s)=1$ for all $s$ if no correlation between marks and environment exists.

We can also consider Equation \ref{eq:marked_def} in an alternate way. By a simple re-arrangement, the denominator can be expressed as one plus the regular correlation function, defined in Equation \ref{eq:corr_func_def}, and the numerator as one plus a ``weighted'' correlation function. This weighted correlation function, defined as $W(s)$, can be calculated using the same estimator as $\xi(s)$, except the $i^{\text{th}}$ SAMI galaxy contributes a weight of $m_i/\overline{m}$ to the relevant radial bins:
\begin{align}
    M(s)&\equiv \frac{\sum (m(\boldsymbol{x})/\overline{m})\mathcal{I}(|\boldsymbol{x}-\boldsymbol{y}|-s)}{\sum\mathcal{I}(|\boldsymbol{x}-\boldsymbol{y}|-s)}\equiv \frac{1+W(s)}{1+\xi(s)}\nonumber\\
    &=\frac{\langle WG(s) \rangle}{\langle SG(s) \rangle},
\end{align}
where $\langle WG(s) \rangle$ are \textit{weighted} pair counts between weighted SAMI ($W$) and GAMA ($G$) galaxies, and $\langle SG(s) \rangle$ are \textit{unweighted} pair counts between SAMI ($S$) and GAMA galaxies. Marked correlation functions with random marks are also calculated, as a check that any signal seen in the real functions is legitimate. Once our correlation functions are calculated, we use the Python emcee package \citep{EMCEE} to fit a function of the form:
\begin{equation}\label{eq:power_law}
    M(s)=1+As^{-m}
\end{equation}
We use this functional form as we expect $M(s)\approx 1$ for all scales other than small scales, where it may deviate according to possible spin-environment and mass-environment relations. Uncertainties are calculated from the $16^{\text{th}}$ and $84^{\text{th}}$ percentiles.
\subsection{Uncertainty Calculation}
\label{sec:errors}
\update{We choose bootstrap re-sampling for our uncertainty estimate, as bootstrap re-sampling is robust and as shown by \cite{FISHER1994}, at worst overestimates error in correlation functions. Importantly, bootstrap uncertainties assume zero correlation between points, which is not strictly true in our case, as single galaxies contribute to multiple pair counts. However, on small scales the data points are largely independent, due to the small number of galaxies contributing to pairs at such small separation. As another test, we also evaluated another two different error estimations. We divided our sample into nine regions, and calculated a correlation function in each region, taking the standard deviation between the regions. We also calculated Poissonian errors for each bin. We found that the bootstrap errors were similar to Poissonian at small scales, and similar to the nine region standard deviation at large scales. Due to this, and our largely independent points at small scales, we use ordinary bootstrap re-sampling as a close approximation of errors.}

As bootstrap re-sampling generally overestimates uncertainties at small scales by a factor $\sqrt{3}$ \citep{CROOM1996}, the number of SAMI galaxies drawn per bootstrap sample was $3N$, where $N$ is the total number of galaxies, as also suggested by \cite{NORBERG2001}. Previous galactic correlation function works have used 10 re-samples \citep{HERMIT1996} or 20 re-samples \citep{Madgwick}. Due to our smaller data set in SAMI, we created 10,000 re-samples for each correlation function.

\section{Results}
\label{sec:results}
\subsection{SAMI Galaxies}
We present marked correlation functions for SAMI Galaxies, with ranked $\lambda_{R_e}$ and stellar mass marks, in Figure \ref{fig:corr_funcs_sami}.

Towards small scales $s$ we find that the marked correlation measurements (orange symbols) are significantly less than 1, i.e., there is a significant anti-correlation of ranked \ifmmode{\lambda_{R_{\rm{e}}}}\else{$\lambda_{R_{\rm{e}}}$}\fi\ with environment. The best-fit power-law (Equation \ref{eq:power_law}, Figure \ref{fig:corr_funcs_sami}a, red line), shows a significant turn below $M(s)=1$, beginning at around $s=1$ Mpc. \update{We find a best-fit value for $A=-0.038^{+0.012}_{-0.013}$, $3.2\sigma$ below zero}. This implies that galaxies with low $\lambda_{R_e}$ start being preferentially located in dense environments at scales of $s\approx 1$ Mpc. As done by \cite{HARKER2006}, we also take one large radial bin out to 1 Mpc to find significance at small scales ($M(<1Mpc)$). \update{We find $M(<1Mpc)=0.925\pm 0.034$, $2.2\sigma$ below 1}. Randomised marks are consistent with $M(s)=1$ at all scales.

In Figure \ref{fig:corr_funcs_sami}b we also detect a significant correlation of ranked stellar mass with environment. The best-fit power-law (red line) shows a significant up turn above $M(s)=1$ towards small scales in $s$, beginning at around $s=1$ Mpc. \update{We find a best-fit value of $A=0.036^{+0.015}_{-0.015}$, $2.4\sigma$ above zero}. This result implies that galaxies with high stellar mass start being preferentially located in dense environments at scales of $\approx 1$ Mpc. \update{We find $M(<1Mpc)=1.058\pm 0.042$, $1.4\sigma$ above zero}. This result is consistent with previous work, as galaxies are well known to cluster according to mass (or luminosity) \citep[e.g.][]{NORBERG2002}. Randomised marks are consistent with $M(s)=1$ at all scales.

\begin{figure*}
 \includegraphics[width=\textwidth]{sami_spin_mass_final_new.pdf}
 \caption{Marked correlation functions for ranked \ifmmode{\lambda_{R_{\rm{e}}}}\else{$\lambda_{R_{\rm{e}}}$}\fi\ (left panel) and ranked stellar mass marks (right panel), with SAMI galaxies (blue data). Randomised marks are plotted as a test of our method (orange data). We plot the best-fit power law (dark red line). One hundred samples from the MCMC chain (light red lines) are plotted, as a representation of the uncertainty in the fits. Contours for the 68$^{th}$, 90$^{th}$ and 95$^{th}$ percentiles are plotted for the real data in $A$ and $m$ space, as explored by the emcee algorithm, in the inset. The red overlay lines in these parameter spaces represent the values of $A$ and $m$ selected at the $50^{\text{th}}$ percentile. Measures of the significance can be found in Table \ref{tab:significance}. \ifmmode{\lambda_{R_{\rm{e}}}}\else{$\lambda_{R_{\rm{e}}}$}\fi\ is negatively correlated with environment, while stellar mass is positively correlated with environment in SAMI galaxies.}
  \label{fig:corr_funcs_sami}
\end{figure*}
\subsubsection{Further Analysis of SAMI Results}

We have found an anti-correlation of \ifmmode{\lambda_{R_{\rm{e}}}}\else{$\lambda_{R_{\rm{e}}}$}\fi\ with environment, and a positive correlation of stellar mass with environment. We can compare these results with different tests. First, we reverse the rank order for \ifmmode{\lambda_{R_{\rm{e}}}}\else{$\lambda_{R_{\rm{e}}}$}\fi\ and compare $M(<1Mpc)$ for these reverse \ifmmode{\lambda_{R_{\rm{e}}}}\else{$\lambda_{R_{\rm{e}}}$}\fi\ ranks and stellar mass ranks. \update{We find $M(<1Mpc)=1.075\pm 0.043$ for reverse \ifmmode{\lambda_{R_{\rm{e}}}}\else{$\lambda_{R_{\rm{e}}}$}\fi\, and $M(<1Mpc)=1.058\pm0.041$ for stellar mass. This gives a difference of $0.017\pm0.059$}. Using this metric, mass and reverse \ifmmode{\lambda_{R_{\rm{e}}}}\else{$\lambda_{R_{\rm{e}}}$}\fi\ rank marks are consistent. However, this simplistic approach of assuming a one-to-one relationship between reversed \ifmmode{\lambda_{R_{\rm{e}}}}\else{$\lambda_{R_{\rm{e}}}$}\fi\ ranks and stellar mass ranks is not entirely supported, as there is a lot of scatter in this relationship \citep{CROOM2021}. In Figure \ref{fig:randoms_and_noslowrot}, we present two methods of testing for the physical driver of our SAMI $\lambda_{R_e}$-environment anti-correlation. Bins are defined in stellar mass of 0.1 dex in width. Galaxies are then assigned a \ifmmode{\lambda_{R_{\rm{e}}}}\else{$\lambda_{R_{\rm{e}}}$}\fi\ value from another random galaxy in their mass bin. This failed to reproduce a negative correlation with environment. The lack of correlation in this case implies that the signal we see between \ifmmode{\lambda_{R_{\rm{e}}}}\else{$\lambda_{R_{\rm{e}}}$}\fi\ and environment cannot be caused purely by a mass-environment relation. Removing slow rotators ($\lambda_{R_e}<0.2$) from the sample did reproduce a negative correlation with environment. This existence of a correlation implies that the observed correlation between \ifmmode{\lambda_{R_{\rm{e}}}}\else{$\lambda_{R_{\rm{e}}}$}\fi\ and environment cannot be attributed solely to slow rotators.

\begin{figure*}
 \includegraphics[width=\textwidth]{assigned_spin_nosr_final_new.pdf}
 \caption{Marked correlation functions for different tests exploring whether our SAMI \ifmmode{\lambda_{R_{\rm{e}}}}\else{$\lambda_{R_{\rm{e}}}$}\fi\ anti-correlation with environment can be explained by mass trends. Selecting a galaxy's \ifmmode{\lambda_{R_{\rm{e}}}}\else{$\lambda_{R_{\rm{e}}}$}\fi\ value at random from 0.1 dex mass bins failed to reproduce a negative correlation (left panel). Removing slow rotators from the SAMI sample ($\ifmmode{\lambda_{R_{\rm{e}}}}\else{$\lambda_{R_{\rm{e}}}$}\fi<0.2$) still reproduces a negative correlation (right panel). Randomised marks are also plotted as a test of our method. Symbols and lines are similar as in Figure \ref{fig:corr_funcs_sami}, but now presented for mass-defined \ifmmode{\lambda_{R_{\rm{e}}}}\else{$\lambda_{R_{\rm{e}}}$}\fi\ and for a sample with slow rotators removed. The anti-correlation of \ifmmode{\lambda_{R_{\rm{e}}}}\else{$\lambda_{R_{\rm{e}}}$}\fi\ with environment is not driven purely by stellar mass, or by slow rotators.}
  \label{fig:randoms_and_noslowrot}
\end{figure*}
\subsection{EAGLE Galaxies}

We show the distribution in \ifmmode{\lambda_{R_{\rm{e}}}}\else{$\lambda_{R_{\rm{e}}}$}\fi\ and stellar mass of SAMI and EAGLE galaxies in Figure \ref{fig:mass_spin_dist}. As the mass distributions are significantly different, we create a set of observations sampled from EAGLE in such a way to match the SAMI mass distribution. Fifty bins spaced in $\log(M/M_{\odot})$ were defined across the mass range for SAMI, and EAGLE galaxies in each bin were drawn randomly until there were three times as many EAGLE galaxies in each bin as SAMI.

We present marked correlation functions for these re-sampled EAGLE galaxies, with ranked $\lambda_{R_e}$ and stellar mass marks, in Figure \ref{fig:corr_funcs_sim_resampled}. We see a significant anti-correlation of ranked $\lambda_{R_e}$ with environment in our marked correlation function for EAGLE (Figure \ref{fig:corr_funcs_sim_resampled}a). The best-fit power-law shows a significant downward trend towards small scales in $s$. $m$ is lower for EAGLE than for SAMI, so the correlation between \ifmmode{\lambda_{R_{\rm{e}}}}\else{$\lambda_{R_{\rm{e}}}$}\fi\ and environment extends to larger scales in EAGLE than in SAMI. We find a best-fit value of $A=-0.020^{+0.003}_{-0.003}$, $6.6\sigma$ below zero. This implies that in the EAGLE mock-observations, low-$\lambda_{R_e}$ galaxies are preferentially located in dense environments, out to scales of the order $\sim$1 Mpc. This is consistent with our SAMI results. We find $M(<1Mpc)=0.966\pm0.008$, $4.2\sigma$ below 1. Randomised marks are consistent with $M(s)=1$ at all scales.

We see a significant correlation of ranked stellar mass with environment in our marked correlation function for EAGLE. The best-fit power-law shows a significant upwards trend towards small scales in $s$. We find a best fit value of $A=0.030^{+0.003}_{-0.003}$, $10\sigma$ above zero. This implies that in the EAGLE mock-observations, high stellar mass galaxies are preferentially located in dense environments, out to scales of the order $\sim10$ Mpc. We find $M(<1Mpc)=1.066\pm0.015$, $4.4\sigma$ above 1. Randomised marks are consistent with $M(s)=1$ at all scales, as expected.

\begin{figure*}
 \includegraphics[width=\textwidth]{sim_spin_mass_final.pdf}
 \caption{Marked correlation functions for ranked \ifmmode{\lambda_{R_{\rm{e}}}}\else{$\lambda_{R_{\rm{e}}}$}\fi\ (left panel) and ranked stellar mass marks (right panel). Symbols and lines are similar as in Figure \ref{fig:corr_funcs_sami}, but now presented for the EAGLE data. \ifmmode{\lambda_{R_{\rm{e}}}}\else{$\lambda_{R_{\rm{e}}}$}\fi\ is negatively correlated with environment, while stellar mass is positively correlated with environment in EAGLE galaxies.}
  \label{fig:corr_funcs_sim_resampled}
\end{figure*}

\begin{deluxetable*}{cclDl}\label{tab:significance}
\tablenum{1}
\tablecaption{Measurements of the significance of the signals in the correlation functions for SAMI and EAGLE data. We present the values of $A$ and $m$ for each data type and mark type, and the value of a single bin below 1 Mpc for each data type and mark type. The significance for the single bin is taken as the distance away from $M(s)=1$, the expected value given no correlation.}
\tablewidth{0pt}
\tablehead{
\colhead{Data Type} & \colhead{Mark Type} &  \colhead{$A$} &
\multicolumn2c{$m$} & \colhead{Single Bin Value}
}
\decimalcolnumbers
\startdata
SAMI Spin Data & Real & $-0.038^{+0.012}_{-0.013}$ & $1.060^{+0.196}_{-0.179}$ & $0.925\pm0.035$ \\[2pt]
  & Random & $-0.008^{+0.005}_{-0.008}$ & $1.850^{+0.412}_{-0.317}$ & $0.981\pm0.038$ \\[2pt]
SAMI Mass Data & Real & $0.036^{+0.015}_{-0.015}$ & $0.713^{+0.320}_{-0.255}$ & $1.058\pm0.041$ \\[2pt]
   & Random & $-0.003^{+0.009}_{-0.011}$ & $0.795^{+0.706}_{-0.512}$ & $1.019\pm0.039$ \\[2pt]
  Re-sampled EAGLE Spin Data & Real & $-0.020^{+0.003}_{-0.003}$ & $0.529^{+0.092}_{-0.094}$ & $0.966\pm 0.008$ \\[2pt]
   & Random & $-0.001^{+0.002}_{-0.003}$ & $0.530^{+0.480}_{-0.332}$ & $1.000\pm 0.008$ \\[2pt]
  Re-sampled EAGLE Mass Data & Real & $0.030^{+0.003}_{-0.003}$ & $0.488^{+0.045}_{-0.046}$ & $1.066\pm 0.015$ \\[2pt]
   & Random & $-0.003^{+0.003}_{-0.003}$ & $0.379^{+0.414}_{-0.225}$ & $1.005\pm 0.015$ \\[2pt]
  Mass Assigned Spin Data, SAMI & Real & $-0.002^{+0.009}_{-0.008}$ & $1.035^{+0.807}_{-0.641}$ & $0.984\pm 0.038$\\[2pt]
   & Random & $-0.014^{+0.011}_{-0.013}$ & $0.792^{+0.536}_{-0.412}$ & $0.989\pm 0.038$\\[2pt]
  No Slow Rotators Spin Data, SAMI & Real & $-0.035^{+0.013}_{-0.014}$ & $1.308^{+0.218}_{-0.183}$ & $0.937\pm 0.037$\\[2pt]
   & Random & $0.005^{+0.014}_{-0.009}$ & $0.812^{+0.662}_{-0.460}$ & $1.008\pm 0.041$
\enddata
\end{deluxetable*}

\section{Conclusion}
We investigate the existence of a kinematic morphology-density relation (KMDR) using the technique of marked correlation functions. Ranking galaxies by their spin (parametrised by $\lambda_{R_e}$) and stellar mass allows us to investigate whether the relation is driven purely by stellar mass or if some residual spin correlation exists. This is the first time that marked correlation functions have been applied to galactic spin. We apply our analysis to both the SAMI observations and EAGLE simulation mock-observations, allowing us to compare SAMI to simulated results. The GAMA dataset is used as the ``background'' in the SAMI analysis for measuring large scale structure, due to GAMA's highly complete spectroscopic data. We estimate uncertainties in all our correlation functions through bootstraps with 10,000 re-samples, and estimate significance by fitting power laws of the form $M(s)=1+As^{-m}$ and looking at single bins below 1 Mpc.

We find the following results:
\begin{itemize}
    \item Mass is positively correlated with environment. Fitting our ranked mass mark correlation function by a power law gave parameters $A=0.036^{+0.015}_{-0.015}$, and $m=0.713^{+0.320}_{-0.255}$. These results are consistent with previous work \citep[e.g.][]{NORBERG2002}.
    \item Spin (parametrised by $\lambda_{R_e}$) is negatively correlated with environment. Fitting our ranked $\lambda_{R_e}$ mark correlation function by a power law gave parameters $A=-0.038^{+0.012}_{-0.013}$, and $m=1.060^{+0.196}_{-0.179}$. This is consistent with previous work that found evidence of a KMDR \citep{ATLAS3DVII,DEUGENIO2013,HOUGHTON2013,SCOTT2014,FOGARTY2014}, and the presence of an environmental effect on spin \citep{CHOI2018}.
    \item Spin anti-correlation with environment is not driven solely by stellar mass. Defining bins in stellar mass of 0.1 dex, galaxies are assigned random \ifmmode{\lambda_{R_{\rm{e}}}}\else{$\lambda_{R_{\rm{e}}}$}\fi\ values corresponding to their mass bin. This failed to reproduce a negative correlation with environment.
    \item The anti-correlation between spin and environment is not driven purely by slow rotators. When slow rotators are removed from the SAMI sample ($\lambda_{R_e}<0.2$), the anti-correlation is found to still be present.
    \item Using mock-observations from the EAGLE simulations we find qualitatively consistent results as compared to observed data. Mass is correlated with environment, and spin is anti-correlated with environment with EAGLE galaxies as well as SAMI galaxies. A larger sample size in EAGLE allows us to see to larger scales than SAMI with high signal to noise, but the overall qualitative result remains the same. Although the relationship extends to larger scales for EAGLE, we know there are offsets between SAMI and EAGLE spin values \citep{VANDESANDE2019}. Investigation of whether these offsets can fully explain the larger scale result is beyond the scope of this paper.
\end{itemize}

Previous work has suggested that the $\lambda_{R_e}$-environment anti-correlation is simply a manifestation of the mass-environment correlation, caused by dynamical friction \citep{BROUGH2017,HOUGHTON2013,CAPPELLARI2016,VEALE2017-2,GREENE2017}. Our results suggest that the \ifmmode{\lambda_{R_{\rm{e}}}}\else{$\lambda_{R_{\rm{e}}}$}\fi\-environment anti-correlation is not as simple as a re-statement of the mass-environment correlation. Physical process such as environmental quenching, mergers and interactions may drive galaxies towards lower \ifmmode{\lambda_{R_{\rm{e}}}}\else{$\lambda_{R_{\rm{e}}}$}\fi\ values in denser environments (see also van de Sande et al., in Prep). The connection between mean stellar population age and environment \citep[e.g.][]{SCOTT2017} combined with the relation between mean stellar age and the dynamical thickness or intrinsic shape of galaxies \citep{vandesande2018} indicates that the star-formation history of a galaxy, its environment, and dynamical properties are closely connected.

Future IFS surveys such as HECTOR \citep{BRYANT2016} will improve the statistics that we find here, in particular at small separation.

\section*{ACKNOWLEDGEMENTS}
The SAMI Galaxy Survey is based on observations made at the Anglo-Australian Telescope. SAMI was developed jointly by the University of Sydney and the Australian Astronomical Observatory (AAO). The SAMI input catalogue is based on data taken from the Sloan Digital Sky Survey, the GAMA Survey and the VST ATLAS Survey. The SAMI Galaxy Survey is supported by the Australian Research Council (ARC) Centre of Excellence ASTRO 3D (CE170100013) and CAASTRO (CE110001020), and other participating institutions. The SAMI Galaxy Survey website is http://sami-survey.org/

JvdS, JJB, JBH, SB, MSO acknowledge support of the ARC (DE200100461, FT180100231, FL140100278, FT140101166, FT140100255). The SAMI instrument was funded by the AAO and JBH through FF0776384, LE130100198. FDE acknowledges funding through the H2020 ERC Consolidator Grant 683184.