Трендовая и гармоническая составляющие

Автор
avatar-yurevyurev
Notebook

Трендовая и гармоническая составляющие (разложение Фурье)

Анализ временных рядов — важный инструмент исследования различных процессов, будь то экологические измерения, технические характеристики оборудования или показатели жизнедеятельности биологического организма.

Один из ключевых этапов обработки временного ряда заключается в выявлении закономерностей, трендов и циклических компонент, скрытых в измеренных данных. Эти компоненты позволяют лучше понимать динамику исследуемых явлений и делать обоснованные прогнозы.

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

  1. Парсинг данных: чтение сырых данных из файла и преобразование их в удобный формат для дальнейшей обработки.

  2. Добавление временной шкалы: расчет секундной отметки каждого наблюдения относительно начала измерений.

  3. Выделение тренда: выявление долгосрочной тенденции изменения уровня CO путём усреднения значений в скользящем окне и последующей аппроксимацией полиномом третьей степени.

  4. Спектральный анализ: применение быстрого преобразования Фурье (БПФ) для выделения значимых частот колебаний.

  5. Реконструкция сигнала: восстановление основных гармонических компонентов, формирующих структуру процесса после удаления тренда.

  6. Итоговая визуализация: сравнение восстановленного сигнала с исходными данными и оценка качества предложенного метода.

В этой работе мы продемонстрируем возможности аналитического подхода для выявления характерных особенностей временных рядов на примере реальных экспериментальных данных.

Описание реализации алгоритма

Основные этапы алгоритма

1. Подготовка данных

Для подготовки данных используется функция parse_line, позволяющая разбивать строку на отдельные числовые значения и приводить их к соответствующим типам данных. После парсинга данные загружаются в таблицу типа DataFrame.

Затем добавляется временная метка в формате секунд, рассчитанная относительно первого замера.

2. Определение тренда

Используется два способа определения тренда:

  • среднее значение в скользящем окне фиксированного размера,
  • аппроксимация полиномом третьей степени.

Первый метод удобен своей простотой и позволяет визуально оценить общий тренд, второй же даёт возможность точной параметрической модели.

3. Спектральный анализ

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

4. Реконструкция сигнала

На основе выделенных гармоник восстанавливается приближённый сигнал, который накладывается на исходную кривую, и тренд для оценки точности восстановления.

In [ ]:
# Функция для парсинга одной строки
function parse_line(line)
    # Заменяем запятые на точки для корректного парсинга чисел
    line = replace(line, "," => ".")
    # Разделяем по пробелам/табам и фильтруем пустые элементы
    elements = filter(!isempty, split(line, r"\s+"))
    # Преобразуем элементы в нужные типы
    co = parse(Float64, elements[1])
    ch = parse(Float64, elements[2])
    co2 = parse(Float64, elements[3])
    o2 = parse(Float64, elements[4])
    lamb = parse(Float64, elements[5])
    n = parse(Int, elements[6])
    nox = parse(Int, elements[7])
    tmas = parse(Float64, elements[8])
    time = Time(elements[9])
    return (CO=co, CH=ch, CO2=co2, O2=o2, lamb=lamb, n=n, NOx=nox, Tmas=tmas, Time=time)
end
Out[0]:
parse_line (generic function with 1 method)
In [ ]:
using DataFrames
# Чтение данных
data_lines = readlines("$(@__DIR__)/data.txt")[2:end]  # Пропускаем заголовок
parsed_data = [parse_line(line) for line in data_lines] # Парсим все строки
df = DataFrame(parsed_data) # Преобразуем в DataFrame
println(first(df, 5)) # Выводим первые 5 строк для проверки
println()
using Dates, FFTW, Statistics
# Добавим колонку с временем в секундах
start_time = df.Time[1]
df[!, :Seconds] = [Dates.value(t - start_time) / 1e9 for t in df.Time]  # переводим наносекунды в секунды
println(first(df, 5)) # Выводим первые 5 строк для проверки
5×9 DataFrame
 Row │ CO       CH       CO2      O2       lamb     n      NOx    Tmas     Time
     │ Float64  Float64  Float64  Float64  Float64  Int64  Int64  Float64  Time
─────┼──────────────────────────────────────────────────────────────────────────────
   1 │    0.64    113.0    13.83     0.75    1.012    820     13      0.0  00:00:00
   2 │    0.64    112.0    13.83     0.75    1.012    770    185      0.0  15:56:25
   3 │    0.64    112.0    13.83     0.75    1.012    790      0      0.0  15:56:26
   4 │    0.64    112.0    13.83     0.75    1.012    800      0      0.0  15:56:27
   5 │    0.64    112.0    13.84     0.74    1.012    990      0      0.0  15:56:28

5×10 DataFrame
 Row │ CO       CH       CO2      O2       lamb     n      NOx    Tmas     Time      Seconds
     │ Float64  Float64  Float64  Float64  Float64  Int64  Int64  Float64  Time      Float64
