Fundamental component enhancement via adaptive nonlinear activation functions

We start by describing the problem in the simplest possible setting: suppose we are given a $1-$periodic signal

and define the support of $S$ to be $S=\left\{k\in\mathbb{N}:a_{k}^{2}+b_{k}^{2}>0\right\}$. Our goal is to recover $\gcd(S)$, the greatest common divisor of all the elements in $S$. If we have recorded two or more periods of the signal, then the $\gcd$ will be bigger than 1. In that case, we call $\gcd(S)$ the fundamental frequency. Mathematically speaking, there is no difficulty: compute the Fourier transform of $f$ and check.

There are several reasons why this seemingly trivial problem is difficult in practice.

1. (1)

   In practice, all Fourier coefficients of $f$ will be distinct from 0; some form of thresholding to decide significance is required.

2. (2)

   The actual fundamental frequency may not be constant, it may change from time to time, or it may not exist at all [1, 4, 10, 11, 13].

3. (3)

   We only have $f$ sampled at finitely many points; moreover, in practice one should always expect some form of noise.

We describe an explicit problem coming from medical signal processing which was partially the motivation for our paper. Our example can be seen in Figure 2(a) and comes from a phonocardiogram (PCG) signal [2]: this is a sound signal created by the vibrations created by the closure of the heart valves. One fundamental component of the signal, one heartbeat, can be seen within each of the two red boxes: it is comprised of two ingredients, S1 and S2 (indicated by blue errors): S1 is due to the atrioventricular valves closing at the beginning of systole and S2 is the consequence of the aortic and pulmonary valves closing at the end of systole.

Due to the heart rate variability, the periods between two consecutive cycles are not fixed: an important but challenging problem is to estimate the heart rate at a local point in time, this amounts to estimating the time-varying period, or the instantaneous frequency. Naturally, one would be inclined to apply time-frequency (TF) analysis [3]. This signal oscillates roughly twice per second, so the frequency is around 2 Hz, and we would expect to see a dominant curve around 2 Hz that represents the fundamental component of the signal. However, we cannot see anything concrete around 2 Hz from the standard spectrogram shown in Figure 2(b) and, in this sense, the fundamental component of the signal does not exist. In short, we need a different solution if we want to estimate the time-varying heart rate.

There is a surprisingly simple and widely applied solution; that is, take the ‘rectification’ of the signal before running any sort of time-frequency analysis. Mathematically speaking, take the signal $\left|f(t)\right|$ where $f(t)$ is the signal of interest. We refer to Figure 3 where the signal from Figure 2 has been analyzed after rectification: we analyze $\left|f(t)\right|$ instead of $f(t)$.

Here, rectification leads to a dramatically different and much more informative spectrogram in which a clearly defined fundamental frequency exists (indicated by the red arrows). This simple solution has been widely applied in practice if we want to estimate the fundamental frequency, for example, in the problem of extracting fetal electrocardiogram (ECG) from the trans-abdominal maternal ECG [9] or extracting f-wave from the ECG signal of a patient with atrial fibrillation [6], which is a far from complete list. To the best of our knowledge, there is no theoretical argument explaining why the rectification trick works to enhance the fundamental component. We believe this to be of substantial theoretical interest.

Open problem. Why is standard time frequency analysis applied to $\left|f(t)\right|$ better able to recover the fundamental underlying frequency? Why does the ‘rectification’ trick work?

Abstracting the trick, we may interpret the rectification trick as the application of a nonlinear activation function $g(t)=|t|$: instead of analyzing $f(t)$, one analyzes $g(f(t))$. Naturally, there are many other choices of the activation function [7], like the commonly used Rectified Linear Unit (ReLU) function $g(t)=\mbox{ReLU}(f(t))$ and other widely applied activation functions from the theory of neural networks that are at our disposal. In practice, ReLU seems to work and work roughly as well as the absolute value. We do not know of any theoretical support for any of these functions and can now formulate the main question that motivates our paper.

Open problem. Which nonlinear activation function $g:\mathbb{R}\rightarrow\mathbb{R}$ leads to the ‘best’ recovery of the fundamental frequency?

One of the main points of our paper is to discuss the problem from the point of view of classical Fourier analysis and to propose, based on that, a somewhat unorthodox choice, provide theoretical support and discuss its behavior in practice.

