Multi-layered Cepstrum for Instantaneous Frequency Estimation

TL;DR

The paper introduces a multi-layered cepstrum method that recursively enhances fundamental frequency estimation, effectively handling noisy and contaminated signals in polyphonic music analysis.

Contribution

It presents a novel recursive cepstrum approach for improved instantaneous frequency estimation in challenging noise conditions.

Findings

Effective in high-pass filter noise scenarios

Improves MF0 saliency in polyphonic music

Demonstrates potential in real-world datasets

Abstract

We propose the multi-layered cepstrum (MLC) method to estimate multiple fundamental frequencies (MF0) of a signal under challenging contamination such as high-pass filter noise. Taking the operation of cepstrum (i.e., Fourier transform, filtering, and nonlinear activation) recursively, MLC is shown as an efficient method to enhance MF0 saliency in a step-by-step manner. Evaluation on a real-world polyphonic music dataset under both normal and low-fidelity conditions demonstrates the potential of MLC.

Figures

Tables

Table 1. Table 1 : Comparison of different layer combinations and testing method on the Bach10 dataset.

Layers	Brute Force									Greedy				SGD (10-fold CV)
$\{L,L+1\}$	${\gamma }_0$	${\gamma }_1$	${\gamma }_2$	${\gamma }_3$	${\gamma }_4$	${\gamma }_5$	P (%)	R (%)	F (%)	${\gamma }_i$	P (%)	R (%)	F (%)	P (%)	R (%)	F (%)
0th & 1st	0.3	1	x	x	x	x	77.72	83.51	80.51	0.33	78.05	83.57	80.71	77.67	78.43	78.05
1st & 2nd	0.3	0.5	1	x	x	x	80.55	87.20	83.74	0.48	81.05	85.23	83.09	79.35	88.10	83.50
2nd & 3rd	0.2	0.6	0.9	1	x	x	79.33	89.93	84.30	0.71	80.28	86.66	83.35	82.61	85.99	84.26
3rd & 4th	0.1	0.9	0.9	0.5	1	x	81.06	90.80	85.66	0.65	80.16	87.08	83.47	81.02	85.83	83.36
4th & 5th	0.1	0.9	0.9	0.7	0.8	1	81.14	90.67	85.64	0.78	80.20	87.27	83.58	81.80	86.22	83.95
5th & 6th	0.1	0.9	0.9	0.7	0.8	0.5	82.35	91.10	86.50	0.70	80.29	87.57	83.77	81.97	82.18	82.07

Equations

$${\mathbf{z}}^{(0)}={\sigma }^{(0)}\left(|\mathbf{Fx}|\right),{\mathbf{z}}^{(l)}={\sigma }^{(l)}\left({\mathbf{W}}^{(l)}\mathbf{F}{\mathbf{z}}^{(l-1)}\right),$$

$${\mathbf{Z}}^{(0)}:={\sigma }^{(0)}\left(\left|\mathbf{X}\right|\right),{\mathbf{Z}}^{(l)}:={\sigma }^{(l)}\left({\mathbf{W}}^{(l)}\mathbf{F}{\mathbf{Z}}^{(l-1)}\right),$$

$${\mathbf{Y}}^{(l_e,l_o)}[k,n]={\mathbf{Z}}^{(l_e)}[k,n]{\mathbf{Z}}^{(l_o)}\left[\lfloor N/k\rceil ,n\right]$$

Full text

**Index Terms— ** Cepstrum, multiple fundamental frequency estimation, time-frequency analysis.

1 Introduction

Multiple fundamental frequency (MF0) estimation is an essential problem in analyzing various multi-component signals including multi-talker speech [1, 2], polyphonic music [3], fetal electrocardiogram (FECG), photoplethysmographic (PPG) signals [4], and IoT sensor data [5], to name but a few. Numerous methods have been proposed to enhance the saliency of F0s in the signal representation, such as the autocorrelation function (ACF) [6], generalized cepstrum [7], nonnegative matrix factorization (NMF), MUltiple SIgnal Classification (MUSIC) [8], combined frequency and periodicity (CFP) [9, 10, 11], and deep learning approaches [12]. The generalized cepstrum, being one of the earliest proposed MF0 estimation methods [7], still gains attention in recent years, by showing its effectiveness in suppressing the unwanted harmonics and in localizing the F0 components by combining it with a spectrum [9, 10, 13]. A recent study further proved that a cepstrum can perfectly localize the instantaneous F0s of a signal with the intrinsic mode function (IMF) type [11], a general type of non-stationary oscillatory signals.

