Размещено на http://www.allbest.ru/

ТРИГОНОМЕТРИЧЕСКАЯ АППРОКСИМАЦИЯ

ДИСКРЕТНЫХ ПОСЛЕДОВАТЕЛЬНОСТЕЙ

Вычисление естественных спектров

М.Я. Хандрос, к.т.н.

Оглавление

Введение

1. Типы тригонометрических полиномов в пакете

2. Числовые характеристики данных и полинома

3. Особенности тригонометрической аппроксимации

4. Алгоритм аппроксимации

5. Проверка функционирования алгоритма

6. Предварительная обработка данных

7. Модификация алгоритма гармонического анализа

8. Практические примеры расчёта

8.1 Тяговое сопротивление плуга

8.2 Крутящий момент

8.3 Анализ солнечных пятен

Заключение

Литература

Введение

Дискретные последовательности могут иметь различную физическую природу. Однако, с математической точки зрения - это набор чисел, упорядоченный по номерам ,i=1..n. Номер соответствует конкретному значению независимой переменной t=(i-1)h, где h - интервал фиксации, n - количество элементов последовательности. Независимая переменная может быть временем, длиной, углом поворота и пр. Интервал фиксации измеряется в годах, днях, секундах, метрах, сантиметрах, градусах, радианах и пр.

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

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

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

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

, (1)

тригонометрический полином аппроксимация

где параметры - неотрицательные действительные числа, t-независимая переменная (аргумент) экспериментальной последовательности. Параметр называется частотой и имеет размерность, обратную аргументу: рад/с, рад/м, рад/день и пр. Для аргумента - угла в радианах параметр безразмерный. Относительно элементов системы (1) тригонометрические полиномы являются обычными алгебраическими полиномами.

В настоящее время основным инструментом для приближения данных

тригонометрическими полиномами является аппарат рядов Фурье. При этом используются конечные отрезки таких рядов, построенные на ограниченном интервале изменения аргумента [0,T]. Особенностью полиномов Фурье является то, что набор частот для них не определяется по исходным данным, а задаётся заранее в зависимости от длины интервала T. Все частоты выбираются кратными основной частоте. Такой подход позволяет для вычисления коэффициентов полинома использовать весьма простые, удобные формулы. Однако, если искусственно назначенных частот в экспериментальном процессе нет, то ряд Фурье сходится очень медленно и подменяет каждую естественную частоту несколькими искусственными.

При таких обстоятельствах для аппроксимации целесообразно использовать полиномы, для которых вычисляются не только коэффициенты при косинусах и синусах, но так же их естественные частоты [1], которые в общем случае несоизмеримы и могут быть любыми неотрицательными действительными числами. Определение естественных частот представляет значительные трудности, так как система уравнений для их определения громоздка, многоэкстремальна и относится к классу задач нелинейного программирования. Из многих локальных экстремумов нужно находить глобальный минимум.

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

Программы для определения параметров тригонометрических полиномов, как и данный документ, находятся в отдельной папке TRIGA. Они написаны на языке компьютерной математики Maple и объединены в пакет. Разработка и отладка программ произведена в релизах 11 и 12 для работы в классическом интерфейсе, но никаких специфичных команд и пакетов именно для этих версий Maple не использовалось. Поэтому программы могут работать и в младших версиях системы Maple. Однако, при их использовании в стандартном интерфейсе возникают трудности с интерпретацией кириллицы в текстовых комментариях.

Пакет программ TRIGA включает пять программ:

vorarb.mws- предварительная статистическая обработка данных;

furpol.mws- классический гармонический анализ с помощью ряда Фурье;

linpol.mws- определение параметров линейного по косинусам и синусам полинома;

multiqvpol.mws- определение параметров квадратичного тригонометрического полинома;

cubpol.mws- определение коэффициентов и частот неполного кубического тригонометрического полинома

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

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

1. Типы тригонометрических полиномов в пакете

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

рассматриваются практические вопросы аппроксимации тригонометрическими полиномами.

Экспериментальные последовательности,которые подлежат аппроксимации, c математической точки зрения являются ограниченными на всём интервале наблюдений. Поэтому непрерывные процессы , которые они представляют, так же можно считать равномерно ограниченными на том же интервале. Они могут иметь чётко выраженный сезонный характер, например, колебания температуры воздуха в определённом географическом пункте имеют в течение года всем известные четыре сезона. В следующем году сезоны повторяются. Однако температуры по дням внутри сезонов не совпадают. Для температурного ряда нельзя подобрать такого значения , чтобы на всей длине записи соблюдалось условие

,t=0..T.

Однако, это условие может соблюдаться приближённо

.

Здесь -- положительное число.