We define, for $\varepsilon>0$ small, an adaptive activation function. Introducing, for $0<\varepsilon<1$, the function $h_{\varepsilon}:[-1,1]\rightarrow\mathbb{R}_{\geq 0}$

we propose to run standard time frequency on the normalized signal

and expect it to enhance the fundamental component of $f$. This activation function is adaptive to the input signal since it depends on the $L^{\infty}$ norm of the input signal. Note that $h_{\varepsilon}$ is also $2\pi-$periodic. So by construction, it inherits the periodicity of $f$. The question is now: can we determine the strength of the fundamental component of $h_{\varepsilon}$? To provide a theoretical analysis this question, we put some assumptions.

1. (1)

   We assume that the function $f:[0,2\pi]\rightarrow\mathbb{C}$ is of the type

   $$0\neq f(t)=\sum_{k=1}^{n}a_{k}e^{im_{k}t}\quad a_{k}\in\mathbb{C},m_{k}\in\mathbb{N},~{}0<m_{1}<m_{2}<\dots<m_{n}.$$

   In particular, this function will have fundamental frequency $\gcd(m_{1},\dots,m_{n})$.

2. (2)

   We assume that the function $g(t)=\left|f(t)\right|$ assumes its global maximum in finitely many points $\left\{t_{1},\dots,t_{m}\right\}\subset[0,2\pi]$ and

3. (3)

   that these maxima are non-degenerate: $g^{\prime\prime}(t_{i})<0$ for all $1\leq i\leq m$.

We can now state our result which shows that, as $\varepsilon\rightarrow 0$, the method will indeed recover the fundamental frequency (which we expect to be 1 in the generic setting).

Remarks.

1. (1)

   We note that for ‘generic’ signals (say, random signals picked from a suitable probability distribution), we expect that there is a unique maximum $t_{1}$: in that case, the leading order expansion is always guaranteed to be of order

   $$\widehat{h_{\varepsilon}}(1)=\frac{1}{\sqrt{\varepsilon}}\frac{\pi\cdot e^{it_{1}}}{\sqrt{-\frac{g^{\prime\prime}(t_{1})}{2\|g\|_{L^{\infty}}}}}+\mathcal{O}\left(\frac{1}{\varepsilon^{1/4}}\right)$$

   and we are guaranteed that the fundamental frequency is recovered.

2. (2)

   If there is more than one maximum, then we only fail at recovering the fundamental frequency if

   $$\sum_{j=1}^{m}\frac{e^{it_{j}}}{\sqrt{-g^{\prime\prime}(t_{j})}}=0.$$

   This is something that we do not expect to happen for generic signal (it corresponds to a precise algebraic identity).

3. (3)

   If the frequencies share a nontrivial greatest common divisor

   $$G=\gcd(m_{1},\dots,m_{n})>1,$$

   then the maxima of $g(t)$ will be $G-$periodic on the unit circle and the sum will vanish. The Theorem can then be applied to the function $g(t)=f(t/G)$ to conclude that the fundamental frequency $G$ will be recovered.

The purpose of this section is to give some perspective on the problem from the point of view of classical Fourier analysis. We note that the intrinsically nonlinear nature of the problem does pose an interesting challenge, nonetheless, there are some natural considerations partially inspired by a classical paper of Rudin [8] that seem like they might be relevant. We recall that our goal is, given a function $f$ of the type

to come with a nonlinear function $g:\mathbb{R}\rightarrow\mathbb{R}$ such that $g(f(t))$ has a clear frequency contribution at $\gcd\left(m_{1},\dots,m_{n}\right)$. The main difficulty in practice is, of course, that we do not actually know the precise formula of $f$ when only $f$ is given, and have no real idea what the actual frequencies $m_{1},\dots,m_{n}$ might be. We simplify things by assuming that the nonlinearity $g$ can be expanded into a Taylor series

Then each individual power applied to $f$ can be properly analyzed: note that

This shows that large powers naturally lead to a function whose frequencies is comprised of sums of the individual frequencies. If we replace powers by powers of the absolute value, then at least for even powers the identity $|x|^{2k}=x^{k}\overline{x}^{k}$ leads to a way of writing

as

There is a nontrivial contribution coming from the coefficients $(a_{k})_{k\in\mathbb{Z}}$, however, if we only consider the support of the frequencies, then a rather clear picture emerges: denoting the set of frequencies that make up $f$ by