─────┼───────────────────────────────────────────────────────────────────────────────────────
   1 │    0.64    113.0    13.83     0.75    1.012    820     13      0.0  00:00:00      0.0
   2 │    0.64    112.0    13.83     0.75    1.012    770    185      0.0  15:56:25  57385.0
   3 │    0.64    112.0    13.83     0.75    1.012    790      0      0.0  15:56:26  57386.0
   4 │    0.64    112.0    13.83     0.75    1.012    800      0      0.0  15:56:27  57387.0
   5 │    0.64    112.0    13.84     0.74    1.012    990      0      0.0  15:56:28  57388.0

Далее переходим к выбору переменной для анализа. Возьмем переменную CO (можно заменить на любую другую):

In [ ]:
i = 1 # Номер столбца для анализа
y = df[:,i]
column_names = names(df) # Получаем имена столбцов
println("Имя столбца для анализа: $(column_names[i])")
t = df.Seconds.-(15*60*60+56*60); t[1] = 0;
println("t: $(t[1:5])") # Выводим первые 5 строк для проверки
Имя столбца для анализа: CO
t: [0.0, 25.0, 26.0, 27.0, 28.0]

Далее построим трендовую составляющую. Используем скользящее среднее для выделения тренда:

In [ ]:
window_size = 15  # размер окна для сглаживания
trend = [mean(y[max(1, i-window_size):min(end, i+window_size)]) for i in 1:length(y)]

# График исходных данных и тренда
plot(t, y, label="Исходные данные", xlabel="Время (с)", ylabel="CO", legend=:topleft)
plot!(t, trend, label="Тренд", linewidth=2)
Out[0]:

Для более точного тренда можно использовать полиномиальную аппроксимацию вместо скользящего среднего:

In [ ]:
using Polynomials
trend_coef = fit(t, y, 3)  # полином третьей степени
trend = trend_coef.(t)

# График исходных данных и тренда
plot(t, y, label="Исходные данные", xlabel="Время (с)", ylabel="CO", legend=:topleft)
plot!(t, trend, label="Тренд", linewidth=2)
Out[0]:

Произведем разложение Фурье для гармонической составляющей.

In [ ]:
# Удаляем тренд для анализа периодических составляющих
detrended = y .- trend

# Вычисляем БПФ
n = length(detrended)
freq = fftfreq(n, 1/mean(diff(t)))  # частоты в Гц
fft_vals = fft(detrended)

# Амплитуды и фазы
amplitude = abs.(fft_vals) ./ n * 2
phase = angle.(fft_vals)

# Оставляем только положительные частоты (симметрия БПФ)
pos_freq = freq[1:div(n,2)+1]
pos_ampl = amplitude[1:div(n,2)+1]

# График амплитудного спектра
plot(pos_freq, pos_ampl, xlabel="Частота (Гц)", ylabel="Амплитуда", 
     title="Амплитудный спектр", legend=false)
Out[0]:

Произведем восстановление основных гармоник. Выберем несколько наиболее значимых частот и восстановим их вклад:

In [ ]:
# Находим 3 наиболее значимые частоты
top_n = 3
sorted_idx = sortperm(pos_ampl, rev=true)
top_freq = pos_freq[sorted_idx[1:top_n]]
top_ampl = pos_ampl[sorted_idx[1:top_n]]
top_phase = phase[sorted_idx[1:top_n]]

# Восстанавливаем гармоники
harmonics = zeros(n)
for (f, a, ϕ) in zip(top_freq, top_ampl, top_phase)
    harmonics .+= a .* cos.(2π * f * t .+ ϕ)
end

# Итоговый график
plot(t, detrended, label="Детрендированные данные", xlabel="Время (с)", ylabel="CO")
plot!(t, harmonics, label="Основные гармоники", linewidth=2)
Out[0]:

Полное разложение: тренд и гармоники.

In [ ]:
plot(t, y, label="Исходные данные", xlabel="Время (с)", ylabel="CO")
plot!(t, trend, label="Тренд", linewidth=2)
plot!(t, trend .+ harmonics, label="Тренд + гармоники", linewidth=2, linestyle=:dash)
Out[0]:

Вывод

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

  • Выделены долговременные колебания уровня CO, представленные явным трендом, определяемым движением транспортного средства и внешними факторами окружающей среды.
  • Определена структура периодических изменений концентрации CO, выраженная через выделение трёх главных частот колебаний.
  • Оценено качество реконструкции исходного сигнала, которое показало хорошую корреляцию между детрендинированными данными и основными гармониками.

Таким образом, разработанный подход оказался эффективным для анализа сигналов типа временно́й зависимости концентраций веществ, позволяя исследователям извлекать полезные сведения даже из зашумленных данных реального эксперимента.

Предложенный алгоритм легко адаптировать для анализа аналогичных временных рядов в иных областях исследований, где важно разделение данных на тренд и высокочастотные колебания.