Быстрое преобразование Фурье
Высокая степень развития рационального интеллекта у человека привела к значительной и печальной потере его интуиции, которая так необходима сегодня для дальнейшего развития.
Жак Моно, "Случайность и необходимость".
Введение
Речь пойдет о давно известных вещах, о затертых, вывалянных в пыли понятиях.
Стандартные методы спектрального оценивания с использованием быстрого вычислительного алгоритма, именуемого быстрое преобразование Фурье (БПФ), основаны на модели представления данных с помощью рядов Фурье, то есть анализируемый процесс полагается состоящим из некоторого набора гармонически связанных синусоид.
В нетехнических областях в течении многих лет используются модели на основе других временных рядов. Юл и Уокер использовали АР-модели для прогнозирования трендов во временных рядах, характеризующих экономические данные. Прони предложил процедуру согласования экспоненциальных моделей с данными экспериментов в области химии газов.
Ряд других моделей появился в других областях, связанных со статистическим и численным анализом. Их терминология и используемые алгоритмы весьма разнообразны, поскольку разрабатывались они в самых различных областях, таких как анализ речевых сигналов, обработка изображений, геосейсмический анализ и радиолокация. Все это многообразие методов и алгоритмов, в сочетании с нелинейными моделями анализа на основе современных вычислительных технологий ждет своего практического применения в различного рода информационных системах.
Рассматриваемый ниже компонент – маленький шаг в этом направлении.
Библиотека компонентов для Delphi, реализующая линейную обработку данных. Включает методы спектрального и последовательного анализа. Библиотека компонентов поставляется с исходными текстами. Способ распространения - бесплатно для некоммерческого применения.
Компонент, реализующий низкочастотную фильтрацию. Позволяет анализировать всплески (выбросы), сглаживать и подавлять шумы в исходном сигнале. Реализован линейный последовательный алгоритм оценивания состояния, обеспечивающий минимизацию среднеквадратичной ошибки (фильтр Калмана).
Компонент работает в двух режимах. Процесс получения текущих оценок в первом режиме получил название скользящего окна. В этом случае по сути решается переопределенная система линейных уравнений стандартными методами линейной алгебры. В принципе, этого достаточно для задачи адаптивной фильтрации. Однако это приводит к значительному объему вычислений, так как коэффициенты фильтра каждый раз пересчитываются с самого начала. Более предпочтителен, с точки зрения эффективности вычислений, второй режим. В этом случае процесс оценивания реализуется как рекуррентная процедура, являющаяся одним из частных случаев алгоритмов фильтрации, известных как фильтры Калмана.
Компонент LowPassFilter
Компонент осуществляет фильтрацию временного ряда в частотной области с использованием классического периодограммного метода спектрального анализа. Для преобразования в частотное представление и обратно используется эффективный вычислительный алгоритм БПФ (быстрое преобразование Фурье).
Компонент может работать в трех режимах: нет фильтрации, фильтрация нижних частот и режим вычитания шума. Присваивая свойству Mode значение mdNoneFiltering, mdLowPassFilter или mdSubtractionNoise соответственно, можно задать эти режимы.
В первом случае все данные проходят со входа на выход без преобразования спектра. Во втором случае на выход проходит низкочастотная часть спектра, указанная в свойстве BandWidth.
В режиме вычитания шума на основе анализа спектра входных данных принимается решение о принадлежности спектральной составляющей к шумовой или к полезной составляющей. На выход проходят лишь составляющие, отнесенные к полезному сигналу.
В дополнение к перечисленному, в зависимости от значения свойства Overshoot (true или false) компонент осуществляет редактирование неправдоподобных значений во входных данных.
Борьба с неправдоподобными значениями
Большинство систем сбора информации вносит в данные ложные значения. Это может происходить по многим причинам, например, в результате технических сбоев, или ошибки оператора, или еще чего-либо. Неправдоподобные значения, возникающие в результате этих сбоев, могут вызвать значительные трудности при последующем анализе или сделать его совсем невозможным.
Хорошо бы на этапе предобработки данных найти и удалить такие значения. К сожалению, отделить неправдоподобные данные от правдоподобных трудно. Общих процедур для автоматического удаления таких данных нет. Существующие схемы не являются полностью удовлетворительными и, в той или иной мере, основываются на некоторых априорных предположениях о характере данных.
Алгоритм подавления неправдоподобных выбросов, реализованный в данном компоненте, основывается на, так называемой, робастной (устойчивой) оценке средней скорости изменения временного ряда. В качестве последней используется медиана.
Если все данные упорядочить в порядке возрастания, то медианой будет то значение, которое в этом ряду будет посередине.
Для пояснения – хрестоматийный пример. Пусть есть данные о росте 100 человек. Определим среднее и медиану. Среднее вычисляется обычным путем, а медиана – при помощи упорядочивания результатов в порядке возрастания (или убывания) и выборе значения, оказавшегося посередине.
Если при вводе значений не было больших ошибок, то эти значения окажутся довольно близкими. Если же при вводе одного из значений произошла ошибка, например на 10 метров, то среднее возрастет на 10 см, в то время как медиана останется практически неизменной.
В алгоритме это свойство робастности медианы используется дважды. Исходный ряд данных $X(n)$ дифференцируется и получается ряд $abs(X′(n))$. Вычисляется медиана полученного ряда, то есть, получаем медиану абсолютных скоростей изменения исходного ряда – $V_x$. Таким же образом вычисляется медиана для ряда $abs(abs(X′(n)) – V_x)$, то есть, получаем медиану абсолютных ускорений изменения исходного ряда – $A_x$.
Свойство компонента SuppressionDegree задает степень подавления выбросов. В зависимости от значения этого свойства (edSmall, edMedium, edLarge) изменения соседних выборок по абсолютной величине ограничиваются соответственно значениями $V_x + 4 * A_x, V_x + 3 * A_x, V_x + 2 * A_x$.
Подавление неправдоподобных выбросов можно включать или отключать, присваивая свойству Overshoot значение true или false соответственно.
Некоторые детали
Известные алгоритмы БПФ имеют наиболее простой вид, когда количество выборок во входных данных есть целая степень двойки. Если это не так, поступают просто: добавляют недостающее количество нулей до ближайшей степени двойки. Более того, чтобы увеличить частотное разрешение спектральных составляющих, например, в два раза, исходный временной ряд дополняют таким же количеством нулей.
При этом вычисляемые спектральные составляющие располагаются в два раза чаще. Понятно, что таким же образом можно и утроить спектральное разрешение и увеличить его в 4 раза и т.д. Чтобы управлять эти процессом служит свойство FrequencyResolution. Допустимые значения – целые числа от 1 до 10.
Необходимо помнить о том, что задавая большое разрешение, нужно быть готовым к тому, что процедура преобразования в частотную область и обратно может занять продолжительное время (особенно при больших выборках исходного временного ряда). Поэтому для контроля процесса рекомендуется отображать ход его выполнения с помощью индикатора выполнения ProgressBar.
Для этого необходимо обрабатывать событие OnStepExecute с параметром типа double, равным проценту выполненной работы. Для инициализации ProgressBar можно привязаться к событию BeforeExecute. Чтобы прервать обработку, существует метод BreakExecute.
Низкочастотная фильтрация
В режиме фильтрации нижних частот можно ограничить верхнюю частоту спектра, пропускаемую на выход, задавая свойство BandWidth от 0 до 1. При этом меняется степень сглаживания исходного процесса.
Пусть, например, имеется временной ряд с данными о ежедневных объемах продаж в течении какого то времени. Чтобы выделить тенденции, характеризующие изменение еженедельных объемов продаж, необходимо отфильтровать исходный ряд, задав BandWidth равным 1/7, ежемесячных – 1/30 и т.д.
Вычитание шума
В этом режиме оценивается распределение спектральных составляющих по амплитуде и полагая, что шумовые составляющие распределены по релеевскому закону вычисляется порог, при превышении которого составляющая проходит на выход фильтра, иначе – отбрасывается.
В процессе обработки спектр упорядочивается по амплитудам, строится гистограмма и вычисляется порог. Число интервалов в гистограмме считается по формуле $1+3.2*log(N)$, где $N$ – число спектральных составляющих. Порог вычисляется по эмпирической формуле и соответствует той точке на гистограмме справа от максимума, где значение гистограммы пересекает уровень $K * \mbox {<максимальное значение>}$, где множитель $K$ можно задать равным 0.1, 0.06 или 0.03, присваивая свойству SubstractionNoiseDegree значения edSmall, edMedium и edLarge соответственно.
Результаты работы в этом режиме тем удовлетворительнее, чем больше отсчетов в исходных данных и выше спектральное разрешение.
Во всех трех режимах строится упорядоченный по амплитудам массив составляющих выходного спектра. Размерность массива возвращается в свойстве HostHarmonicCount. Каждый элемент массива доступен через индексированное свойство HostHarmonicItem. Возвращаемое значение равно периоду гармоники, выраженному в количестве точек исходного временного ряда, укладывающихся на этом периоде.
Файлы
Компонент размещен для бесплатного скачивания на сайте. К нему имеется Help-файл и редактор компонента, позволяющий легко редактировать параметры фильтра на этапе разработки. К компоненту прилагается демонстрационный пример.
Файлы имеют следующие имена:
Package файл – SpectralEstimation.dpk
Компонент – SELowPassFilter.pas
Файл помощи – LowPassFilterHelp.hlp, LowPassFilterHelp.cnt
Редактор компонента – SELpfEditor.pas, SELpfEditDlg.pas, SELpfEditDlg.dfm
Демонстрационный пример – TestLpf.pas, TestLpf.dfm
