Сравнение MATLAB и Engee в задачах спектрального и частотно-временного анализа¶

Открыть пример в Engee

Сравнивая MATLAB и Engee в областях спектрального анализа и частотно-временного анализа, рассмотрим их возможности и подходы для двух задач: прямое и обратное преобразование Фурье и частотно-временной анализ.

Зададим случайный входной сигнал.

Pkg.add(["DSP"])

x_in = randn(300)
plot(x_in)

Спектральное преобразование Фурье (FFT и IFFT)¶

В Engee используются функции из пакета FFTW:

fft: прямое быстрое преобразование Фурье.
ifft: обратное преобразование Фурье.

using FFTW

# Прямое преобразование Фурье
X = fft(x_in)

# Обратное преобразование Фурье
x_e = ifft(X)

300-element Vector{ComplexF64}:
   0.4463051717889283 + 5.921189464667501e-17im
   0.8855411237159746 - 6.706442469221183e-17im
  0.07243544230012243 + 2.1699078972421766e-17im
   0.7807779029413544 + 4.304182606664186e-18im
  -0.7744411506274916 - 1.0928995693991658e-16im
   -1.228942424222238 + 1.4683980316522056e-16im
   0.1428703752318775 + 2.967057684964929e-16im
  -0.5864750847294308 + 9.009751805455231e-17im
  0.05285563932230948 + 2.0991548884401417e-16im
  -0.7790856945787832 + 1.2582026155502825e-16im
   0.7670042670857968 + 2.326449445437145e-16im
   0.9362308439525574 + 3.904630245671109e-17im
   1.2082983266432845 + 1.305234552263469e-17im
                      ⋮
  -1.0917404944568745 + 4.5493154035744685e-17im
  -0.2219188971922492 - 9.021592475224685e-17im
 -0.20553174827279228 - 8.521592275351213e-17im
  -0.9086379856526511 - 3.541162156791164e-18im
  -0.9905704238868754 - 9.622342670683686e-17im
   0.5614164972979249 - 9.672808704078372e-17im
   0.2281289591141986 + 3.328797311281298e-16im
  -0.6088203503665865 - 2.113291501585097e-16im
  0.47762193437252043 - 6.470451606116849e-17im
  -0.9614319626845755 - 2.4417912440022187e-16im
   0.7503811721918421 - 1.5024680408776403e-17im
  -0.5137131256364178 + 1.609145981158678e-17im

В MATLAB имеются встроенные функции для быстрого преобразования Фурье:

fft: прямое быстрое преобразование Фурье.
ifft: обратное быстрое преобразование Фурье.

using MATLAB

# Прямое преобразование Фурье
mat"""X = fft($(x_in));"""

# Обратное преобразование Фурье
x_m=mat"""ifft(X);"""

300-element Vector{Float64}:
  0.4463051717889284
  0.8855411237159747
  0.07243544230012253
  0.7807779029413541
 -0.7744411506274914
 -1.2289424242222382
  0.142870375231878
 -0.5864750847294304
  0.0528556393223091
 -0.7790856945787833
  0.7670042670857973
  0.9362308439525577
  1.2082983266432852
  ⋮
 -1.091740494456875
 -0.22191889719224908
 -0.20553174827279197
 -0.9086379856526511
 -0.9905704238868758
  0.5614164972979251
  0.2281289591141988
 -0.6088203503665864
  0.47762193437252043
 -0.9614319626845759
  0.7503811721918426
 -0.5137131256364174

Сравним результаты и рассчитаем погрешность.

plot(x_in)
plot!(real(x_e))
plot!(x_m)

error_e = abs.(x_in-real(x_e))
error_m = abs.(x_in-x_m)
println("Средняя ошибка Engee: $(sum(error_e)/length(error_e))")
println("Средняя ошибка MATLAB: $(sum(error_m)/length(error_m))")

Средняя ошибка Engee: 1.945579114481788e-16
Средняя ошибка MATLAB: 2.279629031839055e-16

В обоих средах результат прямого преобразования Фурье представляет собой комплексные числа, а обратное преобразование возвращает исходный сигнал с небольшой погрешностью, связанной с численным округлением. Таким образом, вы можете видеть, что синтаксис и логика кода идентичны.

Частотно-временное преобразование Фурье¶

В Engee под такие преобразования нет стандартных функций, и их необходимо реализовывать вручную.

using DSP
# Функция для генерации чирп-сигнала
function chirp(t::AbstractVector{T}, f0::T, t1::T, k::T) where T<:Real
    return cos.(2π * (f0 .* t .+ 0.5 * k .* t.^2))
end

# Окно Хэмминга
function hamming(N::Int)
    n = 0:N-1
    return 0.54 .- 0.46 .* cos.(2π * n / (N - 1))
end

l=500 # Размер окна
window=vcat(hamming(l),zeros(length(x)-l))

x = chirp(0:0.001:1, 0.0, 1.0, 100.0)  # Чирп-сигнал
P = periodogram(x; fs=1000, window=window, nfft=length(x))

# Визуализация спектра
plot(freq(P), 20*log10.(power(P)), xlabel="Frequency (Hz)", ylabel="Power Spectrum (dB)", title="Power Spectrum", linewidth=2)

MATLAB имеет функцию pspectrum, которая позволяет выполнять частотно-временной анализ сигнала. pspectrum автоматически выбирает подходящий метод анализа. Поддерживает различные визуализации (спектрограммы, плотности мощности).

mat"""
% Частотно-временной анализ
x = chirp(0:0.001:1, 0, 1, 100);
pspectrum(x, 1000); % Частотно-временной спектр
saveas(gcf, 'frequency_time_spectrum.jpg'); % Сохраняет текущую фигуру как JPG
"""
load( "$(@__DIR__)/frequency_time_spectrum.jpg" )

Исходя из графиков можно заметить, что основные тенденции поведения сигнала отражены в Engee. Для более детального построения графиков необходимо дополнительно реализовать более детальный расчёт окна для функции периодограммы.

Также давайте сравним корректность сгенерированных входных данных в MATLAB и Engee.

plot(mat"chirp(0:0.001:1, 0, 1, 100)"[:])
plot!(chirp(0:0.001:1, 0.0, 1.0, 100.0))
title!("Sum_error: $(sum(mat"chirp(0:0.001:1, 0, 1, 100)"[:]-chirp(0:0.001:1, 0.0, 1.0, 100.0)))")

Как мы видим, входные сигналы идентичны.

Вывод¶

В MATLAB есть встроенные функции для выполнения этих операций. В Engee под эту задачу можно использовать сторонние библиотеки, такие, как DSP или FFTW для более сложных операций с сигналами. Также важно учесть, что Engee даёт больше гибкости в том, как можно реализовать преобразования. Чтобы упростить выполнение этих операций для реализации ваших задач в Engee, используйте готовые подходы, описанные в данной демонстрации.

Ещё один аспект в сравнении возможностей двух сред говорит о том, что MATLAB оптимизирован под работу с матрицами и сигналами. Однако Engee в целом быстрее в вычислениях, если код написан с учетом особенностей этой среды, и это даёт существенный прирост по скорости выполнения и тестирования алгоримах в больших системах.