we expect that

However, we also expect, for $k$ sufficiently large, that the numbers that can be appropriately written as sums and differences of elements in $M_{n}$ have a nice limiting behavior and (for a suitabye lenient definition of $\lim$)

These heuristic considerations suggest the following rough guidelines:

1. (1)

   The Fourier transform of $g(f(t))$ depends on the Taylor expansion of $g$.

2. (2)

   If the function $g$ is smooth, then the Taylor expansion will have rapidly decaying coefficients and sum sets of the type

   $$\underbrace{(M+M+\dots+M)}_{k~{}\mbox{\tiny times}}\quad\mbox{only contribute for small}~{}k.$$

   A choice like $g(x)=\cos{(x)}$ will not lead to a good detection function.

3. (3)

   This suggests functions $g$ for which the Taylor expansion decays slowly: functions that are discontinuous or have discontinuous first derivative.

4. (4)

   Moreover, introducing the absolute value $g(|x|)$ leads to sum-difference sets which might generally lead to better results.

These principles naturally suggest functions like $g(x)=|x|$ or $g(x)=\mbox{ReLu}(x)$. However, we note that both of these functions are actually continuous which leads to a decay of coefficients in the Taylor expansion. Indeed, considering all these principles, we should focus our attention on functions of the form

for which the $(b_{j})_{j=0}^{\infty}$ decay slowly. A particularly natural candidate is

and this is how we arrive at our proposed construction. As we will show in Section 3, this construction is not only well motivated by the aforementioned arguments but also naturally incorporates a second approach which is related to asymptotic analysis and encapsulated by Theorem 1.

We generate synthetic data to investigate how different activation functions enhance the fundamental component. Fix the sampling rate to be $512$ Hz. First, randomly select an integer $K$ between $5$ and $100$ according to a uniform distribution. Then, randomly select $K$ integers between $2$ and $250$ following the probability density function $Ce^{-(l/100)^{2}}$, where $C$ is the normalization constant, so that the chosen integers have 1 as the greatest common divisor. Denote the selected integers as $\{j_{k}\}_{k=1}^{K}$. Then, set the function $f(t)$ as

where $a_{k}$, $k=1,\ldots,K$, are identically and independently (i.i.d.) sampled from $(0,1]$, and $\phi_{k}$, $k=1,\ldots,K$, are i.i.d. sampled from $(0,2\pi]$. Note that by construction, $f(t)$ is a toy example as that shown in (1) with $\hat{f}(1)=0$ but $1$ as the fundamental frequency. For a given activation function $g(t)$, we evaluate the resulting fundamental component enhancement by evaluating the fundamental component energy ratio, defined as

Then, repeat the above procedure for $10^{5}$ times. The results of different activation functions are shown in Figure 5 as histograms of fundamental component energy ratios. Quantitatively, the median and median absolute deviation of the fundamental component energy ratio over $10^{5}$ times are 0.28% and 0.46% for the rectification, 0.07% and 0.12% for the RELU, 0.29% and 0.46% for $h_{0.2}$, 0.31% and 0.42% for $h_{0.1}$, and 0.33% and 0.33% for $h_{0.05}$.

