Linux. Типовая обработка научных данных

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

Чтобы ввести в курс дела я немного расскажу что и для чего делается. Зайду издалека.

Дело в том, что в мире насчитывается очень много производителей пожарных извещателей. Клепают извещатели все кому не лень! Понятно, что среди такого зоопарка сметчикам и проектировщикам очень сложно разобраться — какой извещатель лучше, а какой хуже. Надо как-то все извещатели приводить к общему знаменателю. Вот по этому-то и был придуман такой способ — создание искусственного пожара, параметры которого строго оговорены. Нет-нет, не надейтесь, дом поджигать никто не собирается!

В качестве эквивалента, создается типовой очаг — это такой металлический тазик, строго заданных размеров. В него наливается спирт (или другая горючая жидкость) в заданном количестве. Потом спирт в этом тазике поджигается, алкаши — рыдайте!

А пожарные извещатели, установленные, на опять-таки — на строго заданном расстоянии (все это определяется ГОСТ-ом), должны сработать за заданное время.

Представили картину?

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

Процесс оцифровки сигнала осуществлялся с частотой 1000 Гц на 10-разрядном АЦП (ATMEGA8). В результате получился файл с чуть более чем одним миллионом (ахренеть!) точек.

Когда горит пламя — свет о него как бы пульсирует, то есть во времени изменяется его яркость. Эти изменения очень незначительные, реально на 3-4 порядка меньше, чем средняя яркость. Поэтому, уловить эти колебания — задача не совсем тривиальная. Следует заметить, что частота колебаний лежит в диапазоне от одного до десятка-полутора Герц. Средина (горб) приходится на частоту 4-7 Гц, и не особенно ярко выражена. У разных веществ горения этот частный пик разный.

Полученные данные были записаны в файл в текстовом виде, каждое измерение (точка) записывалась с новой строки. Файл записывала младшая дочь, она и дала ему незатейливое название — spirt111. В общем так сложилось исторически.

Вот фрагмент этого файла:

94
0 489
0 485
0 482
0 482
0 482
0 483
0 484
0 484
0 485
0 486
0 487
0 490
...
0 501
0 504
0 504
0 503
0 500
0 495
0 489
0 484
0
0 503
0 501
0 498
0 493
0 487
0 483
...

В файле две колонки. В первой колонке — находятся данные по ультрафиолетовому каналу, во второй — по инфракрасному. Данные по обоим каналам синхронные. Данные по УФ каналу записывались так: ноль — соответствовал состоянию «нет УФ излучения», 1 — состоянию «получен импульс с УФ датчика». Данные по ИК каналу представляли значения, получаемые с АЦП микроконтроллера. Средний уровень, соответствующий отсутствию сигнала с ИК датчика был выведен примерно на средину шкалы. Вообще-то средина — это есть значение 512, но у меня получились несколько меньше. Это не проблема! В ИК канале нас не интересует постоянная составляющая. Там главное, чтобы сигнал был примерно симметричен относительно постоянного уровня и с сигналом можно было работать, то есть он не должен быть вылетать за диапазон АЦП и не должен быть слишком маленьким, чтобы не оказаться задавленным шумами квантования. Мне удалось добиться соблюдения этих требований не сразу, первые две записи оказались «провальными», но потом я научился и эта запись, которая описывается в этой статье, — почти идеальная.

Правда, по приезду домой мы все-таки обнаружили маленькую ошибочку. Файл начал записываться не с первой колонки, а почему-то со второй. Более того, примерно на 50-ой миллисекунде после начала записи (на 50-ои отсчете) оказалось пропущено значение из второй колонки. В прочем, это не составляет никаких проблем, так как файл текстовый и легко подправляется руками (Linux-way — рулит!). А потом, эти дефекты записи лежат вне времени, когда горел типовой очаг пожара.

В других файлах эта ошибка не появлялась, я счел, что она могла возникнуть из-за проблем связи стенда с компом по USB-каналу. Огневые испытания — это несколько нервозно-эмоциональная обстановка, и вполне могло случиться так, что сначала оказался запущенным стенд, а потом запустили комп на считывание с USB-порта. То есть, могло произойти элементарное переполнение буфера, в результате чего оказались потерянны некоторые данные. Но поскольку, эти данные являются «кадрами» до возникновения пожара, то никакой ценности они не представляют.

В результате «ручной коррекции» был создан файл spirt111.1. Правило исследовательских работ гласит — НИКОГДА НЕ ИЗМЕНЯЙ ПОЛУЧЕННЫХ ДАННЫХ! Если что-то нужно скорректировать, то создай копию и делай с ней, что хочешь!