Функции, удовлетворяющие последнему условию, называются почти периодическими [3]. Числа называются почти периодами. Если величина невелика по сравнению с наибольшим модулем элемента, то запись такого процесса похожа на полигармонический процесс. Сумма гармонических составляющих с несоизмеримыми частотами представляет собой почти периодическую функцию [4,5]. Согласно теореме об аппроксимации, для почти периодических функций можно подобрать c требуемой точностью конечный тригонометрический полином [6]

,

где могут быть любыми положительными, отличными друг от друга действительными числами.

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

В качестве порождающей функции для тригонометрических полиномов пакета служит алгебраический полином многих переменных

, (2)

где m - целое чётное, трёх видов --постоянные коэффициенты,- переменные.

Количество одночленов в полиноме даётся формулой

. (3)

Для перехода к тригонометрическому полиному необходимо следующим образом обозначить переменные в выражении (2)

, ,

, ,

.. .. .. .. (4)

, .

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

,..,

Когда аппроксимирующий полином является квадратичной функцией, он имеет такой вид

+..++..

++..++.. (5)

++..++..

++..+

… … … … … …

+.

Здесь

.

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

.

В простейшем полиноме число слагаемых равно шести, а в полиноме с десятью частотами - шестидесяти шести.

Учитывая громоздкость такого полинома и затраты времени на расчёты, реальное количество частот, которое можно определить с помощью такого полинома не превышает десяти.

Важнейшим частным случаем тригонометрического полинома (2-4) является линейный полином

. (6)

При рассмотрении отдельной составляющей полинома (6) индекс j можно опускать. Тогда

. (7)

Для удобства сопоставления с отрезком ряда Фурье в выражении (6) коэффициенты при косинусах можно обозначить через a , при синусах через b. В таком случае

. (8)

Если в уравнении (8) положить , где T=(n-1)h -длина полного интервала фиксации экспериментальной последовательности, получаем полином Фурье для этого интервала.

Координатные функции ортогональны на отрезке [0,T], что существенно упрощает вычисление коэффициентов полинома. Однако при этом необходимо помнить, что такие частоты в процессе x(t) могут в принципе отсутствовать. В этом случае для аппроксимации с заданной точностью может потребоваться большее число гармоник, чем для случая, когда частоты достаточно точно определяются по экспериментальным данным.

Если соотношения частот /,k=1..m -- рациональные числа, то полином представляет собой периодическую функцию, если при этом отношения /, j=1..m равны j, то полином - классический ряд Фурье. Если же соотношения частот - действительные числа общего вида, то полином -- почти - периодическая функция.

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

. (10)

На рис.1. изображены две составляющие этого полинома. Эти составляющие имеют более гибкую и разнообразную форму, чем гармоники линейного полинома.

Рис.1. Вид составляющих неполного кубического полинома с разными частотами и коэффициентами.

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

В заключение данного параграфа ещё раз остановимся на свойствах почти периодических функций. Все рассмотренные здесь тригонометрические полиномы относятся к классу почти периодических функций, так как сумма двух и более периодических функций с действительными частотами является почти периодической функцией. Действительные частоты, найденные по экспериментальным данным, мы называем естественными[1]. Полином (8) для случая, когда его частоты -- целые кратные одного и того же числа, кроме того является периодической функцией и одновременно рядом Фурье.

2. Числовые характеристики данных и полинома

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

.Средние значения:

, (11)

где - экспериментальные данные, - соответствующие значения полинома. Суммы квадратов отклонений:

. (12)

В дальнейшем для краткости эти суммы будем называть вариациями. Остаточная сумма (вариация) используется для оценки точности аппроксимации.

Стандартные отклонения:

, (13)

гдеколичество ограничений. Для и величина . Для величина

равна количеству параметров полинома. Стандартные отклонения являются

мерой разброса конкретных данных.

Центральные моменты:

. (14)

Как и предшествующие числовые характеристики, моменты вычисляются и для экспериментальных данных, и для вычисленных дискретных значений полинома.

Асимметрия эксцесс:

. (15)

Коэффициенты и используются для приближённого выбора эмпирической функции распределения и выдвижения гипотез [7,8]. Для каждой такой функции определяется зависимость и в координатах (,) строятся специальные диаграммы.

На этих диаграммах различным типам распределений соответствуют точки, линии и области между линиями. Например, нормальному распределению соответствует точка с координатами (0,3),экспоненциальному --(4,9), Релея--(0.398,0.245). -распределения связаны формулой

,

и этому распределению на диаграмме соответствует прямая линия.

Размах

. (16)

Размах служит ещё одной мерой разброса для данных и полинома.

Взаимная корреляционная функция[9].

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

. (17)

Тогда дискретная корреляционная функция

. (18)

является коэффициентом корреляции между полиномом и данными.

ТочностьРазмещено на http://www.allbest.ru/

