Frequency Estimation Method Using Sinusoidal Constraint Differential Equation (SCDE) and Python Implementation

Introduction

Conventional frequency estimation methods like FFT (Fast Fourier Transform) and autocorrelation functions are widely used. However, these methods share a common challenge: it is difficult to accurately estimate the frequency when multiple sounds are playing simultaneously from short signals (e.g., tens of milliseconds). FFT faces a trade-off between temporal and frequency resolution, and for short signals, frequency resolution deteriorates, often making it impossible to separate closely spaced frequency components.

In this article, I will introduce one approach to tackle this challenge: the Sinusoidal Constraint Differential Equation (SCDE) [1-3]. I will also present its Python implementation.

Relationship between Sine Waves and Differential Equations: Basics of SCDE

The core of SCDE lies in the simple properties of sine waves. The sine wave we are familiar with can be expressed in the following form:

$\begin{align*} x(t) = A \sin(\omega t + \phi), \end{align*}$

where $A$ is the amplitude, $t$ is time, $\omega$ is the angular frequency, and $\phi$ is the initial phase. Differentiating this sine wave twice with respect to time yields the following 2nd-order differential equation:

$\begin{align*} \frac{d^{2}}{dt^{2}} x(t) = -A\omega^{2} \sin(\omega t + \phi) = -\omega^{2} x(t) \end{align*}$

Rearranging the right side to the left side, we obtain the following 2nd-order differential equation:

$\begin{align*} \frac{d^{2}}{dt^{2}} x(t) + \omega^{2} x(t) = 0 \end{align*}$

This is the most basic form of the "Sinusoidal Constraint Differential Equation (SCDE)." What is noteworthy is that this equation does not depend on the amplitude $A$ or the initial phase $\phi$ . This suggests the possibility of extracting only the frequency $\omega$ directly from the signal.

In SCDE, we seek $\omega$ such that this differential equation "holds as true as possible," meaning we find $\omega$ that minimizes the squared residual of the given objective function. By differentiating this objective function with respect to $\omega$ and setting it to 0, a simple algebraic equation concerning $\Omega = \omega^{2}$ is obtained, and by finding its roots, the frequency can be estimated.

The specific form of the algebraic equation is shown below. First, $J(\omega)$ is rearranged into an expression containing $\Omega$ . To explicitly indicate that the objective function is a function of $\Omega$ , it is written as $J(\Omega)$ .

$\begin{align*} J(\Omega) &= \int_{\Gamma} \left[ \left(\frac{d^{2}}{dt^{2}} x(t)\right)^{2} + \Omega \left(\frac{d^{2}}{dt^{2}} x(t)\right) x(t) + \Omega^{2} x(t)^{2} \right] dt \\ &= \int_{\Gamma} \left(\frac{d^{2}}{dt^{2}} x(t)\right)^{2} dt + \Omega \int_{\Gamma} \left(\frac{d^{2}}{dt^{2}} x(t) \right) x(t) dt + \Omega^{2} \int_{\Gamma} x(t)^{2} dt \end{align*}$

Next, differentiate with respect to $\Omega$ .

$\begin{align*} \frac{J(\Omega)}{ d \Omega} = \int_{\Gamma} \left(\frac{d^{2}}{dt^{2}} x(t)\right) x(t) dt + 2\Omega \int_{\Gamma} x(t)^{2} dt \end{align*}$

Here, the covariance defined by the following equation is introduced.

$\begin{align*} S_{i, j} = \int_{\Gamma} \left[ \frac{d^{2i}}{dt^{2i}} x(t) \frac{d^{2j}}{dt^{2j}} x(t) \right] dt \end{align*}$

The algebraic equation for finding the optimal $\Omega$ is given by:

$\begin{align*} \frac{J(\Omega)}{ d \Omega} = 2 S_{0, 0} \Omega + S_{1, 0} = 0 \end{align*}$

This equation is a 1st-order equation for $\Omega$ , so it can be easily solved. The finally estimated frequency is:

$\begin{align*} \hat{\omega} = \sqrt{\Omega} = \sqrt{-\frac{ S_{1, 0} } { 2 S_{0, 0} }} \end{align*}$

Estimating Frequencies from Multiple Independent Sine Waves with SCDE

While basic SCDE is effective for a single sine wave, real sounds contain multiple frequency components. SCDE can be extended to handle sums of multiple independent sine waves. Generally, for $N$ sine waves, the signal is written as:

$\begin{align*} x(t) &= A_1 \cos (\omega_1 t + \phi_1) + A_2 \cos (\omega_2 t + \phi_2) + \cdots + A_N \cos (\omega_N t + \phi_N) \\ &= \sum_{i = 1}^{N} A_i \cos(\omega_i t + \phi_i) \end{align*}$

This signal satisfies a high-order differential equation obtained by chaining the differential operators that each sine wave satisfies:

$\begin{align*} \left(\frac{d^{2}}{dt^{2}} + \omega_1^{2}\right)\left(\frac{d^{2}}{dt^{2}} + \omega_2^{2}\right) \cdots \left(\frac{d^{2}}{dt^{2}} + \omega_N^{2}\right) x(t) = 0 \end{align*}$

Indeed, all terms from $A_1 \cos (\omega_1 t + \phi_1)$ to $A_N \cos (\omega_N t + \phi_N)$ satisfy the equation based on the above differential operator. Expanding the left hand side, it can be rearranged into the following form:

$\begin{align*} \alpha_{0} D^{2N} x(t) + \alpha_{1} D^{2N - 2} x(t) + \cdots + \alpha_{N-1} D^{2} x(t) + \alpha_{N} x(t) = 0 \end{align*}$

Here, $D^{k} = d^{k} / d t^{k}$ is introduced as a symbol. Each $a_n$ can be written as a fundamental symmetric polynomial of $\omega_1^{2}, \ldots, \omega_N^{2}$ . Conversely, $\omega_n$ can be obtained by using information from $\alpha_{0}$ to $\alpha_{N}$ (expressed as "information" with the meaning that they may not be obtained algebraically).

Similar to the case of a single sine wave, to find $\alpha_{n}$ such that the above equation holds as true as possible, we consider minimizing the objective function defined by squaring the left side:

$\begin{align*} J (\alpha_0, \alpha_1, \ldots, \alpha_N) = \int_{\Gamma} \left[ \sum_{n=0}^{N} \alpha_{n} D^{2(N-n)} x(t) \right]^{2} dt \end{align*}$

By solving the system of linear equations obtained from the extremum conditions for $\alpha_1, \ldots, \alpha_N$ , each $\alpha_{n}$ can be determined. The result of differentiating $J (\alpha_0, \ldots, \alpha_N)$ with respect to each $\alpha_1, \ldots, \alpha_N$ and setting it to 0 is:

$\begin{align*} S_{N,N-1} &+ \alpha_1 S_{N-1,N-1} + \dots + \alpha_N S_{0, N-1} = 0 \\ S_{N, N-2} &+ \alpha_1 S_{N-1,N-2} + \dots + \alpha_N S_{0, N-2} = 0 \\ &\vdots \nonumber \\ S_{N, 0} &+ \alpha_1 S_{N-1, 0} \hspace{5.1mm} + \dots + \alpha_N S_{0,0}\hspace{5.5mm} = 0 \end{align*}$

However, since $\alpha_{0}$ is always a constant from the expansion, the extremum condition for $\alpha_0$ is omitted.The definition of $S_{i, j}$ has already been stated, but it is reiterated here.

$\begin{align*} S_{i, j} = \int_{\Gamma} \left[ D^{2i} x(t) D^{2j} x(t) \right] dt \end{align*}$

Rearranging the equations slightly, we get:

$\begin{align*} \textbf{B} \textbf{a} = - \textbf{b} \end{align*}$

where $\textbf{a} = \left( \alpha_1, \alpha_2, \ldots, \alpha_N\right)^{\top}$ , $\textbf{b} = \left( S_{N,N-1}, S_{N,N-2} \ldots, S_{N, 0} \right)^{\top}$ , and

$\begin{align*} \textbf{B} = \begin{pmatrix} S_{N-1,N-1} & \cdots & S_{0,N-1} \\ \vdots & \ddots & \vdots \\ S_{N-1, 0} & \cdots & S_{0, 0} \end{pmatrix}. \end{align*}$

$\textbf{B}$ corresponds to the covariance matrix. The optimal are then given by:

$\begin{align*} \textbf{a} = - \textbf{B}^{-1} \textbf{b} \end{align*}$

Once the coefficients of the objective function are obtained, they are fundamental symmetric polynomials, $\omega_1^{2}, \ldots, \omega_N^{2}$ can be estimated as the solutions of the following $N$ -the order equation:

$\begin{align*} \Omega^{N} + (-1)^{1} \alpha_1 \Omega^{N-1} + (-1)^{2} \alpha_1 \Omega^{N-2} + \cdots + (-1)^{N-1} \alpha_N = 0 \end{align*}$

Algebraic solutions exist up to 4th-order equations, but 5th-order and higher equations cannot be solved algebraically. In the implementation, all equations are solved numerically.