One issue which is relatively less discussed in MF0 estimation is the so-called convolutional noise in signals [14]. Convolutional noise is ubiquitous; it occurs when signals are transmitted through an intermediate channel or device. An MF0 algorithm usually fails in the presence of convolutional noise, for example, in the processing of audio signals played by a smartphone, as the response of a smartphone speaker typically behaves like a high-pass filter whose cutoff frequency is around 500 Hz [15, 9]. Such a high-pass filter suppresses the true F0 peaks, as similar to the case of the “missing fundamental” effects in low-pitched musical signals [16]. In such cases, to identify the true F0 peaks merely from the interleaved harmonics of multiple components is challenging.

In this paper, we generalize the notion of cepstrum to resolve the above-mentioned issue in MF0 estimation. The proposed multi-layered cepstrum (MLC) refines the saliency of the F0s in a recursive manner. In a nutshell, the MLC takes the operation of cepstrum (i.e., Fourier transform, filtering and nonlinear activation) repeatedly: for a time-domain input signal $\mathbf{x}\in\mathbb{R}^{N}$, $l\geq 0$, the $l^{\textrm{th}}$-layer output of MLC is:

[TABLE]

where $\mathbf{F}\in\mathbb{C}^{N\times N}$ is the $N$-point discrete Fourier transform (DFT) matrix, $\mathbf{W}\in\mathbb{R}^{N\times N}$ is a high-pass filter, and $\sigma(\cdot)$ is an element-wise power-scaled rectification unit (see Section 3). In the following sections, we will explain how the idea of MLC is motivated, how it works in MF0 estimation, and its structural resemblance to multi-layered perceptrons (MLP). We will then focus on combining two MLCs lying in both the frequency and time domains to achieve MF0 estimation by employing the CFP approach. Experiments on both synthetic and real data demonstrates the effectiveness of MLC in MF0 estimation in both normal and contaminated conditions, and indicate positively the fact that MLC performs better than other methods based on a ‘shallow’ DFT.

2 Motivation of MLC

The cepstrum $\hat{\mathbf{z}}:=\mathbf{F}^{-1}\sigma(\mathbf{F}\mathbf{x})$ of $\mathbf{x}$ [17, 18] is a classic method to estimate the F0 of $\mathbf{x}$. A logarithm cepstrum employing $\sigma(x):=\log x$ is commonly used in single F0 estimation, while the generalized cepstrum (GC) employing $\sigma(x):=x^{\gamma}$, $0<\gamma\leq 2$, is found more applicable in MF0 estimation, probably because it is numerically more stable than using the logarithm function and is more robust to noise [19, 20]. Cepstrum-based MF0 estimation is based on an assumption that the fast-varying components in a spectrum (e.g., harmonic sequences) are relevant to F0s, while the slow-varying ones (e.g., spectral envelope) are irrelevant to F0s. Saliency of F0s can be enhanced in the quefrency111Since the quefrency has the same unit as time, in this paper, the terms quefrency and time are used interchangeably. domain by taking another Fourier transform and a high-pass filter with a cutoff quefrency at $q_{c}$ on a nonlinearly-stretched spectrum, where the nonlinear operation is to fit the magnitude of the spectrum to humans’ perception scale.

To motivate the MLC, we revise the above-mentioned assumption to that 1) the fast-varying and periodic components in a spectrum are relevant to F0s while others are not, and that 2) such periodic components still remain periodic (and therefore fast-varying) after a Fourier transform, which is a basic property of DFT. These assumptions allow us to consider repeating the operation of cepstrum on $\mathbf{x}$ to iteratively discard the slow-varying or aperiodic components and to purify the F0-relevant components, which should remain after a number of filtering processes. This idea has been preliminarily verified in [21], where the generalized cepstrum of spectrum (GCoS), being the $2^{\text{nd}}$layer output of the MLC, behaves as a more succinct feature indicating the F0 peaks than a spectrum. This paper is an extension of [21] to a multi-layered time-frequency analysis tool.

Admittedly, the above-mentioned idea is also motivated by deep learning, which catches wide attention recently in MF0 estimation [22, 12, 23, 24, 25]. Not only achieving state-of-the-art results, deep learning also provides inspiring observation when solving the MF0 estimation problem. For example, results in [12] show that the learned networks for melody tracking resemble traditional pitch detection methods: the first layer behaves like a spectral analyzer, while the second layer behaves like a comb filter. It is again intriguing to consider a multi-layered network for MF0 estimation, where the well-known pitch detection blocks such as the Fourier transform and filters are taken directly as subnetworks.

3 The MLC-CFP method

