Engee 文档
Notebook

光谱普罗尼

Prony 频谱是一种基于自回归 (AR) 模型的随机过程功率谱估计方法。该方法由法国电气工程师 René Prony 于 1969 年提出。其基本思想是利用自回归参数将时间序列近似为正弦分量和噪声之和。普罗尼频谱分析法尤其适用于提取噪声信号中的频率成分，在信号处理和雷达系统中得到了广泛应用。

Engee 没有使用 Prony 变换进行频谱分析的内置函数，但您可以在 Julia 中自行实现该算法，或使用可用库完成此类任务。

本演示的目的就是实现这种算法。下面是使用 Prony 方法实现功率谱估计算法的示例。

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

接下来，我们实现 prony_spectrum 函数。该函数对给定信号进行普罗尼变换。它需要三个参数：信号本身（signal）、模型阶次（order）和采样频率（fs）。函数返回两个量：频率及其相应的幅度，它们代表信号频谱的组成部分。这一过程包括形成自相关矩阵、求解线性方程组以获得自回归模型系数、找到特征方程的根并计算出频率和振幅。

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)

接下来，我们设置输入的噪声信号。

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]:

现在，让我们计算并绘制普罗尼频谱图。

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]:

结论

在本例中，我们演示了如何在Engee中使用自定义频谱计算函数构建 Prony 频谱。