аппроксимации. Точность аппроксимации оценивается по относительной квадратичной ошибке и корреляционной функции:

. (19)

Допустимая ошибка задаётся. Как правило, достаточно точности 0.05. Уменьшение допустимой точности приводит к увеличению числа слагаемых полинома, последние из них могут быть чисто шумовыми.

Корреляционная функция определяет не только разброс взаимных отклонений, но и взаимное поведение полинома и последовательности данных. Если коэффициент корреляции близок к 1, то всюду, где последовательность возрастает , растёт и полином, а там, где значения последовательности уменьшаются, убывает и полином. Иными словами при близком к 1 значении коэффициента интервалы возрастания и убывания полинома и дискретной последовательности почти совпадают.

Числовые характеристики в пакете TRIGA вычисляются не только для полинома в целом, но и для отдельных составляющих и коэффициентов, если их рассматривать в качестве случайных величин.

Особенности тригонометрической аппроксимации

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

(20)

В данной формуле -- граничная частота Котельникова - Найквиста [4,10,11], K-поправочный коэффициент, который определяется свойствами исследуемого процесса и объекта. Например, динамические системы различной природы всегда обладают определённой инерционностью и на частоты выше определённой пороговой частоты практически не реагируют. Поэтому при расчётах можно рекомендовать величину K брать в пределах от 1 до 4.

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

, (21)



где и --составляющие полинома и их частоты, - вектор значений коэффициентов j-й составляющей. Неизвестные входят в выражение (21) внутри функций синуса и косинуса.

Критерий (21) представляет собой относительную квадратичную погрешность аппроксимации. Чем меньше величина , тем лучше приближение. Относительная форма критерия удобна тем, что для неё легко выбирать допустимое значение погрешности . Весьма часто величина принимается равной 0.05, т.е. при расчётах требуется, чтобы 0.05.

Необходимо обратить внимание, что при определении составляющих полинома значения , уже известны. Поэтому величина является функцией наборов переменных. Величину m считаем заданной, хотя на первом этапе вычислений она тоже неизвестна. Нужно найти такие значения переменных, при которых достигается минимум критерия (21). В регрессионном анализе [12] поступают следующим образом. Дифференцируют равенство (21) по каждой неизвестной и полученные выражения приравнивают нулю. В результате получается система нормальных уравнений:

=0

=0 (22)

В нелинейном случае для решения системы задаются начальными значениями переменных и применяют различные итеративные численные методы[12,13]. Успешное решение задачи зависит от сложности уравнений, особенностей используемого итеративного метода и удачного выбора начальных значений. В тригонометрической аппроксимации эти приёмы не подходят из-за специфических особенностей функции (21) и системы (22).

Система уравнений очень громоздка. При десяти составляющих для определения линейного полинома она содержит 40 неизвестных, а для определения кубического полинома --70. Система существенно нелинейна и не поддаётся линеаризации. Функция (21) многоэкстремальна, что делает невозможным выбор подходящей начальной точки. Даже, если решение задачи с конкретным m будет найдено, ещё необходимо проверить условие 0.05. Если условие не выполнено, то необходимо значение m увеличить и произвести пересчёт параметров полинома.

4. Алгоритм аппроксимации

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

. (23)

До начала основных вычислений необходимо задать допустимую погрешность ,

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

Специфика тригонометрического полинома такова, что удаётся разделить определение частот и коэффициентов на три последовательных этапа: выбор частоты, расчёт коэффициентов полинома, проверка значения критерия. Выбор частоты осуществляется в соответствии с методом статистических испытаний

[14,15], в интервале генерируется некоторое количество равномерно распределённых случайных чисел. Из этого набора последовательно выбираются частоты. Для каждой частоты обычным линейным методом наименьших квадратов МНК [11,16] вычисляются коэффициенты полинома.

В матричной форме это делается следующим образом. Вычисляется матрица

(24)

Матрица V имеет размерность [n,3] и составлена для линейного полинома типа (7), для кубического полинома соответствующая матрица имеет 6 столбцов и содержит строки вида

,i=1..n . (25)

Синусы и косинусы в матрице с целью экономии машинного времени вычисляются по рекуррентным формулам:

, i=1,..,(n-1),

, (26)

.

В матричных вычислениях экспериментальная дискретная последовательность рассматривается как вектор - столбец x с n компонентами и для хранения данных вводится дополнительный столбец X=x. Коэффициенты определяются c помощью матричного уравнения

X=VC , (27)

где C--вектор-столбец искомых коэффициентов. Для определения коэффициентов методом наименьших квадратов в Maple предусмотрена специальная функция. Найденный таким образом вектор коэффициентов C обеспечивает выход на минимум критерия (23) при определённой частоте .