We formulate the MLC of a non-stationary signal $\mathbf{x}$ in the time-frequency or time-quefrency domain. Given a window function $\mathbf{h}\in\mathbb{R}^{N}$ and hop size $H$, the short-time Fourier transform (STFT) of $\mathbf{x}$, $\mathbf{X}:=\mathbf{X}[k,n]$, is represented as $\mathbf{X}:=\mathbf{X}[k,n]=\sum^{N-1}_{m=0}\mathbf{x}[m+nH]\mathbf{h}[m]e^{-\frac{j2\pi km}{N}}$, where $\mathbf{X}\in\mathbb{C}^{N\times M}$, $M$ is the number of frames, $k$ is the frequency index, and $n$ is the time index. Denote $|\mathbf{X}|$ the element-wise absolute value of $\mathbf{X}$. The $l^{\textrm{th}}$-layer short-time MLC is:

[TABLE]

where $\mathbf{Z}^{(l)}\in\mathbb{R}^{N\times M}$, $\mathbf{Z}^{(l)}:=\mathbf{Z}^{(l)}[k,n]$ for even $l$, $\mathbf{Z}^{(l)}:=\mathbf{Z}^{(l)}[q,n]$ for odd $l$, and $q$ is the quefrency index. $\sigma^{(l)}$ is a nonlinear, element-wise power function such that for $\gamma_{l}>0$, $\sigma^{(l)}\left(x\right)=x^{\gamma_{l}}$ for $x>0$, and $\sigma^{(l)}\left(x\right)=0$ for $x\leq 0$. With the $N$-point DFT matrix $\mathbf{F}$, the frequency resolution of $\mathbf{Z}^{(l)}[k,n]$ is $f_{s}/N$ per bin, and the time resolution of $\mathbf{Z}^{(l)}[q,n]$ is then $1/f_{s}$ per bin. $\mathbf{W}^{(i)}\in\mathbb{R}^{N\times N}$ represents a high-pass filter; it filters out the elements lying below a cutoff frequency (quefrency) $i_{c}$. More precisely, $\mathbf{W}^{(l)}$ is a diagonal matrix such that $\mathbf{W}^{(l)}[i,i]=1$ when $i>i_{c}$ and $i<N-i_{c}$, while $\mathbf{W}^{(l)}[i,i]=0$ otherwise. For even $l$, $i_{c}:=k_{c}$ defines the cutoff frequency at $f_{c}=k_{c}f_{s}/N$, and for odd $l$, $i_{c}:=n_{c}$ defines the cutoff quefrency at $q_{c}=n_{c}/f_{s}$. Also note that since $\mathbf{Z}^{(l)}[:,n]$ is even for $l>2$, we use $\mathbf{F}$ instead of $\mathbf{F}^{-1}$ when transforming from the frequency domain to the time domain. Many of the basic pitch detection functions are essentially special cases of (1). For example, $\mathbf{Z}^{(0)}$ is the STFT of $\mathbf{x}$, and $\mathbf{Z}^{(1)}$ is known as the autocorrelation function (ACF) for $\gamma_{0}=2$ and $\gamma_{1}=1$. (1) also resembles a DNN structure with zero bias vector: $\mathbf{W}^{(l)}\mathbf{F}$ resembles the fully-connected layer and $\sigma^{(l)}$ resembles the activation function.

To enhance the true F0 components and suppress the high-order harmonics, consider the last two output layers, $\mathbf{Z}^{(l_{e})}$ and $\mathbf{Z}^{(l_{o})}$, of an MLC, where $l_{e}$ is even and $l_{o}$ is odd. $\mathbf{Z}^{(l_{e})}$ is in the frequency domain and $\mathbf{Z}^{(l_{o})}$ is in the quefrency domain. Based on the CFP approach, we utilize the duality between them to enhance the saliency of F0 by employing a nonlinear mapping to $\mathbf{Z}^{(l_{o})}$ from the quefrency domain to frequency domain, and then multiplying it by $\mathbf{Z}^{(l_{e})}$ [13, 9, 10, 11]:

[TABLE]

where the time-quefrency representation $\mathbf{Z}^{(l_{o})}[k,n]$ is nonlinearly mapped into the frequency domain, and $\lceil\cdot\rfloor$ is the rounding function. Since $\mathbf{Z}^{(l_{o})}$ suppresses the harmonic peaks in $\mathbf{Z}^{(l_{e})}$ while $\mathbf{Z}^{(l_{e})}$ suppresses the sub-harmonic peaks in $\mathbf{Z}^{(l_{o})}$, the resulting CFP representation $\mathbf{Y}^{(l_{e},l_{o})}:=\mathbf{Y}^{(l_{e},l_{o})}[k,n]$ has only true F0 peaks enhanced. The nonlinear mapping of $\mathbf{Z}^{(l_{o})}$ in (2) can be implemented by a filterbank following a custom-defined indexes of $k$; see [10, 21] for details. We name the above process as the MLC-CFP algorithm hereafter.

