Документация Engee
Notebook

Инженерные данные и сигналы

Импорт библиотек

Подключим библиотеки для работы с данными, сигналами и визуализацией. StatsBase и Statistics обеспечат статистические функции. CSV, XLSX и MAT позволят читать данные из соответствующих форматов. WAV - для аудиофайлов. Images - для загрузки изображений. DSP и FFTW - для спектрального анализа.

In [ ]: 
import Pkg
Pkg.add(["Statistics", "StatsBase", "StatsPlots", "CSV", "DataFrames", "Dates", "XLSX", "MAT", "WAV", "Images", "FileIO", "DSP", "FFTW", "DelimitedFiles"])
In [ ]: 
using Statistics, StatsBase, StatsPlots
using CSV, DataFrames, Dates, XLSX, MAT, WAV, Images, FileIO
using DSP, FFTW, DelimitedFiles
include("media_player.jl")
Out[0]:
base64encode_image (generic function with 1 method)

Чтение CSV-файла с данными энергопотребления

Начнём с самого распространённого формата инженерных данных - CSV. Прочитаем файл с показателями энергопотребления сталелитейного предприятия. Функция CSV.read возвращает DataFrame - табличную структуру, где столбцы доступны по именам. Сразу выведем размерность, чтобы понять масштаб данных.

In [ ]: 
df_energy = CSV.read("Steel_Industry.csv", DataFrame)
println("Размерность данных: ", size(df_energy))
println("Имена столбцов: ", names(df_energy))

Размерность данных: (35040, 11)
Имена столбцов: ["date", "Usage_kWh", "Lagging_Current_Reactive_Power_kVarh", "Leading_Current_Reactive_Power_kVarh", "CO2_tCO2_", "Lagging_Current_Power_Factor", "Leading_Current_Power_Factor", "NSM", "WeekStatus", "Day_of_week", "Load_Type"]

Первичный осмотр данных

Первое правило инженера при получении данных - посмотреть на них глазами. Выведем первые и последние строки, типы столбцов. Это поможет обнаружить проблемы: пропуски, некорректные типы, неожиданные значения.

In [ ]: 
println("Первые 10 строк:")
first(df_energy, 10)
In [ ]: 
println("\nТипы данных столбцов:")
[names(df_energy) eltype.(eachcol(df_energy))]

Типы данных столбцов:
Out[0]:
11×2 Matrix{Any}:
 "date"                                  String31
 "Usage_kWh"                             Float64
 "Lagging_Current_Reactive_Power_kVarh"  Float64
 "Leading_Current_Reactive_Power_kVarh"  Float64
 "CO2_tCO2_"                             Float64
 "Lagging_Current_Power_Factor"          Float64
 "Leading_Current_Power_Factor"          Float64
 "NSM"                                   Int64
 "WeekStatus"                            String7
 "Day_of_week"                           String15
 "Load_Type"                             String15

Преобразование даты и времени

В инженерных данных временные метки часто хранятся как строки. Преобразуем столбец date в формат DateTime. После преобразования можем вычислить диапазон дат и средний шаг измерений. Добавим новый столбец DateTime в DataFrame.

In [ ]: 
df_energy.DateTime = DateTime.(df_energy.date, dateformat"dd/mm/yyyy HH:MM")
println("Диапазон: ", minimum(df_energy.DateTime), " — ", maximum(df_energy.DateTime))

step_sec = mean([Dates.value(dt)/1000 for dt in diff(df_energy.DateTime)])
println("Средний шаг: ", round(step_sec, digits=2), " сек")

Визуализация временного ряда

График - лучший способ понять поведение системы. Построим график энергопотребления за весь период. По оси X - время, по Y - киловатт-часы. Видим цикличность и выбросы.

In [ ]: 
plot(df_energy.DateTime, df_energy.Usage_kWh,
     xlabel="Время", ylabel="кВт·ч",
     title="Энергопотребление",
     linewidth=0.5, legend=false, size=(900, 350))

Чтение данных из MAT-файла

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

In [ ]: 
подшипник = matread("Bearing.mat")
println("Ключи в файле: ", keys(подшипник))
println("Размер сигнала normal: ", size(подшипник["normal"]))

Преобразование MAT-данных в вектор

Данные из MAT-файла хранятся как матрицы размером 1×N. Преобразуем их в одномерные векторы с помощью функции vec. Теперь это обычный массив чисел Julia, с которым можно выполнять математические операции.