Затем по матричным формулам вычисляются вариации для отдельной составляющей:

, (28)

где Т--знак операции транспонирования матрицы

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

быстрее, чем просто по формулам (12). Если все найденные параметры запоминаются во вспомогательных переменных:

. (29)

После этого выполняется уточняющий этап данной стадии вычислений. В диапазоне

(30)

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

Значения параметров составляющей, которые соответствуют наименьшему значению критерия, сохраняются в матрице результатов B. Ниже приведены элементы j-й строки матрицы:

(31)

Параметры каждой составляющей хранятся в отдельной строке. В пятый столбец помещаются амплитуды установленных составляющих.

Затем проверяется условие и, если оно выполняется, то вычисления продолжаются, определяются разности

. (32)

Полученные разности принимаются за новые значения экспериментальной последовательности для следующей ступени. Аналитическое выражение полинома последовательно накапливается в переменной y в соответствии с рекуррентной суммой

. (33)

Предварительно переменную y полагают равной нулю.

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

5. Проверка функционирования алгоритма

Для проверки правильности работы алгоритма и точности расчёта было решено создать искусственную последовательность данных с заранее известными спектральными характеристиками. Для этого был использован

линейный тригонометрический полином

(33)

Этот полином соответствует формуле (8) .

В полиноме количество слагаемых m=3; набор частот =;

векторы коэффициентов.

Матрица спектра содержит следующие значения:

.

Последний столбец - это значения амплитуд составляющих полинома.

В соответствии с формулой (33) для дискретных значений

при n=205 и h=0.1 была вычислена последовательность значений полинома. Эти значения были приняты за исходные данные для последующих расчётов. Расчёты произведены с помощью программы linpol. В результате был получен такой естественный спектр:

Естественный спектр процесса, m=3



psi Sy[j] a[j] b[j] Am[j]

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

Как видим, параметры истинного и вычисленного спектра очень близки между собой, а вычисленный полином почти полностью совпадает с исходными данными (рис.2)

Рис.2. Сопоставление контрольного и вычисленного естественных полиномов.

Вычисленный программой linpol естественный спектр с наложенным на него истинным спектром представлен на рис.3.



Рис.3. Естественный спектр экспериментальной последовательности с наложенным истинным спектром („крестик“-- истинный спектр,“ромбик“--вычисленный спектр ).

Теперь для сравнения по этим же данным построим полином Фурье. В пакете TRIGA это можно осуществить с помощью программы furpol.Полученный Фурье-спектр представлен на рис.3.



Рис.3.Контрольный Фурье спектр с наложенными линиями истинного спектра (`ромбик`- вычисленный спектр, ?+?- истинный спектр).

Как видим, в спектре Фурье появилось много ложных частот, хотя истинные частоты среди них хорошо различимы.

Чтобы наглядно проследить влияние длины записи на параметры спектра, уменьшим количество данных в последовательности до 105 и рассмотрим полученный Фурье спектр (рис.4). Здесь мы видим, что ложные частоты стали значимее, ошибка в определении истинных спектральных линий увеличилась, хотя они всё ещё неплохо различимы. Требуемая точность приближения данных полиномом достигается за счёт использования 17 составляющих, из которых 14 ложных. Все амплитуды близких к истинным частот колебаний несколько меньше истинных амплитуд. За счёт этого высвободилась энергия для ложных частот.



Рис.4. Спектр Фурье при 105 отсчётах с наложенными линиями истинного спектра (обозначены «+»).

Для сравнения ниже представлен естественный спектр, полученный в тех же условиях (рис.5). Для представления данных с той же точностью потребовалось восемь составляющих, из которых пять ложных. Истинные спектральные линии просматриваются хорошо, хотя ошибка в определении частот и амплитуд несколько увеличилась. Тут не лишне напомнить, что программа определения естественного спектра при каждом варианте расчёта даёт несколько отличные результаты. Из этих вариантов можно выбрать наилучший, т.е. такой, который при заданной точности найдёт меньшее количество составляющих. Причина появления ложных составляющих состоит в том, что приближение данных полиномами мы осуществляем во временной области, а о характере приближения судим в частотной области. Однако приближать данные сразу в частотной области мы не можем, так как они заданы именно во временной области.

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

. (34)

Рис.5. Естественный спектр контрольных данных при 105 отсчётах.

В нашем случае, чтобы между частотой 0.956 рад/с и нулём не появлялась ложная частота, длина записи должна быть не менее 13с.



6. Предварительная обработка данных

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

вычисление числовых статистических характеристик;

центрирование и нормирование данных;

обнаружение и удаление тренда;

корректировка выбросов;

сглаживание данных.

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

