Engee 文档
Notebook

Prony的光谱

Proni频谱是一种基于自回归模型（AR）估计随机过程功率谱的方法。 该方法由法国电气工程师Rene Proni于1969年提出。 基本思想是使用自回归参数将时间序列近似为正弦分量和噪声的总和。 Proni方法的频谱分析对于隔离噪声信号中的频率分量特别有用，并积极用于信号处理和雷达系统。

Engee没有使用Prony变换构建谱的内置函数，但在Julia中，您可以自己实现算法或使用可用的库来解决类似问题。

本演示的目的是实现这样的算法。 以下是使用Proni方法估计功率谱的算法的实现示例。

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

接下来，我们实现prony_spectrum函数。 该函数对给定信号执行Prony转换。 它需要三个参数：信号本身（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]:

现在让我们执行Proni谱的计算和构建。

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中构建频谱。