Implementation

The Python code is available from the below. Enjoy!

A demonstration of Frequency Estimation Based on SCDE. · GitHub

Overview of the implemented functions:

apply_bandpass_filter: A function that applies a Butterworth bandpass filter.
get_kth_even_derivative : A function that calculates the $2i$ -th order derivative $D^{2i} x(t)$ . Numerical differentiation is implemented using the centered difference method.
calc_covmat_mn: A function that calculates each component of the covariance matrix $\textbf{B}$ . Since the signal length shortens with each numerical differentiation, trimming is necessary to calculate the sum of products (numerical integration) of high-order derivatives while considering the signal length.
solve_poly: A function that solves the $N$ -th order equation for frequency estimation.
estimate_frequencies_scde: A function that solves the system of linear equations for coefficients $\alpha_{n}$ using the results of covariance matrix calculations by calc_cov_mn. Subsequently, it solves the $N$ -th order equation based on these coefficients using solve_poly to estimate the final frequencies. In this experiment, $N=3$ .
main : A function for conducting verification experiments.

Experiments

Experiments were conducted to estimate the frequencies of a signal composed of three sine waves: 440 Hz, 460 Hz, and 480 Hz. The sampling frequency was 44.1 kHz. Signal lengths were of two types: 40 ms and 2 ms. For the 40 ms case, a pattern where white noise was added to the signal before estimation was also attempted. The noise frequency band at that time covered the three frequencies, from 430 Hz to 490 Hz. Noise bandwidth limitation was achieved using a Butterworth bandpass filter.

Experimental Resutls

The results of running the script are as follows.

--- Test without noise ---
True frequencies: [440. 460. 480.] Hz
Estimated frequencies: [439.92795114 459.91768358 479.90645314] Hz
Estimation errors: [-0.07204886 -0.08231642 -0.09354686] Hz

--- Test with band-limited noise ---
True frequencies: [440. 460. 480.] Hz
Estimated frequencies with noise: [439.86240043 459.95433428 479.8400905 ] Hz
Estimation errors: [-0.13759957 -0.04566572 -0.1599095 ] Hz

--- Test with very short signal duration (2ms) ---
True frequencies: [440. 460. 480.] Hz
Estimated frequencies for 2ms signal: [440.91214252 460.90604637 480.98345977] Hz
Estimation errors: [0.91214252 0.90604637 0.98345977] Hz

"Test without noise" shows the results of estimating a sum of three sine waves without adding noise. The estimation error is quite small, indicating successful pitch estimation.

"Test with band-limited noise" shows the results of estimation after adding band-limited white noise. Even with a noise level of 10dB, the estimation results are good.

"Test with very short signal duration" shows the results of estimation for a very short signal, 2ms long. Although the error increases compared to the 40ms case, the error remains around 1Hz.

Conclusion

In this article, I introduced the pitch estimation method using SCDE and its implementation. Based on simple verification experiments, pitch estimation seems to have been successful.

SCDE can also be applied to signals with harmonic structures such as speech, but each frequency is estimated independently, and the harmonic structure between frequencies is not considered. By imposing harmonic constraints on SCDE, further improvement in estimation accuracy can be expected. It has also been pointed out that SCDE has low robustness against noise. In Reference [2], a method based on Finite-Order Harmonic Constraint Differential Equation (FOHCDE), which introduces harmonic constraints into SCDE, is proposed. In the next article, I will introduce the FOHCDE.

References

[1] K. Yamada, Y. Masuyama, Y. Wakabayashi and N. Ono, "Simultaneous Frequency Estimation for Three or More Sinusoids Based on Sinusoidal Constraint Differential Equation," 2022 Asia-Pacific Signal and Information Processing Association Annual Summit and Conference (APSIPA ASC), pp. 975-978, 2022. https://ieeexplore.ieee.org/document/9980228

[2] K. Yamada, Y. Masuyama, K. Yamaoka and N. Ono, "Fundamental Frequency Estimation Based on Finite-Order Harmonic Constraint Differential Equation," 2023 Asia Pacific Signal and Information Processing Association Annual Summit and Conference (APSIPA ASC), pp. 868-872, 2023. https://ieeexplore.ieee.org/document/10317393

[3] K. Yamada, Y. Masuyama, K. Yamaoka, N. Ueno and N. Ono "Multiple Pitch Estimation Based on Finite-Order Harmonic Constraint Differential Equation," IEICE Technical Report, vol. 123, no. 401, pp. 315-320, 2024. https://ken.ieice.org/ken/paper/20240301scCv/