Определение числовых характеристик выполняется в соответствии с формулами (11 -16), операция центрирования заключается в вычитании среднего значения из всех элементов последовательности. Если среднее значение велико по абсолютной величине, то такая операция повышает точность последующих вычислений. Нормирование данных осуществляется по формуле (17) и кроме центрирования включает деление на стандартное отклонение. Нормирование позволяет сопоставлять последовательности с различными физическими единицами элементов.

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

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

(35).

Для контрольной последовтельности этот полином имеет вид

,

совместный график тренда и контрольной последовательности представлен на рис.6. Такой же график, но для укороченной последовательности показан на рис.7. Как видим, для исходной последовательности тренд мал и не повлиял на вычисление истинных частот спектра. В укороченной последовательности влияние низкочастотного шума усилилось. Это привело к искажению спектра и появлению ложных частот.

Полиномиальный тренд удаляется путём вычитания вычисленных значений тренда из соответствующих значений последовательности. Однако это нужно делать осторожно, так как тренд зависит не только от внутренних тенденций в процессе, но и от длины записи. Выше мы в этом убедились. Определение и удаление тренда можно выполнить с помощью программы vorarb пакета TRIGA. окончательное решение об удалении тренда принимается пользователем присваиванием соответствующего значения управляющей переменной utr.



Рис.6. График контрольной последовательности с наложенным на него

кубическим трендом.

Рис.7. Полиноминальный тренд укороченной контрольной последовательности

Под выбросом мы понимаем аномальное, не согласующееся с остальными элементами, последовательности значение. Для проверки значимости выбросов разработаны специальные статистические критерии. В программе vorarb используется критерий наибольшего абсолютного отклонения [17,c.547]. Критерий основан на статистике

. (36)

Критическое значение критерия с доверительной вероятностью 0.95 вычисляется по формуле

(37)

и зависит от числа элементов последовательности n. При > соответствующее значение признаётся выбросом. Корректировка выбросов осуществляется медианным сглаживанием[18].

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

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

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

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

В программе vorarb процедура сглаживания включает два этапа. Сначала выполняется медианное скользящее сглаживание по пяти точкам [18]

. (38)

Если все наборы по пять чисел упорядочить по возрастанию, то медианой будет третье число. Последовательность медиан {v_i} будет короче исходной последовательности на четыре элемента, для краевых значений сглаживание выполняется по формулам:

(39)

.

На следующем этапе строится последовательность {z_i} по таким соотношениям:

, (40)

.

Для остальных программ пакета последовательность {z_i} рассматривается как исходная, т.е. {x_i}={z_i}.

7. Модификация алгоритма гармонического анализа

Разложение в ряд Фурье или классический гармонический анализ оказываются полезными при решении многих практических задач [19-21].

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

Для дискретных последовательностей непосредственное применение интегральных формул Фурье при определении коэффициентов не представляется возможным. Поэтом они вычисляются численными методами по следующим формулам:

,

, (41)

,

где m - наибольшее целое число, которое меньше или равно n/2,. Для дальнейшего использования оставляют k первых коэффициентов.

Формулы (41) соответствуют методу левых прямоугольников. Используются также и более точные методы: трапеций и Симпсона. Однако, сравнительный анализ показал, что для наших экспериментальных данных они не дают ощутимого увеличения точности .

Выполнение гармонического анализа в данном пакете осуществляется с помощью программы furpol. Первая модификация заключается в том, что после вычисления каждой пары коэффициентов проверяется по критерию (23) относительная погрешность аппроксимации. Если достигнуто значение критерия, которое меньше заданной погрешности , то вычисление коэффициентов Фурье прекращается. Таким образом, во многих случаях удаётся сократить время вычислений и количество данных, подлежащих анализу. Следующая важная модификация касается случаев, когда точно известен основной период исследуемого процесса , и последовательность содержит несколько таких периодов. B качестве примеров можно привести записи колебаний индикаторного давления в цилиндре двигателя внутреннего сгорания и крутящего момента в зависимости от угла поворота коленчатого вала. Угол принято измерять в градусах. Период колебаний индикаторного давления в 4-х тактном двигателе составляет 720°. Период колебаний крутящего момента 8-ми цилиндрового двигателя составляет 90°. Гармонический анализ момента приводится ниже.

Формальное применение разложения (41) приводит к появлению большого количества незначимых частот в спектре процесса, которые затемняют истинную картину. Для устранения этого явления, необходимо просмотр частот выполнять не от величины , где - длина записи, а от величины , гдеТ1-известный период колебаний. Дальше каждая гармоническая составляющая,проверяется на значимость вклада в общую вариацию. Если, (42)

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

8. Практические примеры расчёта

8.1 Тяговое сопротивление плуга