В настоящей работе (в этой статье) меня не интересует ни синхронность по каналам, ни данные по УФ каналу. Поэтому мне нужно откинуть первый столбец. Для этого я воспользуюсь командой cut:

$ cut -d' ' -f2 spirt111.1 > spirt111.1-2

В результате у меня получился файл с одним столбцом:

489
485
482
482
482
483
484
484
485
486
487
490
494
498
502
504
504
502
499
...

Отлично! Теперь было бы не плохо посмотреть на полученные данные в виде графика. Воспользуемся программой MS Excel.

Я пошутил. Excel сразу умрет под тяжестью миллиона точек. Поэтому воспользуемся орудием профессионалов — да поможет нам gnuplot:

$ gnuplot

После запуска gnuplot в интерактивном режиме (и находясь в его среде) набираем незатейливую команду для построения графика:

gnuplot> plot "spirt111.1-2" w l

Наберите команду plot и (не нажимая пока кавычек) введите имя входного файла. Здесь, точно так же как и в shell, действует клавиша Tab, которая поможет вам быстро набрать имя файла. Затем обрамите кавычками имя файла. Это делать в общем-то необязательно, если в имени файла содержаться пробелы. Но привычка обрамлять имена файлов кавычками — это хорошая привычка.

Потом введите две буквы — ‘w’ и ‘l’ (маленькая латинская «эль»), — это сокращение «with line», что означает — рисовать график в виде линии. Количество типов графиков в gnuplot — очень огромное, но сейчас нам лучше отобразить наш график в виде линии.

Через продолжительное время (на моем стареньком компе это заняло примерно минуту) на экране появится прекраснейшая картина маслом — почти Васнецов или Репин:

spirt111.1-2

(Все графики кликабельны)

На графике хорошо видно, что ИК колебания начали увеличиваться примерно с 20-ой секунды, одна секунда — это 1000 точек (отсчетов). Это и есть начало возгорания. Очаг прогорел примерно до 960 секунды, что мы видим резкий спад колебаний.

Если вам режет слух использование слова «примерно», то gnuplot позволяет легко растянуть (масштабировать) график и получить точное значение. График «тянется» как по оси абсцисс, так и по оси ординат.

Есть еще два примечательных момента на графике. Начну со второго. На графике можно видеть достаточно большой выброс (иголку) примерно на 70 секунде (отсчет 70000). Ни одна живая душа не земле не знает — что это такое, но я по секрету скажу — это наш дорогой сотрудник (Евгений Яковлевич) решил щелкнуть зажигалкой рядом со стендом во время записи. К счастью, зажигалка не вспыхнула. А я на него зашипел… не важно. В общем, вот так, слегонца, Евгений Яковлевич увековечил себя на графике.

Первый же момент — оказался весьма неприятный. Посмотрите на на начало и конец графика — туда, где должен отсутствовать сигнал с ИК датчика. Вместо тонкой линии (отсутствия сигнала) на графике видим широкую «шлангу», размах амплитуды которой примерно 20-25 единиц.

А воспользуюсь масштабированием и увеличу маленький кусочек начального участка графика:

spirt111.1-2.50hz

Что мы видим? Мы видим четкий периодический, почти синусоидальный сигнал, почти неизменной амплитуды. Период сигнала — ровно 20 мс. Кто-нибудь уже догадался, что это такое?

Да, все правильно! Это ни что иное, как сетевая 50-герцовая наводка от сети. В этот раз на огневых испытаниях нам не очень повезло. Видимо, где-то рядом работала какая мощная энергоустановка, которая гадила помехами. До проведения этих испытаний, у меня запись сигнала проходила очень хорошо. Экранирование, защита стенд от помех работала и помех не было. Но в этот раз присутствовала не просто помеха, а мега-помеха. Это мы обнаружили уже дома.

Что делать? Не ехать же нам снова за двести километров на полигон, не проводить же заново дорогие испытания?

Первое, что приходит в голову — «прогнать» по полученным данным усредняющий фильтр с П-образной импульсной характеристикой. Ширина фильтра должна быть равна периоду помехи, то есть равна 20-ти отсчетам. Усреднив значения по близлежащим 20 отсчетам мы исключим полностью 50-герцовые колебания.

Для этого дела нужно обратиться к MathCad-у… Шучу! Давайте воспользуемся Питоном!

Я написал короткую программку, которая усредняет значения по 20 точкам:

#! /usr/bin/env python
# coding: utf-8

