НАЦИОНАЛЬНЫЙ ИССЛЕДОВАТЕЛЬСКИЙ ЯДЕРНЫЙ УНИВЕРСИТЕТ «МИФИ»
Кафедра информатики и процессов управления (№17)
Дисциплина «Прикладное ПО»
для студентов факультета КиБ,
3-й курс, 5-й семестр.
Методические указания
Занятие 2
Статистическая обработка данных и анализ сигналов.
Содержание
Статистическая обработка данных и анализ сигналов. 1
Индивидуальные варианты.. 2
Исходные данные. 2
Загрузка матриц в Matlab (Octave) 4
Считывание из файла. 4
Извлечение векторов из статистического массива-матрицы.. 5
Логическая индексация. 6
Булевы выражения. 6
Фильтрация некорректных данных при помощи логической индексации. 6
Задание для самостоятельного решения №1. 7
Условие. 7
Оценка. 7
Пояснение. 7
Дополнительно. Анализ сигналов. 8
Индивидуальные варианты
№ | Город | № | Город | № | Город | ||
1 | Архангельск | 11 | Калининград | 21 | Магадан | ||
2 | Астрахань | 12 | Калуга | 22 | Мурманск | ||
3 | Белгород | 13 | Кемерово | 23 | Новгород | ||
4 | Брянск | 14 | Киров (Вятка) | 24 | Новосибирск | ||
5 | Владимир | 15 | Кострома | 25 | Омск | ||
6 | Волгоград | 16 | Краснодар | 26 | Оренбург | ||
7 | Вологда | 17 | Красноярск | 27 | Орел | ||
8 | Воронеж | 18 | Курган | ||||
9 | Иваново | 19 | Курск | ||||
10 | Иркутск | 20 | Липецк |
Исходные данные
В качестве примера будут использоваться данные, полученных из архива наблюдений за погодой.
Зайдите на сайт http://rp5.ru/ и скачайте архив погоды за январь 2015 года (а если данных нет – то на другой месяц) на МЕТЕОСТАНЦИИ в соответствии с индивидуальным номером своего варианта (номер студента в списке группы). !! Не берите архив аэропорта, там данные иной структуры.
Из полученного архива *.gz извлеките файл XLS и откройте его в Excel для некоторой предварительной обработки. Обратите внимание, что файл имеет сложную структуру, поэтому его нельзя непосредственно использовать для импорта данных в MATLAB при помощи простейшей команды load. (Можно использовать более сложные инструменты MATLAB, но лучше выполнить предварительную обработку):
Данные в таком виде не могут быть восприняты командой «load» без предварительной обработки, поскольку Matlab (Octave) спроектированы в основном для работы с числовыми матрицами. Основными препятствиями для загрузки являются: наличие текстовых строчек, наличие заголовков столбцов, имеются пустые ячейки (без данных), присутствуют нечисловые параметры (форма облаков, структура снега,…), Дата+Время присутствуют в таблице в виде текста. Предупреждение: не следует «входить» внутрь ячейки с датой, поскольку Excel после этого изменит содержание ячейки.
Выполним следующие действия:
· Удалим первые 7 строчек (включая заголовки столбцов), оставим только непосредственно табличные данные.
· Удалим все столбцы, кроме первых шести, поскольку уже в седьмом столбце присутствуют нечисловые данные.
· И еще необходимо заполнить пустые клетки значением NaN, чтобы Matlab (Octave) понимал, что это НеЧисло.
Выделите прямоугольную область с данными и используйте команду «Заменить все», отметив «Ячейки целиком»:
Результат сохраните в виде текстового файла с разделителями табуляция под именем файла meteo.dat, поместите его в рабочую директорию Matlab (Octave).
Загрузка матриц в Matlab (Octave)
Считывание из файла
Команда load считывает текстовые файлы, содержащие численные данные. Текстовые файлы должны быть сформированы в виде прямоугольной таблицы чисел, отделенных пробелами, с равным количеством элементов в каждой строке. Команда:
> load meteo.dat
прочитает этот файл и создаст переменную meteo, содержащую эту матрицу с архивом погоды. Посмотрите в окне просмотра переменных содержание этого массива:
31,01 21 -3,500 752,6 753,4 -0,5000 100
31,01 18 -4,200 753,1 753,9 0,0 100
31,01 15 -3,500 753,1 753,8 0,6000 96
31,01 12 -2,100 752,5 753,3 1,400 89
31,01 9 -2,000 751,1 751,9 1,0 92
В третьем столбце – температура. В четвертом и пятом – давление на разных уровнях, потом точка росы и влажность. Во втором столбце – время, в первом столбце – дата, воспринятая при импорте, как вещественное число.
Чтобы привести дату к более разумному виду, применим операцию округления к первому столбцу:
> meteo(:,1)=round( meteo(:,1));
Посмотрите в окне просмотра переменных, как изменился этот массив:
31 21 -3,500 752,6 753,4 -0,5000 100
31 18 -4,200 753,1 753,9 0,0 100
31 15 -3,500 753,1 753,8 0,6000 96
31 12 -2,100 752,5 753,3 1,400 89
31 9 -2,000 751,1 751,9 1,0 92
И сделаем еще одну операцию (сортировку строк), чтобы данные стали иметь общепринятый вид. Сортировка строк в матрице осуществляется командой sortrows(A, [m,n,k]) , где A – матрица, m n k и далее – номера столбцов, которые используются при сортировке, то есть сортируем по значению в столбце m, а в случае равенства значений сортируем по значению в столбце n, и так далее. Если [m,n,…] не указаны, то подразумевается, что используется [1,2,…], то есть сортируем по первому столбцу (в нашем случае – дата), а при равенстве дат – сортируем по второму столбцу (время). Для сортировки данных по дате и времени введите команду:
> meteo=sortrows(meteo);
Посмотрите в окне просмотра переменных, как изменился этот массив:
1 0 -13,40 749,0 749,8 0,9000 93
1 3 -12,40 749,2 750.0 0,2000 93
1 6 -11,10 749,4 750,1 0,2000 94
1 9 -9,800 748,9 749,6 -0,5000 95
1 12 -9,200 748,8 749,5 -0,1000 96
Теперь данные отсортированы. Можно увидеть данные графически, например, температуру:
> plot(meteo(:,1),meteo(:, 3),'+')
На графике видно, что каждый день было снято несколько измерений. Число измерений за день может быть не везде одинаково.
Извлечение векторов из статистического массива-матрицы
Построим вектор из второго столбца матрицы (это время наблюдения time, в часах), и из третьего столбца, где записана температура (temp):
> time=meteo( : , 2 );
> temp=meteo( : , 3);
Логическая индексация
Булевы выражения
В Matlab (Octave) можно использовать логические (булевы) выражения, имеющие значения: истина (1) или ложь (0). Попробуйте два простых примера:
> 5 > 3
> 5 >= 10
Вектора и матрицы также могут иметь булевы элементы. Ниже пример вектора, содержащего булевы значения. В данном случае анализируется, чтобы синус от индекса был больше косинуса:
> sin(1:10) > cos(1:10)
Если вектор из булевых значений используется в качестве индекса, то это называется логической индексацией. Она означает, что из массива выбираются только те элементы, для которых индекс имеет значение Истина.
В качества примера построим вектор t из значений температуры в полдень (12 часов). time == 12 – это условие равенства, результатом проверки условия будет булевый вектор, который используется в качестве индекса
> t = temp( time == 12 );
> t'
Построим график температуры в полдень по дням:
> plot(t)
Фильтрация некорректных данных при помощи логической индексации
Предположим, присутствуют неверные данные, по причине отказа датчика, или ошибок при передаче информации. Смоделируем такую ситуацию, преднамеренно изменим температуру 15 января и сотрем значения за 20 января:
> t(15) = 100 ;
> t(20) = NaN ;
При помощи логической индексации проверяем, все ли данные присутствуют, и отбрасываем пустые значения:
> t = t( isfinite(t) );
Функция isfinite() является истиной для всех конкретных значений и ложью для NaN и Inf. (В некоторых реализациях она называется finite() ).
При помощи логической индексации и встроенных статистических функций
mean() – среднее,
std() – среднеквадратичное отклонение
можно отфильтровать результаты и отбросить заведомо некорректное измерение за 15 января:
> t = t( abs(t-mean(t)) <= 3*std(t) );
Теперь данные можно считать корректными и можно рассчитать среднюю дневную температуру за январь, дисперсию и другие важнейшие статистические значения:
> mean(t)
> std(t)
> max(t)
> min(t)
> median(t)
Задание для самостоятельного решения №1
Условие
Рассчитать среднесуточную температуру по дням и построить ее график.
Оценка
Максимальный балл: 6.
Пояснение
Используя изученные приемы работы можно рассчитать среднюю суточную температуру. Например, за первое число: выбираем все записи, у которых дата совпадает с датой в первой строке, и усредняем:
> t1=meteo((meteo(1,1)==meteo(:,1)),3)
> mean(t1)
Для расчета среднесуточной температуры в каждый из дней воспользуемся циклом. Предварительно запомним величину сдвига в поле дат, чтобы не обращаться всякий раз к ячейке meteo(1,1):
> base=meteo(1,1) - 1
> for i=1:31 Sr(i)=mean(meteo(base+i==meteo(:,1), 3)); end
Можно построить график среднесуточной температуры:
> plot(Sr)
Дополнительно. Анализ сигналов
Рассмотрим композицию из трех синусоидальных сигналов разной частоты:
> N=1024 ;
> x= 1:N ;
> y1=sin(x);
> y2=sin(x/3);
> y3=sin(x/8);
Сложим эти сигналы, предположим, мы слушаем музыку, где одновременно звучат эти три ноты (гармоники) :
> y = y1 + y2 + y3 ;
Matlab (Octave) имеет встроенную функцию быстрого преобразования Фурье (fast Fourier transform) fft(), которое позволяет провести спектральный анализ полученного сигнала (звука) :
> S = fft(y);
> freq = (1:N/2)/N;
> power = abs(S(1:N/2));
> plot(freq, power), grid on
> xlabel('Frequency')
> title('Spectral analysis')
На графике явно видно три пика, соответствующих трем гармоникам. Например, значение частоты ω=0.16 примерно соответствует первой из них: sin(2π ω x) ≈ sin(x).
Предположим, музыка нам понравилась, и мы увеличили громкость в наушниках (или в усилителе), например, в 4 раза
> y = 4*y;
Повторим исследование, посмотрим, как изменился спектральный состав нашего звука:
> S = fft(y);
> power = abs(S(1:N/2));
> plot(freq, power)
На наше счастье, ничего не изменилось.
Но это только потому, что наушники были хорошие, был запас по мощности усиления.
А теперь посмотрим, что произойдет, если увеличение громкости произойдет нелинейно. Пусть физические параметры усиления ограничены значением: 4. Обрежем наш сигнал по планке:4:
> y = min(y, 4);
Повторим исследование, посмотрим, как изменился спектральный состав нашего звука:
> S = fft(y);
> power = abs( S(1:N/2));
> plot(freq , power)
На графике должны отразиться посторонние шумы.