В качестве исходных данных для примера была использована осциллограмма, записанная при полевых испытаниях пахотного сельскохозяйственного агрегата в составе трактора Т-74 и плуга П-5-35МГ. На осциллограмму записывались тяговое сопротивление плуга, крутящий момент и угловая частота вращения двигателя СМД-14А. Последующая обработка данных показана для тягового сопротивления, которое измеряется в кН. Отметки времени на осциллограмме проставлялись в виде тонких вертикальных линий через 0.1 с.

Значения тягового сопротивления были сняты с осциллограммы через 0.1с. при помощи мерной линейки с ценой деления 1мм и записаны в журнал, а в последствии перенесены на машинный носитель. Полученный текстовый файл выглядит следующим образом:

55 53 57 51 55 54 52 50 50 48 49 48 46 46 45 51 49 51 54 40

46 54 47 42 48 45 45 54 55 50 50 52 52 52 51 54 44 49 46 47

52 49 47 54 57 45 44 47 46 48 50 43 43 49 46 47 54 54 55 45

48 47 52 54 52 50 54 52 54 57 57 53 49 57 52 47 54 50 48 50

48 48 53 50 49 48 49 50 47 44 50 49 52 56 52 51 56 61 54 56

57 48 47 50 50 49 51 55 59 59 53 52 57 57 48 55 52 57 56 52

57 59 57 60 53 51 52 52 57 45 50 51 56 55 59 54 58 64 66 63

60 65 63 55 64 67 56 57 56 52 46 52 45 37 41 48 45 47 45 40

46 45 47 44 48 45 45 48 52 52 52 54 55 49 51 48 50 52 54 53

49 54 55 57 55 56 55 61 56 53 55 52 52 57 61 53 49 48 48 44

46 48 50 51 54 51 48 47 51 54 55 49 46 49 51 49 46 45 49 50

46 38 40 45 43 44 45 777

0.1 0.3345 0

Данные в файле отделены друг от друга пробелами. В файле 13 строк, в каждой строке кроме двух последних по 20 чисел. Количество чисел в строке устанавливается заранее и учитывается в операторе программы readdata. В предпоследней строке данных может быть меньше, в данном случае их 8, последнее число «777» означает конец строки и всей экспериментальной последовательности. Последняя строка содержит служебную информацию: шаг фиксации данных по времени и два масштабных коэффициента для перевода мм в кН. Второй масштабный коэффициент в данном случае равен нулю, поскольку базовая линия отсчёта и нулевая линия для тягового сопротивления совпадают. Для угловой скорости эти линии не совпадают. Отсчёты выполняются не от нуля, и второй коэффициент равен 140 рад/с.

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

Имена и смысл переменных при расчётах:

NF-дескриптор текстового файла данных без расширения, единицы измерений данных [ed] и независимой переменной [et], nf-количество данных в строках ,

K-коэффициент для определения максимальной частоты диапазона поиска, cod-код конца данных, h-шаг по аргументу между дискретными значениями данных.

Имена файла для сохранения результатов расчёта и файла исходных данных:

Значения переменных:

.



Количество строк и общее количество данных n в файле: nsd=13, n=227.

Шаг по независимой переменной , масштабные коэффициенты :

.

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

.

Значения данных после масштабирования:



Определены статистические характеристики исходных данных: среднее значение, дисперсия, стандартное отклонение, размах, оценки асимметрии и эксцесса, вариация. При расчёте они выглядят так:

.

Тренд определялся методом наименьших квадратов в виде полинома третьей степени. В результате расчёта найдено такое уравнение тренда

.

График приведен ниже (рис.8).

Рис 8. Тяговое сопротивление плуга с наложенным трендом.

Удаление тренда привело к снижению дисперсии на 15.2%. Статистические характеристики данных после удаления тренда:

Далее по критерию наибольшего абсолютного отклонения [17] проверялось наличие выбросов в последовательности. По данному критерию выбросов не найдено. Далее выполнялось сглаживание тягового сопротивления плуга (рис.9).

Рис.9. Медианное сглаживание тягового сопротивления плуга.

На рис.9 отчётливо просматриваются высокочастотные колебания, которые наложились на основные колебания сопротивления почвы. Эти колебания могут быть вызваны вибрацией тензозвена, установленного между плугом и трактором. Как видим, медианное сглаживание легко отфильтровывает эти колебания. Статистические характеристики последовательности после сглаживания:

.

Взаимная корреляционная функция исходных и сглаженных данных

.

Результаты работы программы предварительной обработки передаются остальным программам пакета.

Аппроксимация тягового сопротивления линейным тригонометрическим полиномом выполнена программой linpol из пакета TRIGA. Последовательный процесс аппроксимации показан ниже