In [ ]: 
норма = vec(подшипник["normal"])
ролик = vec(подшипник["roller"])
println("Тип данных: ", typeof(норма))
println("Длина сигнала: ", length(норма))
println("\nПервые 5 отсчётов нормального сигнала: \n", норма[1:5])

Чтение аудиофайла WAV

Акустические данные часто хранятся в аудиоформатах. Функция wavread возвращает массив отсчётов и частоту дискретизации. В отличие от MAT, здесь данные сразу представлены как вектор.

In [ ]: 
сигнал_wav, fs_wav = WAV.wavread("Подшипник_1.wav")
println("Частота дискретизации: ", fs_wav, " Гц")
println("Длительность: ", round(length(сигнал_wav)/fs_wav, digits=2), " сек")

Воспроизведение аудиофайла в Engee

Воспроизведём акустический сигнал подшипника с помощью встроенного аудиоплеера.

In [ ]: 
media_player("$(@__DIR__)/Подшипник_1.wav", mode="audio")

Спектральный анализ через FFT

Быстрое преобразование Фурье переводит сигнал из временной области в частотную. Это ключевой инструмент анализа вибраций и акустики. Центрируем сигнал, применяем FFT, берём модуль для амплитудного спектра. Построим только первую половину (частоты от 0 до Найквиста).

In [ ]: 
сигнал = норма[1:4000] .- mean(норма[1:4000])
спектр = abs.(fft(сигнал))
N = length(сигнал)
частоты = (0:N÷2-1) * 10000 / N

plot(частоты, спектр[1:N÷2],
     xlabel="Частота, Гц", ylabel="Амплитуда",
     title="Акустический спектр подшипника",
     linewidth=1, legend=false, size=(800, 350))

Сжатие данных через FFT

Одна из практических задач - сжатие данных. Посмотрим на спектр и определим, в каком диапазоне сосредоточена основная энергия сигнала. Затем оставим только несколько ключевых частотных компонент - это и будет сжатое представление. Вместо 4000 отсчётов сигнала можно хранить 20 чисел.

In [ ]: 
амплитуды = спектр[1:N÷2]
энергия_общая = sum(амплитуды.^2)
кум_энергия = cumsum(амплитуды.^2) / энергия_общая

plot(частоты, кум_энергия,
     xlabel="Частота, Гц", ylabel="Накопленная энергия",
     title="Кумулятивная энергия спектра",
     linewidth=2, legend=false, size=(800, 350))
hline!([0.95], linestyle=:dash, label="95% энергии")

Выбор информативных частотных полос

Разобьём частотный диапазон на 20 равных полос и вычислим среднюю амплитуду в каждой. Это даст компактное представление сигнала - 20 чисел, которые можно хранить в таблице.

In [ ]: 
кол_полос = 20
полосы = range(0, 5000, length=кол_полос+1)
признаки_fft = zeros(кол_полос)

for i in 1:кол_полос
    маска = (частоты .>= полосы[i]) .& (частоты .< полосы[i+1])
    if sum(маска) > 0
        признаки_fft[i] = mean(амплитуды[маска])
    end
end

bar(1:кол_полос, признаки_fft,
    xlabel="Номер полосы", ylabel="Средняя амплитуда",
    title="Сжатое представление сигнала (20 полос)",
    legend=false, size=(700, 300))

Чтение изображения

Прочитаем файл с изображением. Отобразим тип изображения и его размер.

In [ ]: 
изобр = Images.load("image.png")
println("Тип изображения: ", typeof(изобр))
println("Размер: ", size(изобр))

Отобразим изображение.

In [ ]: 
изобр

Преобразование изображения в числовую матрицу

Преобразуем изображение в матрицу чисел.

In [ ]: 
каналы = Float64.(Images.channelview(изобр))
каналы = permutedims(каналы, (2, 3, 1))

Изображение представлено трёхмерным массивом 954×923×4. Первые два измерения это 954 строки и 923 столбца пикселей. Третье измерение — четыре цветовых канала: красный, зелёный, синий и альфа-канал прозрачности. Каждое значение находится в диапазоне от 0 до 1, где для цветовых каналов 0 означает отсутствие цвета, 1 — максимальную интенсивность, а для альфа-канала 1 означает полную непрозрачность.

Чтение текстового файла с данными

Текстовые файлы - простейший формат хранения. Прочитаем txt-файл с числовыми данными о температуре окружающей среды, определим количество значений, диапазон, и выведем первые 5 значений.