4 Numerical results

We illustrate the effectiveness of the MLC-CFP method in emphasizing F0 peaks of a simulated signal contaminated with both additive and convolutional noise. Given two components, $\mathbf{x}_{1}$ and $\mathbf{x}_{2}$, both are sampled at $f_{s}=1$ kHz. $\mathbf{x}_{1}$ is a square wave having 20% duty cycle, and has a constant F0 at $2$ Hz. $\mathbf{x}_{2}$ is a frequency-modulated sawtooth wave with positive ramp and vertical drop, and has time-varying F0 represented as $f_{0}[n]=2.5+\cos(2\pi n/10f_{s})$. To simulate the effect of convolutional noise, $\mathbf{x}_{1}$ is passed through a 10th-order Butterworth low-pass filter, and $\mathbf{x}_{2}$ is passed through a 10th-order Butterworth high-pass filter. Both filters have a cutoff frequency at 10 Hz. We consider the 2-layer MLC of two sample signals: the first one is $\mathbf{x}=\mathbf{x}_{1}+\mathbf{x}_{2}$, which is with convolutional noise only, while the second one, $\tilde{\mathbf{x}}=\mathbf{x}_{1}+\mathbf{x}_{2}+\mathbf{n}+\mathbf{d}$, is further contaminated with additive pink noise and impulse noise, where $\mathbf{n}$ is additive pink noise making the SNR of $\tilde{\mathbf{x}}$ be 10 dB and $\mathbf{d}$ is an impulse at the 80th second.

Fig. 1 shows the STFT ($\mathbf{Z}^{(0)}$), GC ($\mathbf{Z}^{(1)}$), GCoS ($\mathbf{Z}^{(2)}$), and CFP representation ($\mathbf{Y}^{(2,1)}$) of $\mathbf{x}$ and $\hat{\mathbf{x}}$, respectively. The parameters $[\gamma_{0},\gamma_{1},\gamma_{2}]=[0.24,0.6,1]$ are used. For the STFT of both samples, the F0 peaks and the low-order harmonics of $\mathbf{x}_{2}$ are all unseen, and only the first 5 multiples of $\mathbf{x}_{1}$ are visible due to low-pass filtering. To capture the true F0 of $\mathbf{x}_{2}$ from the STFT is very challenging. However, such information can be found in the GC: although high-pass filtering noise attenuates the true F0 component in the spectrum, it does not change the period of the signal. Such a phenomenon can also be seen in the case with pink noise. The true F0 trajectory then appears in the GCoS with high resolution; this is because the nonlinear scaling of the GC provides the mechanism to down-mixing the frequency when it is transformed back to the frequency domain. Interestingly, we can see that although the true F0 of $\mathbf{x}_{2}$ is missed, the resulting F0 trajectory of $\mathbf{x}_{2}$ has resolution even higher than the one of $\mathbf{x}_{1}$. Since most of the harmonic peaks of $\mathbf{x}_{1}$ are eliminated, its spectrum cannot support better resolution when transformed to the quefrency domain. Finally, the CFP representation computed through (2) clearly indicates the true F0 trajectories of $\mathbf{x}_{1}$ and $\mathbf{x}_{2}$ with less interference.

5 Real-data Experiments

We apply the MLC-CFP method for MF0 estimation of polyphonic music data, and investigate the effects of 1) the number of layers, 2) critical parameters $\gamma_{i}$ of each layer, and 3) the presence of convolutional noise. We used the Bach10 dataset, containing ten quartets of four different instruments [26] for the evaluation of multi-pitch estimation (MPE). To see the behavior of the algorithm in a straightforward way, we consider three methods to optimize parameters $\{\gamma_{i}\}^{L}_{i=0}$ for an $L$-layer MLC-CFP algorithm combining $\mathbf{Z}^{(L-1)}$ and $\mathbf{Z}^{(L)}$:

1. the brute-force method that searches for the optimal parameters exhaustively for each $\gamma_{i}$. The search range of every $\gamma_{i}$ is from 0.1 to 0.9 with step size of 0.1.

2. the greedy method that searches for the optimal value of $\gamma_{l}$ based on the optimal values of $\{\gamma_{i}\}^{l-1}_{i=0}$. The search range is from 0.01 to 0.99 with step size of 0.01, since its light computation allows us to use finer grids. For the post-processing of the brute-force and the greedy methods, we use the CFP-based algorithm reported in [21], as a simplified version of [9] without too detailed hand-crafted rules.