Здесь идентификаторы программы имеют следующий смысл: psi-почти частота, Sx.-вариация исходных данных на соответствующей ступени, Sr-остаточная вариация,Sy-вариация составляющей,Sd-накопленная сумма вариаций всех найденных составляющих, Er-текущая относительная погрешность аппроксимации. Как видим относительная погрешность разложения экспериментальной последовательности на 16 составляющих вида(8) равна 0.05.

Найденному тригонометрическому полиному соответствует следующий спектр. В матрице спектра по столбцам записаны следующие параметры: почти частота, вариация составляющей, коэффициенты полинома и амплитуда. Графическое представление амплитудного естественного спектра показано на рис.10. Из рисунка хорошо видно, что в отличие от спектра Фурье расстояние между спектральными линиями здесь неодинаково. По характеру спектра можно предположить, что исследуемая функция представляет собой сумму из трёх не связанных между собой процессов с диапазонами частот 0..3,4..6 и 7.8 рад/с.

Естественный спектр процесса ,m=16,j=1..16

psi Sy[j] a[j] b[j] Am[j]

Уравнение тригонометрического полинома представляет собой сумму из 32-х гармоник с естественными частотами и свободного члена

Рис.10. Амплитудный естественный спектр тягового сопротивления агрегата.

Точность аппроксимации легко проследить по графику (рис.11).



Рис.11.Аппроксимация тягового сопротивления пахотного агрегата линейным тригонометрическим полиномом с естественными частотами.

Из рисунка 10 хорошо видно, что в отличие от спектра Фурье расстояние между спектральными линиями здесь неодинаково. По характеру спектра можно предположить, что исследуемая функция представляет собой сумму из трёх практически не связанных между собой процессов с диапазонами частот 0..3,4..6 и 7.8 рад/с.

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

Естественный кубический спектр процесса ,m=12,j=1..12



psi Sy[j] Sd[j] № п./п.

Из общей вариации исходных данных 325 кН² кубическим полиномом с 12-ю составляющими объяснено 310.8 кН² или 95.6%. Нормированный естественный спектр, определённый по кубическим составляющим представлен на рис.12.

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



Рис.12.Нормированый естественный спектр на базе кубического тригонометрического полинома. Нормировка по формуле Sn[j]=Sy[j]/Sd, Sd=310.76.

Для сравнения по этим же данным построим спектр Фурье(рис.13). Одинаковая с естественным полиномом точность 5% достигается при 19 гармонических слагаемых. Шаг по частоте для всех спектральных линий одинаков и равен 0.278рад/сек. Структура процесса просматривается слабо. Характер аппроксимации показан нарис.14.

Рис.13. Гармонический спектр тягового сопротивления плуга.



Рис.14. Аппроксимация тягового сопротивления отрезком ряда Фурье из 19 слагаемых.

8.2. Крутящий момент

Дискретная последовательность значений крутящего момента двигателя зафиксирована на интервале [0,360°] угла поворота вала с шагом 10°. Входной файл

имеет следующий вид

465 1660 2050 2185 1495 885 555 255 170

465 1660 2050 2185 1495 885 555 255 170

465 1660 2050 2185 1495 885 555 255 170

465 1660 2050 2185 1495 885 555 255 170

465 777

10 0.001 0.0

В четырёх первых строках находится по 9 значений, в предпоследней строке - одно значение и код конца последовательности. В предпоследней строке помещён шаг и масштабные коэффициенты для перевода Нм в кНм. Общее количество данных n=37.

Ниже приводится полная длина записи, основной период колебаний момента, разрешающая способность по частоте, основная частота, максимальная частота и допустимое количество частот

.

Числовые характеристики данных таковы:

При использовании модифицированной процедуры расчёта получены полином и спектр Фурье:

,

Спектр момента,m=2

Частота коэффициенты амплитуда вклад сумма вкладов

.

Относительная ошибка аппроксимации 0.0143, остаточная вариация 0.2845 при исходной вариации 19.85. Для аппроксимации момента двигателя достаточно двух гармонических составляющих.

Сравним эти данные с результатами программы linpol вычисления естественного спектра. Здесь мы имеем следующую процедуру вычислений

Программе вычисления естественного спектра не даётся информация о периоде колебаний момента, она его находит сама с некоторой погрешностью. Так, при периоде 90° основная частота момента 0.0698. Найденная программой linpol основная частота равна 0.0704. Относительная погрешность естественной аппроксимации 0.0192, т.е. больше, чем у ряда Фурье.

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



Рис.15.Аппроксимия момента двигателя гармоническим полиномом с двумяслагаемыми.

8.3 Анализ солнечных пятен

Одним из интересных приложений алгоритма естественной аппроксимации является поиск периодичности числа солнечных пятен. Данные о среднем количестве солнечных пятен в январе по годам с 1841 по 1984 были взяты из книги [22] и помещены в текстовый файл C:/MAPLE/TRIGA/DAT/sp.txt:

24.0 20.4 13.3 9.4 25.7 38.7 62.6 159.1 156.7 78.0

75.5 68.4 41.1 15.4 12.3 0.5 13.7 39.0 83.7 81.5

62.3 63.1 48.3 57.7 48.7 31.6 0.0 15.6 60.9 77.3

88.3 79.5 86.7 60.8 14.6 14.3 24.4 3.3 0.8 24.0

36.4 45.0 60.6 91.5 42.8 29.9 10.3 12.7 0.8 5.3

13.5 69.1 75.0 83.2 63.3 29.0 40.6 30.2 19.5 9.4

0.2 5.2 8.3 31.6 54.8 45.5 76.4 39.2 56.7 26.4

3.4 0.3 2.3 2.8 23.0 45.3 74.7 96.0 48.1 51.1

31.5 11.8 4.5 0.5 5.5 71.8 81.6 83.5 68.9 65.3

14.6 12.112.3 3.4 18.9 62.8 132.5 98.4 80.3 50.5

45.6 35.612.4 3.7 18.5 47.6 115.7 108.5 119.1 101.6

59.9 40.7 26.5 0.2 23.1 73.6 165.0 202.5 217.4 146.3

57.9 38.7 19.8 15.3 17.5 28.2 110.9 121.8 104.3 111.5

91.3 61.5 43.4 27.6 18.9 8.1 16.4 51.9 166.6 159.6

114.0 111.2 84.3 57.0 777

1 1 0

В предпоследней строке код 777 обозначает конец основных данных, шаг фиксации 1, первый масштабный коэффициент 1, второй 0. Общее количество данных 144. График среднего числа солнечных пятен совместно с выявленным трендом представлен на рис.16.

Рис.16. График колебаний числа солнечных пятен за 144 года совместно с выявленным трендом этой величины.



Как видно из графика минимум солнечной активности приходился на период с 1885 года по 1901 год. Наличие тренда в экспериментальной последовательности свидетельствует о том, что в солнечной активности имеется период, значительно превосходящий период наблюдений.

Анализ спектра последовательности пятен (рис.17) показывает, что наиболее интенсивные колебания приходятся на диапазон периодов 9.0 - 12.1 лет.

Рис.17.Естественный спектр солнечных пятен за период от 1841 до 1984 годов.

Наиболее существенное колебание имеет период 10.7 лет, следующее по значимости колебание за пределами отмеченного диапазона имеет период 21.7 года. Существенные колебания заметны в области со средним периодом 5.4 года. Могут представлять интерес колебания с периодами 54.5 и 71.1 года.



Заключение

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

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

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

Частотный спектр, установленный с помощью естественного полинома также является естественным . Понятие естественного спектра введено Е.В.

Дмитриевым в работе [1] для полиномов с гармоническими составляющими. В данном документе это понятие используется в несколько расширенном формате и применяется не только для гармонических, но и для более сложных составляющих.

Наилучшую точность и частотное разрешение разработанные алгоритмы могут обеспечить когда исследуемый процесс предсталяет собой

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

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

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

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

Литература

1. Дмитриев Е.В. Гармонические спектры и аппроксимация коротких сигналов. Воронеж: http://short-signal-sp.ru, 2005.

2. Хедли Дж. Нелинейное и динамическое программирование. М.:МИР,1967.-506с.

3. Гутер Р.С. и др. Элементы теории функций. Функции действительного переменного. Приближение функций. Почти - периодические функции. М.: ФИЗМАТГИЗ,1963.-244 с.

4.Бендат Дж., Пирсол А. Измерение и анализ случайных процессов. - М.: МИР,1971.408с.

5. Анго А. Математика для электро - и радиоинженеров. М .:1965.780с.

6. Почти периодическая функция. БСЭ. http://bse.sci-lib.com/article092130.html

7. Большев Л.Н., Смирнов Н.В. Таблицы математической статистики. М.: Наука, Гл. ред. физматлит,1983.

8. Хан Г., Шапиро С. Статистические модели в инженерных задачах. М.: МИР,1969.

9. ТЕХНИЧЕСКАЯ КИБЕРНЕТИКА. Теория автоматического регулирования. Книга 2. М.: Машиностроение, 1967.-637с.

10. Котельников В.А. О пропускной способности эфира и проволоки в электросвязи. - Успехи физических наук Т.176,номер 7,2006,с.762-769.

11. Виленкин С.Я. Статистические методы исследования систем математического моделирования. М.: СОВЕТСКОЕ РАДИО,1967.-200с.

12. Дрейпер Н., Смит Г. Прикладной регрессионный анализ. М.:Статистика,1973. -392с.

13. Бард Й. Нелинейное оценивание параметров. М.: Статистика,1979.-349с.