Engee documentation
Notebook

Spectrum Prony

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

Engee does not have inbuilt functions for spectral analysis using Prony transform, but you can implement the algorithm yourself in Julia or use available libraries for such tasks.

The purpose of this demonstration is to implement such an algorithm. Below is an example of implementation of the algorithm for power spectrum estimation using the Prony method.

In [ ]: 
# Установите необходимые пакеты
Pkg.add("LinearAlgebra")
using LinearAlgebra

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

In [ ]: 
# Функция для преобразования Прони
function prony_spectrum(signal::Vector{Float64}, order::Int, fs::Float64)
    N = length(signal)
    if order >= N
        error("Порядок модели должен быть меньше длины сигнала.")
    end

    # Формирование вектора автокорреляции
    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 # Решение задачи линейной регрессии для коэффициентов AR-модели
    a = vcat(1.0, -a)  # Коэффициенты AR модели
    roots_a = roots(a) # Нахождение корней характеристического уравнения
    roots_a = roots_a[abs.(roots_a) .< 1.0] # Оставляем только корни внутри единичной окружности

    # Рассчитываем частоты и амплитуды
    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 set the input noisy signal.

In [ ]: 
fs = 1000.0  # Частота дискретизации
t = 0:1/fs:1-1/fs  # Временная шкала
# Сигнал из нескольких синусоид
signal = sin.(2π * 50 * t) + 0.5 * sin.(2π * 120 * t)
# Добавление шума
noisy_signal = signal + 0.05 * randn(length(signal))

plot(t, noisy_signal)
Out[0]:

Now let's calculate and plot the Prony spectrum.

In [ ]: 
order = 10  # Порядок модели
frequencies, amplitudes = prony_spectrum(noisy_signal, order, fs)

# Визуализация
plot(frequencies, amplitudes, seriestype=:stem, yaxis=:log, xlabel="Частота (Гц)", ylabel="Амплитуда",
     title="Спектр Прони", marker=:circle)
Out[0]:

Conclusion

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