3. the stochastic gradient descent (SGD) method that minimizes the binary cross-entropy between $\mathbf{Y}^{(L)}$ (with 88 log-frequency filterbanks and a sigmoid function as the output layer) and the ground-truth piano roll by performing SGD over every $\gamma_{i}$. We apply non-repeated 10-fold cross-validation on the dataset: for each time one music piece is chosen as the test data and others as training data. The resulting F-scores are obtained by summing all the counts of the true positives, false positives, and false negatives. We initialize $\gamma_{0}$ to 0.24, $\gamma_{1}$ to 0.6, and all other $\gamma_{i}$s to 1, and use the SGD optimizer with 0.1 learning rate, a mini-batch size of 256, and train the network up to 40 epochs.

For other parameters, we compute 7939-point DFT using Blackman-Harris window, the hop size is 10 ms, the cutoff frequency $f_{c}$ is 27.5 Hz (frequency of A0), and the cutoff quefrency $q_{c}$ is 0.24 ms (period of C8). Since there are many task-dependent parameters, all the source codes will be announce after this paper is accepted for reproducibility.

Table 1 lists the optimal precision (P), recall (R), F-score (F), along with $\gamma_{i}$ for every pairs $\left[\mathbf{Z}^{(L-1)},\mathbf{Z}^{(L)}\right]$ for all optimization approaches. When using the brute-force method, the resulting P, R and F all increase consistently when the number of layers increases. Note that the recall increases dramatically with higher-order layer because the weak F0 peaks are enhanced, similar to the simulation example in Fig. 1. When using the greedy method, we observe that although the performance is sub-optimal, it still increases as $L$ increases. When using gradient descent, there is no consistent improvement in F-scores related to the number of layers, possibly because the hand-crafted rules such as harmonic selection [9] could not be applied to this model. However, it still gives competitive performance for $L=3$. More specifically, when compared to other MF0 estimation methods, all the three methods in our discussion outperform the state-of-the-art matrix factorization approaches such as constrained NMF and probabilistic latent component analysis (PLCA) (see Table III in [9] for more details). The brute-force method also outperform the CFP result reported in [9], in which the parameters for the hand-crafted rules are also based on brute-force searching.

Fig. 2 shows the resulting F-scores on the Bach10 dataset with high-pass degradation, which is simulated by a 4th-order Butterworth high-pass filters with cutoff frequency from 10 to 1000 Hz. The parameters used here are those found from the brute force method; see Table 1. Fig. 2 clearly shows that with more layers the results are better: even under a severe cutoff frequency at 1 kHz, a 6-layer MLC still gives 72.74% F-score, 25% higher than the one of a 1-layer MLC. Even more interestingly, results in some degradation conditions are better than those in the normal condition. For example, pairs $\mathbf{Z}^{(5)}$ and $\mathbf{Z}^{(6)}$ achieve 87.10% F-score with cutoff frequency at 100 Hz; this is probably because the simulated degradation suppresses the low-frequency noise incidentally.

To summarize, our experiments show that 1) in the MLC-CFP algorithm, a deep DFT do performs better than a shallow DFT, and 2) MLC is highly robust to high-pass noise.

6 Discussion and conclusion

We have introduced the MLC, explained its physical insight in signal processing, and verified its potential in detecting multiple F0s through simulation and real-world data under normal and noisy conditions. We also emphasize that the MLC is not only useful in the signals contaminated with high-pass convolutional noise, but also useful in all signals which inherently have weak F0 magnitude. For example, the F0 of bassoon in the Bach10 dataset is usually weaker than other instruments (violin, clarinet and saxophone), but experiments show that this issue is addressed by the MLC. From the result in the filter test we also suggest that in practice using MLC with a little high-pass filter can get better result.

One risk of MLC is that it may generate fake F0 peaks. For example, two F0s at $a$ Hz and $b$ Hz may produces a cross term at $|a-b|$ Hz due to effect of nonlinear activation. This cross term cannot be suppressed and may be even enhanced because it appears in both the time-domain and frequency-domain representations. To restrict the strength of cross terms, the results indicate a possible strategy: setting a small $\gamma_{0}$ and larger $\gamma_{i}$s at the higher layers. Our pilot study also shows that even fake cross terms are produced, they are mostly weak unless all $\gamma_{i}$s approach zero.

Developing a more efficient way to optimize the parameters efficiently on a large-scale dataset with complicated sources and harmonic structures will be another focus in the future. From the current finding, we summarize that MLC-CFP works well in general and is relatively easy to optimize for $L\leq 3$, but a great potential exists in a deeper MLC.