We generate a semi-real dataset and compare the proposed activation function with $\epsilon=0.2$, $\epsilon=0.1$ and $\epsilon=0.05$ with the usually applied rectification and RELU activation function. We consider two databases. The first one is the noninvasive trans-abdominal ECG database from the PhysioNet Computing in Cardiology Challenge 2013 (CinC2013)¹¹1https://physionet.org/content/challenge-2013/1.0.0/. There are in total 75 recording, each recording contains 4 channels, and each record lasts for 1 minute. The second one is the Taiwan Integrated Database for Intelligent Sleep (TIDIS)²²2https://tidis.org. There are in total 20 whole night polysomnogram recordings, each recording contains 2 channels, and each record lasts for about 6 hours. For each recording in CinC2013, we detect the maternal R peaks by applying the standard R peak detection algorithm to the mean of 4 channels; for each recording in TIDIS, we apply the same standard R peak detection algorithm to the mean of two channels. For the $k$-th channel in the $i$-th recording, we randomly generate an 1-periodic function, denoted as $s_{ik}^{(0)}(t)$, by averaging $10$ randomly selected cardiac cycles. Then we intentionally remove the fundamental component via the Fourier transform, and denote the resulting 1-periodic function as $s_{ik}(t)$ that is sampled at 512 Hz. Next, we apply the $j$-th activation function to $s_{ik}(t)$, and denote the resulting signal as $g^{(j)}_{ik}(t)$. To evaluate how much the fundamental is activated, we record the energy ratio of the fundamental component defined in (4), denoted as $r^{(j)}_{ik}$. For the $k$-th channel in the $i$-th recording, the above procedure is repeated $100$ times. The histograms of all collected energy ratios for different activation functions are shown in Figure 6. We see that the proposed activation function gives a delta-like signal with the peak centered at the maximal value. Hence, the power spectrum is flatter compared with both rectification and RELU. This explains why the median of the energy ratio is smaller if we apply the proposed activation function since the proposed activation function tends to flatten the spectrum, particularly when $\epsilon$ is small. We shall mention that the main difference between these two databases is that the cardiac cycle is recorded from the abdomen in the CinC2013 database, and from the chest in the TIDIS database. Thus the cardiac cycles are of different morphology, and hence different energy ratio of the fundamental component. Note that after the averaging, the fetal ECG impact is reduced, and we obtain a reasonably clean maternal ECG cycles. See Figure 7 for more details.

To have a complete picture of the fundamental enhancement, we repeat the same data preparation procedure above, except removing the fundamental component; that is, we simulate the practical situation that the fundamental component may or may not exist, and its strength may or may not be weak. In addition to record the energy ratio, also denoted as $r^{(j)}_{ik}$ by a bit abuse of notation, after applying the activation function, we also record the energy ratio before, denoted as $r^{(0)}_{ik}$. In Figure LABEL:Figure6simu, we plot the 2-dim histogram of the distribution of $\{(\log(r^{(0)}_{ik}),\log(r^{(j)}_{ik})/\log(r^{(0)}_{ik}))\}$. By definition, when $\log(r^{(j)}_{ik})/\log(r^{(0)}_{ik}))>1$, we obtain an enhancement of the energy ratio of the fundamental component. The portion of cardiac cycles that the fundamental component strength is enhanced is recorded in the title. It is clear that all activation functions enhance the energy ratio of the fundamental component when it is weak.

Note that at the first glance, our proposed activation function seems to provide a worse result due to the smaller portion of enhanced cardiac cycles. However, we should emphasize that due to the “flattening” nature of our activation function, the energy ratio of the fundamental component is smaller.

Next, we consider the PCG signal as a real example. In Figure 8, we show the spectrogram of the PCG signal shown in shown in Figure LABEL:Figure3PCG(a) composed with RELU (left) and $h_{0.1}$. We could see that both RELU and $h_{0.1}$ provide the fundamental component information; that is, there is an identifiable curve around 2Hz. To have a quantification of this finding, we analyze the maternal PCG signals in the Shiraz University Fetal Heart Sounds Database³³3https://physionet.org/content/sufhsdb/1.0.1/. There are in total 92 recordings. See Figure 9(a) for a PCG signal different from that in Figure 2(a). It is clear that the spectrogram provides limited information about the heart rate. See Figure 9(b) and recall Figure 2(b) for an illustration. To quantify the fundamental component enhancement performance, we carry out the following steps. First, run the de-shape algorithm [5] to determine the instantaneous frequency (IF) of the PCG. See Figure 9(c-d) for the de-shape spectrogram and the detected IF. Denote the estimated IF as $\phi^{\prime}(t)$. The detected heart rate is confirmed by reading and hearing the signal. Second, determine the energy ratio around the band $[\phi^{\prime}(t)-0.2,\phi^{\prime}(t)+0.2]$ Hz; that is,

where $|V|^{2}$ is the spectrogram of the signal of length $T$ with the sampling rate $1/U$. Here, since our focus is the heart rate, we downsample the signal to $U=100$ Hz. Over the 92 recordings, the mean and standard deviation of $R$ determined from the original signal, the rectified signal, the RELU activated signal and $h_{0.1}$ are $0.03\%\pm 0.01\%$, $2.33\%\pm 0.56\%$, $1.64\%\pm 0.41\%$, and $1.25\%\pm 0.28\%$. See Figure 9 for an example of the overall procedure.