'''Произвести сглаживание входных данных путем усреднения.

В сигнале (входных данных) присутствует помеха сетевой частоты 50Гц.
Период помехи равен 20 мс. Каждый отсчет соответствует промежутку в 1 мс. Таким образом период помехи составляет 20 отсчетов.
'''

def average(buf):
  '''Вычислить средне-арифметическое значение массива buf
  '''
  sum = 0.0
  for val in buf:
    sum += val
  result = sum / len(buf)
  return result

def smooth20(src):
  '''"Пробежаться" по исходному массиву данных src и усреднить
  значения соседних 20 отсчетов.
  '''
  buflength = 20  # количество отсчетов для усреднения
  res = []        # выходной массив данных, 50ти-герцовая помеха подавлена
  for i in xrange(len(src) - buflength):
    buf = src[i:(i + buflength)]
    res.append(average(buf))
  return res

#=============================================================================
file_in = open("spirt111.1-2", "r")
strIR = file_in.readlines()  # Строковый массив входных данных
file_in.close()

ir = []  # Целочисленный массив входных данных
for s in strIR:
  ir.append(int(s))

smoothed = smooth20(ir)  # Сглаженный массив

for i in smoothed:
  print i

Теперь запустим эту прогу, а ее выход направим в файл ir-stat.smoothed:

$ ./s2.py > spirt111.1-2.smoothed

Десяток-другой секунд и мы имеем отфильтрованный сигнал. Давайте снова запустим gnuplot и посмотрим на график. К стати, вы в курсе, что gnuplot тоже хранит историю команд, которые были набраны в интерактивном режиме? Эту историю можно пролистывать с помощью клавиш управления курсором — «вверх», «вниз».

Прокрутите историю вверх и измените имя файла с исходными данными:

gnuplot> plot "spirt111.1-2.smoothed" w l

В результате появится «чистенький» график:

spirt111.1-2.smoothed

Мы видим, что фильтр не только подавил 50-герцовые колебания, то он еще и задавил все высокочастотные — те, которые находятся выше 50 Гц. Я снимал данные с частотой 1000 кГц, поэтому, в данных могли присутствовать все частоты аж до 500 Гц. Частоты ниже 50 Гц не затронуты. Зато мы видим — мелкое хулиганство Евгения Яковлевича тоже подавлено.

Хорошо, но это еще не всё! Публикуя эту статью, я ставил задачу, показать людям способ работы с экспериментальными данными в Линуксе. Мне в последнее время очень часто говорят, что подобные работы лучше производить не из-под Линукса, в в Виндовсе. Честно говоря, я не понимаю, чем это желание продиктовано. Скорее всего только тем, что люди не владеют знаниями Линукса и советуют только то, что сумели как-то освоить.

В Линуксе всё получается намного лучше, чем под Виндой. А самое главное — все процессы ты сам контролируешь. Ты понимаешь — что ты делаешь, и как оно делается.

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

Честно говоря, этот график имеет не много практического значения в свете разработки пожарных извещателей. Мне просто было интересно обработать довольно-таки мощный массив данных и посмотреть на его статистику.

График плотности вероятности показывает с какой вероятностью… или как часто возникает (встречается) то или иное значение уровня сигнала. Вот так, не глядючи на график плотности вероятности, на основе графика исходных данных можно сказать, что наиболее часто будут возникать (выпадать) значения в области 490-495 единиц. Значения ближе к нулю и к уровню 900 будут практически отсутствовать. То есть их вероятность возникновения будет очень мала.

Для вычисления плотности вероятности и написал еще одну простую программку на Питоне:

#! /usr/bin/env python
# coding: utf-8

'''Вычисление плотности вероятности на основе файл с данными

Исходные данные находятся в текстовом файле в виде столбца чисел
Плотность вероятности -- это массив шириной 1000 значений
'''

#=============================================================================

file_in = open("spirt111.1-2.smoothed.200-500", "r")
strIR = file_in.readlines() # Данные в виде строк текста
file_in.close()

# Конвертировать текстовые данные в числа
ir = []
for s in strIR:
  ir.append(int(float(s)))

# Создать массив вероятности
a = []
for i in xrange(1000):
  a.append(0)

# Подсчитать вероятность значений
for j in ir:
a[j] += 1

# Вывести массив вероятности
for i in a:
  print i

Программа забирает значения из входного файла, имя которого прописано в самой программе, и выдает на стандартный выход значения созданного массива плотности вероятности. Ширина этого массива равна ширине (или размаху) входных значений исходного сигнала.

Получить на консоль 1000 чисел — в этом нет ничего толкового. Давайте сохраним их в файле:

$ ./s1.py > ir-stat.smoothed

