Engee documentation
Notebook

The spectrum of Prony

The Proni spectrum is a method for estimating the power spectrum of a random process based on an autoregressive model (AR). The method was proposed in 1969 by French electrical engineer Rene Proni. The basic idea is to approximate the time series as a sum of sinusoidal components and noise using autoregression parameters. Spectral analysis by the Proni method is particularly useful for isolating frequency components in noisy signals and is actively used in signal processing and radar systems.

Engee does not have built-in functions for constructing spectra using the Prony transform, but in Julia you can implement the algorithm yourself or use available libraries to solve similar problems.

The purpose of this demonstration is to implement such an algorithm. The following is an example of the implementation of an algorithm for estimating the power spectrum using the Proni method.

In [ ]: 
# Install the necessary packages
Pkg.add("LinearAlgebra")
using LinearAlgebra

Next, we implement the prony_spectrum function. This function performs a Prony conversion for a given signal. It takes three arguments: the signal itself (signal), the order of the model (order), and the sampling frequency (fs). The function returns two values: frequencies and their corresponding amplitudes, which represent the components of the signal spectrum. The process includes forming an autocorrelation matrix, solving a system of linear equations to obtain the coefficients of an autoregressive model, finding the roots of a characteristic equation, and calculating frequencies and amplitudes.

In [ ]: 
# The function for converting the
function prony_spectrum(signal::Vector{Float64}, order::Int, fs::Float64)
    N = length(signal)
    if order >= N
        error("The order of the model must be less than the length of the signal.")
    end

    # Formation of the autocorrelation vector
    X = zeros(Float64, N - order, order)
    for i in 1:(N - order)
        X[i, :] = signal[i:(i + order - 1)]
    end
    y = signal[(order + 1):N]

    a = X \ y # Solving the linear regression problem for the coefficients of the AR model
    a = vcat(1.0, -a)  # AR model coefficients
    roots_a = roots(a) # Finding the roots of a characteristic equation
    roots_a = roots_a[abs.(roots_a) .< 1.0] # We leave only the roots inside the unit circle

    # Calculating frequencies and amplitudes
    frequencies = angle.(roots_a) * fs / (2 * π)
    amplitudes = abs.(roots_a)

    return frequencies, amplitudes
end
Out[0]:
prony_spectrum (generic function with 1 method)

Next, we will set the noisy input signal.

In [ ]: 
fs = 1000.0  # Sampling rate
t = 0:1/fs:1-1/fs  # Timeline
# A signal consisting of several sinusoids
signal = sin.(2π * 50 * t) + 0.5 * sin.(2π * 120 * t)
# Adding noise
noisy_signal = signal + 0.05 * randn(length(signal))

plot(t, noisy_signal)
Out[0]:

Now let's calculate and plot the Proni spectrum.

In [ ]: 
order = 10  # The order of the model
frequencies, amplitudes = prony_spectrum(noisy_signal, order, fs)

# Visualization
plot(frequencies, amplitudes, seriestype=:stem, yaxis=:log, xlabel="Frequency (Hz)", ylabel="The amplitude",
     title="The spectrum of Prony", marker=:circle)
Out[0]:

Conclusion

In this example, we have demonstrated how to build a spectrum in Engee using a custom spectrum calculation function.