In [ ]: 
температура = vec(readdlm("температура.txt"))
println("Количество значений: ", length(температура))
println("Диапазон: ", minimum(температура), " ... ", maximum(температура), " °C")
println("Первые 5 значений: ", температура[1:5])

Чтение данных из Excel

Excel-файлы широко распространены в промышленности. Аналогично прочитаем xlsx-файл c данными о влажности окружающей среды.

In [ ]: 
xlsx_file = XLSX.readxlsx("влажность.xlsx")
лист = xlsx_file[1]
влажность = Float64.(лист["A"][2:end])
println("Количество значений: ", length(влажность),)
println("Диапазон: ", minimum(влажность), " ... ", maximum(влажность), " %")
println("Первые 5 значений: ", влажность[1:5])

Преобразование набора данных

Преобразуем набор данных, когда часть исходных столбцов не нужна, но требуется добавить новые данные из внешних источников. Удалим столбцы NSM, WeekStatus, Day_of_week и Load_Type, а затем добавим столбцы температуры и влажности, загруженные ранее.

In [ ]: 
select!(df_energy, Not([:NSM, :WeekStatus, :Day_of_week, :Load_Type]))
println("Столбцы после удаления: \n\n", names(df_energy))

Добавим ранее загруженные переменные температуры и влажности как новые столбцы. Количество строк должно совпадать.

In [ ]: 
df_energy.Temperature_C = температура
df_energy.Humidity_Pct = влажность

Отобразим первые 10 строк обновлённого набора данных.

In [ ]: 
println("Первые 10 строк:")
first(df_energy, 10)

Сохраним преобразованный набор данных в новый CSV-файл. Исходный файл Steel_Industry.csv остаётся без изменений - все преобразования сохранены в новом файле Сталелитейные_данные.csv.

In [ ]: 
CSV.write("Сталелитейные_данные.csv", df_energy)

Преобразование MAT-файла

Структура данных в MAT-файле также может быть представлена по-разному. Предположим, нам поступило 4 файла с данными подшипников с испытательного стенда:

  • H_1_0.mat – полностью исправный подшипник;
  • B_11_1.mat – неисправность роликового элемента;
  • I_1_1.mat – неисправность внутреннего кольца;
  • O_6_1.mat – неисправность внешнего кольца.

Но для модели анализа данных требуется иная структура – один MAT-файл с четырьмя ключами "normal", "roller", "inner", "outer", где:

normal – полностью исправный подшипник;

roller – неисправность роликового элемента;

inner – неисправность внутреннего кольца;

outer – неисправность внешнего кольца.

Загрузим первый файл.

In [ ]: 
норм_тест = matread("H_1_0.mat")

Данный файл имеет иную структуру данных, отличающуюся от требуемой. Отобразим данные.

In [ ]: 
норм_тест = норм_тест["H_1_0"]

Файл имеет один ключ "H_1_0, данные в котором является матрицей. Нам известно, что акустические данные расположены во втором столбце. Другие данные нам не нужны. Выполним извлечение этих данных.

In [ ]: 
норм_тест = норм_тест[:,2]

Откроем остальные MAT-файлы.

In [ ]: 
ролик_тест = matread("B_11_1.mat")
внут_тест = matread("I_1_1.mat")
внеш_тест = matread("O_6_1.mat")

Выполним извлечение матриц.

In [ ]: 
ролик_тест = ролик_тест["B_11_1"]
внут_тест = внут_тест["I_1_1"]
внеш_тест = внеш_тест["O_6_1"]

Далее аналогично выполним извлечение второго столбца из каждой матрицы.

In [ ]: 
ролик_тест = ролик_тест[:,2]
внут_тест = внут_тест[:,2]
внеш_тест = внеш_тест[:,2]

Теперь у нас есть четыре вектора с акустическими данными: норм_тест ролик_тест, внут_тест, внеш_тест. Их необходимо собрать в один MAT-файл обозначенной ранее структуры. Создадим словарь с нужными ключами.

In [ ]: 
данные = Dict(
    "normal" => норм_тест,
    "roller" => ролик_тест,
    "inner"  => внут_тест,
    "outer"  => внеш_тест)

Сохраним эти данные в новый MAT-файл.

In [ ]: 
matwrite("Подшипник_тест.mat", данные)

Прочитаем созданный файл.

In [ ]: 
matread("Подшипник_тест.mat")

Структура созданного файла соответствуем требуемой. Данные готовы для дальнейшего использования.