а полученный файл посмотрим в виде графика с помощью gnuplot. Команды я не буду приводить, я думаю, что вы и сами уже можете их составить. Вот что у меня получилось:

ir-stat.smoothed

Интересно, не правда ли! Теперь, если тоже самое проделать с исходными данными, которые не были подвергнуты 50-герцовому фильтру, то можно увидеть вместо одного серединного пика — два.

ir-stat

Ну, хорошо. На временнОм графике мы видим, что у нас есть несколько областей ИК сигнала:

— две области в начале и в конце записи, когда полезный сигнал от очага пожара отсутствовал.
— область разгорания очага пожара (примерно 1-2 минуты)
— область стабильного горения (если быть честным, то нифига оно не стабильное! — Оно ОТНОСИТЕЛЬНО стабильное. Мы видим, что процесс горения медленно увеличивался до момента, пока не выгорел весь спирт. При горении бензина — картина еще поганее. Только кто об этом из сертификаторов в курсе!)
— область затухания, когда пламя пошло на спад (минута или около того)

Все эти области попали в график плотности вероятности, и там оставили свой след. Так например, области, когда отсутствовал сигнал, дали на графике плотности вероятности узкие и одновременно высокие выбросы. У нефильтрованного сигнала — это два «уха», а у фильтрованного — один высокий пик.

А что если мы отсечем все области «нестабильного» горения? Давайте вырежем участок сигнала с 200-ой по 500-ю секунды.

Для того чтобы откинуть начало и конец данных воспользуемся командами head и tail. Команда head возьмет из входного файла первые 500000 строк и выдаст их на стандартный выход. (Остаток строк файла будет проигнорирован этой командой.) Мы перехватим этот выход, и «скормим» его команде tail. Команда tail проигнорирует первые строки из входного потока, а на свой выход выдаст последние (считая от конца) 300000 строк. Вы помните, что одной секунде сигнала соответствует 1000 строк?

Усеченный таким образом поток данных мы запишем в файл spirt111.1-2.smoothed.200-500:

head -500000 spirt111.1-2.smoothed | tail -300000 > spirt111.1-2.smoothed.200-500

Вот как выглядит средний участок потока:

spirt111.1-2.smoothed.200-500

Видно, что амплитуда колебаний сигнала на всем протяжении времени измерения не сильно меняется. Я вот сейчас смотрю на этот график и думаю — похоже, этот сигнал  можно даже помещать в учебники как пример квазистационарного шумового сигнала.

Теперь, с помощью программы s1.py получим плотность вероятности и запишем её в файл ir-stat.smoothed.200-500. Далее, с помощью gnuplot посмотрим, что у нас получилось:

ir-stat.smoothed.200-500

Ну что ж, не плохая работа! Почти Гауссовская колоколо-образная кривая!

А теперь было бы интересно еще посмотреть на все три графика плотности вероятности, но совмещенные на одной плоскости. Для этого в gnuplot введем следующую команду:

gnuplot> plot "ir-stat" w l, plot "ir-stat.smoothed" w l, plot "ir-stat.smoothed.200-500" w l

ir-stat-3

Конечно, если покопаться, то можно сделать эту работу еще более красиво. Но дело не в том — оптимально или не оптимально написаны программы на Питоне, и на сколько коряво написаны. Моя задача была в том, чтобы вообще поработать с данными и получить результаты. Вторая задача — поделиться полученным опытом с вами.

Эти программы не коммерческие. Они написаны далеко не идеально. Их можно еще улучшать и совершенствовать. Исходный текст открыт и не требует лицензирования. Совершенно нет никакой проблемы в том, что бы изменить имя входного файла и снова запустить программу. Мы не пижоны, мы исследователи. Свои инструменты мы создаем и подправляем по мере необходимости. Наша задача получить результат обработки научных данных, а не написать универсальную прогу на все случаи жизни. Почувствуйте разницу в подходе к созданию программы — цели разные. Абсолютно разные!

Я буду считать, что я не зря потратил свое время, если вы чему-то научились благодаря этой статье.

Реклама

Добавить комментарий

Заполните поля или щелкните по значку, чтобы оставить свой комментарий:

Логотип WordPress.com

Для комментария используется ваша учётная запись WordPress.com. Выход / Изменить )

Фотография Twitter

Для комментария используется ваша учётная запись Twitter. Выход / Изменить )

Фотография Facebook

Для комментария используется ваша учётная запись Facebook. Выход / Изменить )

Google+ photo

Для комментария используется ваша учётная запись Google+. Выход / Изменить )

Connecting to %s