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

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

ГОСУДАРСТВЕННЫЙ КОМИТЕТ СВЯЗИ, ИНФОРМАТИЗАЦИИ И ТЕЛЕКОММУНИКАЦИОННЫХ ТЕХНОЛОГИЙ РЕСПУБЛИКИ УЗБЕКИСТАН

ТАШКЕНТСКИЙ УНИВЕРСИТЕТ ИНФОРМАЦИОННЫХ ТЕХНОЛОГИЙ

Диссертация

на соискание академической степени магистра

5A330202 - Информационные и мультимедийные технологии

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

Мирзаев Хабибулла Хамидуллаевич

Ташкент - 2013

Введение

Развитие информационно-коммуникационных технологий (ИКТ) является одним из основных факторов благосостояния и экономического роста страны. Сегодня ИКТ становится одним из основных приоритетов государственной политике Узбекистана.

В постановлении президента Республики Узбекистана «о мерах по дальнейшему внедрения и развитию информационно - коммуникационных технологии» (собрание законодательства Республики Узбекистана, 2012 г., № 13, ст. 139) одной из основных задач дальнейшего внедрения и развитие информационно - коммуникационных технологий, в частности, является программа мер по коренному и качественному улучшению функционирование национальной информационной-поисковой системы, увеличению количество ее пользователей.

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

Растущая потребность в современных цифровых системах обработки данных с более чем одной частотой дискретизации привела к развитию новой подобласти цифровой обработки сигналов, известной как обработка данных при нескольких скоростях (multirate processing). При такой обработке данных используются две основные операции: децимация и интерполяция, позволяющие эффективно чередовать скорости передачи данных. Децимация уменьшает частоту дискретизации, эффективно сжимая данные и оставляя только необходимую информацию. Интерполяция, наоборот, увеличивает частоту дискретизации [13].

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

Актуальность данной работы. Преобразование частоты дискретизации неизбежно при обработке сигналов в реальном масштабе времени, когда два процессора, работающие на разных частотах, должны обмениваться цифровыми сигналами. Изменения частоты дискретизации необходимо и при обработке цифровой аудио сигнала, когда надо согласовать частоту дискретизации со стандарта Audio CD 44,1 кГц на стандарт DVD 44,0 кГц.

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

Цель и задачи исследования

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

Для достижения поставленной цели в диссертационной работе решаются следующие задачи:

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

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

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

- создание прикладных программ изменения частоты дискретизации аудиосигналов для применения в сигнальных процессоров.

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

Научная новизна проведенных исследований заключается в том, что:

- в результате исследования разработан алгоритм перевода сигнала в область полиномиального представления и нахождения коэффициентов полинома в базисной системе Адамара;

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

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

Структура и объём магистерской диссертационной работы

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

Глава 1. Изучение современных методов цифровой обработки сигналов

1.1 Методы изменения частоты дискретизации

Простой, но слишком наивный подход к изменению частоты дискретизации цифрового сигнала -- это вернуть его к аналоговому виду, а затем повторно оцифровать его с иной частотой [13].

Рис. 1.1. Процесс изменения частоты дискретизации классической методом.

В этом методе дискретный сигнал с частотой дискретизацией fs₁ передаётся в блок ЦАП (цифро-аналоговый преобразователь), и в этом блоке дискретный сигнал переводится в аналоговую форму f(t). После этого аналоговый сигнал f(t) передается в блок АЦП и там сигнал дискретизируется с нужной частотой дискретизацией fs₂.

При таком подходе изменения частоты дискретизации, теряется много времени за счет того, что сигнал передается с блоков ЦАП и АЦП. Процессу преобразования цифровой-аналоговый-цифровой присущи такие ошибки, как ошибки квантования и наложения, которые значительно искажают сигнал [13]. дискретизация алгоритм сигнал цифровой

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

При передискретизации отсчёты сигнала, соответствующие одной частоте дискретизации, вычисляются по имеющимся отсчётам этого же сигнала, соответствующим другой частоте дискретизации (при этом предполагается, что обе частоты дискретизации соответствуют условиям теоремы Котельникова). Идеальная передискретизация эквивалентна восстановлению непрерывного сигнала по его отсчётам с последующей дискретизацией его на новой частоте [9, 13, 17].

Точное вычисление значения исходного непрерывного сигнала в определённой точке производится следующим образом:

где -- i-й отсчёт сигнала, -- момент времени, соответствующий этому отсчёту, -- циклическая частота дискретизации, -- интерполированное значение сигнала в момент времени t.

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

где h(t) -- импульсная характеристика соответствующего восстанавливающего фильтра. Вид этого фильтра выбирается в зависимости от задачи [7].

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

- децимация с целым коэффициентом (уменьшение частоты дискретизации в целое число раз);

- интерполяция с целым коэффициентом (увеличение частоты дискретизации в целое число раз);

- изменение частоты дискретизации в рациональное (M/N) число раз (этот случай можно рассматривать как комбинацию двух предыдущих) [7, 16].

Децимация с целым коэффициентом.

На рис. 1.2. представлена блок-схема процесса децимации сигнала x(n) с целым шагом М. На ней изображены цифровой фильтр защиты от наложения спектров h(k) и схема сжатия (компрессор) частоты дискретизации, обозначенная стрелкой, направленной вниз, и коэффициентом децимации М. Компрессор частоты дискретизации снижает частоту дискретизации с Fs до Fs/M. Для предотвращения наложения при более низкой частоте выборки используется цифровой фильтр, предварительно ограничивающий полосу входного сигнала до FS/2M. Таким образом, сигнал х(п) вначале ограничивается по полосе. Снижение частоты дискретизации достигается за счет того, что из каждых М выборок фильтруемого сигнала w(n) отбрасывается М-1 выборка. Вход и выход процесса децимации связаны следующим соотношением:

На рис. 1.3 описанный процесс иллюстрируется для простого случая М=3, где из каждых трех выборок х(п) отбрасываются две. Отметим, что, по сути, децимация -- это операция сжатия данных [9,12]

Рис. 1.2. Блок схема децимации с шагом M

Рис. 1.3. Иллюстрация децимации с шагом M=3 во временной области.

Спектральное представление процесса децимации приведено на рис. 1.4. где предполагается, что на вход подается широкополосный сигнал х(n). На рис. 1.4. б пунктиром изображены зеркальные компоненты сигнала, которые привели бы к наложению, если бы входной сигнал х(n) не был ограничен по полосе перед децимацией [9, 12, 13].

Интерполяция с целым коэффициентом

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

Рис. 1.4. Спектральная интерпретация децимации сигнала с частоты 6 до частоты 2 кГц.

Для данного сигнала х(п) с частотой дискретизации Fs процесс интерполяции увеличивает частоту дискретизации в L раз, т.е. до LF_s. Схема интерполятора приведена на рис. 1.5. Он состоит из экспандера частоты дискретизации, обозначенного стрелкой, направленной вверх, и коэффициентом интерполяции L, который показывает, во сколько раз увеличивается частота дискретизации [9, 12].

Рис. 1.5. Блок схема интерполяции с шагом L=3

Рис. 1.6. Иллюстрация интерполяции с шагом L=3 во временной области.

Для каждой выборки сигнала х(п) экспандер вводит (L--1) нулевую выборку, в результате чего формируется новый сигнал w(m) с частотой дискретизации LF_S. Далее этот сигнал пропускается через фильтр нижних частот для удаления зеркальных частот, введенных при увеличении частоты, и получается сигнал у(т).

Рис. 1.7. Спектральная интерпретация интерполяции сигнала с 2 до 6 кГц

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

Процесс интерполяции во временной области иллюстрируется на рис.1.6. для простого случая L = 3. Обратите внимание на то, что каждая выборка х(n) порождает три выходные выборки (две нулевые выборки вводит экспандер) [9, 13].

Интерпретация данного процесса в частотной области приведена на рис.1.7. Функции X(f), W(f) и Y(f)) представляют собой частотные характеристики сигналов х(n), w(m) и у(т) соответственно. H(f) -- это амплитудная характеристика фильтра подавления зеркальных частот. Данный фильтр необходим для удаления зеркальных компонентов, обозначенных в W(f) пунктиром. На данном этапе стоит отметить, что процессы децимации и интерполяции дуальны, т.е. взаимно обратны. Данное свойство дуальности означает, что интерполятор можно легко получить из эквивалентного дециматора наоборот.

Преобразование частоты дискретизации с нецелым шагом

В некоторых ситуациях часто бывает нужно изменить частоту дискретизации в нецелое число раз. Пример -- цифровое аудио, где может требоваться передача данных с одного запоминающего устройства на другое, причем системы поддерживают разные частоты дискретизации. В частности, это передача данных с системы воспроизведения компакт- дисков (44,1 кГц) на цифровую аудиопленку (digital audio tape -- DAT) (48 кГц). Для выполнения требуемого преобразования частоту дискретизации компакт-диска следует увеличить в 48/44,1 раза [13].

На практике подобные нецелые множители представляются рациональным числом, т.е. отношением двух целых чисел, скажем L и М, таких что L/М максимально близко к желаемому множителю. Необходимое же преобразование частоты дискретизации производится в два этапа: 1) интерполяция данных с шагом L и 2) децимация с шагом М (рис. 1.8, а). Необходимо, чтобы процесс интерполяции предшествовал децимации, поскольку в противном случае при децимации исчезнут некоторые необходимые частотные компоненты. В приведенном выше примере “с компакт-диска на цифровую аудиопленку” преобразования частоты в 48/44,1 раза можно достичь следующим образом: интерполировать с шагом L -- 160, а затем провести децимацию с шагом М -- 147, т.е. вначале скорость данных компакт-диска увеличивается в L = 160 раз до 7056 кГц, а затем уменьшается в М = 147 раз до 48 кГц.

Поскольку два фильтра нижних частот на рис. 1.8, a h1(k) и h2(k) соединены каскадно и имеют общую частоту дискретизации, их можно объединить, и тогда получится один обобщенный конвертер частоты дискретизации, изображенный на рис. 1.8, б. Если М > L, операция, производимая конвертером, называется децимацией с нецелым шагом, а если М < L -- интерполяцией. Если М -- 1, действие обобщенной системы сводится к простой интерполяции с целым шагом, описанной ранее, а при L = 1 оно сводится к децимации с целым шагом [9, 13].



Рис. 1.8. Иллюстрация интерполяции с рациональным шагом (L = 3, М = 2)

На рис. 1.8, в иллюстрируется интерполяция с шагом 3/2. Частота дискретизации вначале увеличивается в 3 раза (к каждой выборке х(n) добавляются две нулевые выборки), а затем сигнал пропускается через фильтр нижних частот, результат -- v(i).

Далее фильтрованные данные сокращаются в два раза, т.е. из каждых двух выборок v(i) остается одна. Иллюстрацией данного процесса в частотной области служит рис. 1.9. Входной сигнал х(n) с частотой дискретизации 2 кГц вначале ускоряется в 3 раза до 6 кГц, фильтруется с целью устранения зеркальных частот, которые иначе вызвали бы наложение, а затем замедляется в 2 раза до частоты 3 кГц [16].



Рис. 1.9. Спектральная интерпретация увеличения частоты дискретизации сигнала с 2 кГц в 3/2 раз

Недостатком данного метода является необходимость фильтрации сигнала на повышенной в M раз частоте дискретизации, что требует значительных вычислительных ресурсов. При этом соответствующая частота может во много раз превосходить как исходную, так и окончательную частоту передискретизации, особенно если M и N -- близкие большие числа. Так, например, при передискретизации звукового сигнала с 44100 Гц до 48000 Гц этим методом необходимо увеличить частоту дискретизации в 160 раз до 7056000 Гц и затем уменьшить её в 147 раз до 48000 Гц. Таким образом, в данном примере вычисления приходится производить на частоте дискретизации более 7 МГц [17].

1.2 Анализ классических методов обработки сигналов во временной области

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

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

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

Приближение функций интерполяционными полиномами

Интерполяция, интерполирование -- в вычислительной математике способ нахождения промежуточных значений величины по имеющемуся дискретному набору известных значений [5, 7].

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

Существует также близкая к интерполяции задача, которая заключается в аппроксимации какой-либо сложной функции другой, более простой функцией. Если некоторая функция слишком сложна для производительных вычислений, можно попытаться вычислить её значение в нескольких точках, а по ним построить, то есть интерполировать, более простую функцию. Разумеется, использование упрощенной функции не позволяет получить такие же точные результаты, какие давала бы первоначальная функция. Но в некоторых классах задач достигнутый выигрыш в простоте и скорости вычислений может перевесить получаемую погрешность в результатах [6, 8].

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

Задача интерполяции функции, интерполяционные полиномы.

Пусть на отрезке [a,b] задана функция ѓ(x). Задача интерполяции (или интерполирования ) состоит в построении функции g(x), совпадающей с заданной ѓ(x) в некотором наборе точек {x₀,x₁,...,x_n} из отрезка [a,b] (эти точки называются узлами интерполяции), т.е. должны выполняться условия:

g(x_k)=y_k, k=0,1,...,n,

где y_k - известные значения функции ѓ(x) в точках x_k. Функция g(x) называется интерполянтом функции ѓ(x).

Пример интерполяции с четырьмя узлами приведен на следующем 1.1 рисунке, из которого видно, что узлы интерполяции не обязательно должны располагаться равномерно на отрезке [a,b][6, 19].

Если ѓ(x) табличная функция, скажем полученная из эксперимента, т.е. известны только ее значения y_k в точках x_k, то, вообще говоря, о качестве полученного приближения судить трудно. Однако, если значения ѓ(x) могут быть вычислены в любой точке отрезка [a,b], то в этом случае можно исследовать качество получающегося приближения, например найдя максимальное уклонение функции g(x) от функции ѓ(x). На качество приближения сильное влияние оказывает количество и расположение узлов, а также гладкость функции ѓ(x).

Рис.1.10. Иллюстрация интерполяции с четырьмя узлами.

Мы рассмотрим только линейную интерполяцию, т.е. такую, при которой функция g(x) разыскивается в виде линейной комбинации некоторых функций

где для k=0,1,...,n: ц_k(x) - заданные функции, а a_k - искомые коэффициенты[4, 19].

Ясно, что из постановки задачи интерполяции (т.е. из совпадения значений интерполянта g(x) и интерполируемой функции ѓ(x) в точках x_k) следует, что коэффициенты a_k определяются из решения следующей системы линейных алгебраических уравнений:

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



Очень часто в качестве системы функций ц_k(x) выбирают полиномы, например, степени x, именно:

Тогда соответствующий интерполянт называют интерполяционным полиномом. Существование и единственность интерполяционного полинома гарантируется, если все узлы интерполяции x_k различны.

Действительно, определитель системы линейных алгебраических уравнений для нахождения коэффициентов a_k является определителем Вандермонда, который, как известно, равен выражению 1.8 и, следовательно, отличен от нуля в случае, когда все узлы x_k различны. Поскольку матрица системы невырождена, то решение системы существует и единственно. Итак, интерполяционный полином существует и единственный [3, 4, 19].

Можно рассуждать и по-другому. Предположим, что есть два интерполяционных полинома g_k и h_k степени n такие, что для произвольного набора значений выполняются равенства g(x_k) = h(x_k) = y_k для всех k=0,1,...,n, т.е. для n-ой точки. Тогда их разность h_k - g_k является полиномом степени не выше n, но обращается в ноль в n-ой точке. По известной теореме алгебры у полинома степени n не может быть больше чем n корней, следовательно, h_k - g_k ? 0 и h_k ? g_k. Единственность установлена. А так как полином единственный, то у соответствующей системы линейных алгебраических уравнений есть только одно решение (для произвольной правой части). Из результатов линейной алгебры известно, что у системы может быть либо бесконечное число решений при некоторых правых частях, либо единственное, для произвольной правой части. Последнее как раз и выполняется [4, 6].

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

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

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

Интерполяционный полином Лагранжа

Итак, мы ищем полином L_n(x) степени не выше n, значения которого совпадают со значениями y_k заданной функции ѓ(x) в узлах x_k, где k=0,1,...,n и все узлы различны.

Одним из способов записи интерполяционного полинома является форма Лагранжа. Предположим, что для k=0,1,...,n функции Ф_n(x) являются полиномами степени n, которые обладают следующим свойством

Тогда полином 1.11 будет как раз тем, который нам и нужен, поскольку это полином степени не выше n и L_n(x_k) = y_k для всех k=0,1,...,n.

Функции Ф_n(x) строятся легко. Действительно, функция 1.12 является полиномом степени n, который обращается в ноль для всех x_j не равных x_k. В точке x_k она принимает значение



Итак, интерполяционный полином в форме Лагранжа имеет вид

или в развернутой форме:

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

Необходимо найти интерполяционный полином по формуле Лагранжа.

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

Рис. 1.11. График найденного полинома

Образование интерполяционных многочленов Лагранжа связано с большой вычислительной работой. Также велика вычислительная работа при получении значения L_n(x) для какого-то фиксированного значения x. Если даже мы имеем интерполяционный многочлен Лагранжа, построенный по значением x₀, x₁, …, x_n, то это мало помогает нам при построении интерполяционного многочлена Лагранжа в других точках измерения заданной функции [7].

Интерполяционный полином Ньютона

Форма Лагранжа записи интерполяционного полинома, не является единственно возможной. Кроме этого, она имеет определенные неудобства. Так предположим, задан набора точек (x_k, y_k) k = 0,1, …, n и по нему построен интерполяционный полином в форме Лагранжа 1.17

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

Важно понимать, что каким бы способом не строился интерполяционный полином по заданному набору точек, результат всегда получается один и тот же (с точностью до ошибок, возникающих при округлении вещественных чисел), поскольку через заданные точки проходит ровно один полином n -ой степени [7, 19].

Для вывода интерполяционного полинома в форме Ньютона понадобится понятие о разделенных разностях.

Пусть нам заданы различные точки x₀, x₁, …, x_n и функция f(x). Первой разделенной разностью функции f(x) в точках x_i, x_i+1 называется величина

Разделенные разности более высоких порядков определяются рекуррентно. Например разделенная разность второго порядка в точках x₀, x₁, x₂ определяется через две разделенные разности первого порядка [x₀, x₁]f и [x₁, x₂]f следующим образом:

Разделенная разность второго порядка в точках x₀, x₁, x₂ определяется так:

В общем случае разделенная разность s-го порядка определяется через две разделенные разности (s-1)-го порядков

Разделенные разности удобно представить в виде треугольной таблицы

Разделенные разности обладают рядом важных свойств. Одно из них -- симметрия, т.е. разделенная разность [x₀, x₁, …, x_m] не меняет своего значения при любой перестановке x₀, x₁, …,x_m. Это свойство следует, например, из представления разделенной разности в виде.

Получим выражение для интерполяционного полинома Ньютона для интерполяции функции f(x) по набору точек x₀, x₁, …, x_n, воспользовавшись формой Лагранжа.

Интерполяционный полином L_n(x) для набора точек x₀, x₁, …, x_n может быть записан в виде Ln(x) = L₀(x) + (L₁(x) - L₀(x)) + …+ L_n-1(x) - L_n(x)) (очевидно, что в правой части все слагаемые, кроме L_n(x) сокращаются). Здесь -- L_n(x) интерполяционный полином, построенный для f(x) по набору точек x₀, x₁, …, x_k. Он записывается в форме Лагранжа следующим образом:

Ясно также, что L₀(x) = f(x₁) т.к. интерполяционный полином нулевой степени, построенный по одной точке, есть константа, значение которой равно f(x₁)[6,8].

Наша цель -- установить следующие равенства:

Где

Т.е., надо доказать, что

Тогда интерполяционный полином L_n(x) приобретает следующую форму

которая называется формой Ньютона интерполяционного полинома или интерполяционным полиномом Ньютона с разделенными разностями.

Преимущества записи интерполяционного полинома Ньютона очевидны. Если в точке x вычислено значение интерполяционного полинома L_n(x), построенного по набору точек , а затем к этому набору добавилась еще одна точка x_n+2, то для получения значения нового интерполяционного полинома L_n+1(x) (степень которого на единицу выше) в точке достаточно добавить к значению L_n(x) слагаемое

Интерполяционные формулы для равноотстоящих узлов

Если узлы интерполяции отстоят друг от друга на одинаковое расстояние, т.е. x₂ -x₁ = x₃-x₂=…=x_n+1-x_n=h то для интерполяционных полиномов в форме Лагранжа и Ньютона можно получить специальные формулы, упрощающие вычисления

Рассмотрим, например, интерполяционный полином в форме Лагранжа

Введем обозначение

Тогда выражение для принимает вид



Проводя сокращения в числителе и знаменателе и предполагая, что то соответствующий x(t) является одним из узлов интерполяции и значение интерполяционного по) получим, что

В полученной записи интерполяционного полинома коэффициентне зависит ни от интерполируемой функции, ни от набора узлов. Такие коэффициенты можно вычислить один раз перед вычислением интерполяционного полинома и записать в массив, так же как и значения [4, 7]

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

Необходимо найти интерполяционный полином Ньютона.

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



Рис. 1.12. График найденной функции полином Ньютона

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

Метод наименьших квадратов

Метод наименьших квадратов (МНК) -- один из базовых методов регрессионного анализа для оценки неизвестных параметров регрессионных моделей по выборочным данным. Метод основан на минимизации суммы квадратов остатков регрессии.

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

Если набор экспериментальных данных получен со значительной погрешностью, то не имеет смысла использовать интерполяцию Лагранжа полиномами и сплайнами для обработки результатов. В этом случае необходимо провести аппроксимирующую кривую, которая не проходит через экспериментальные точки, но в то же время отражает исследуемую зависимость и сглаживает "выбросы", возможные за счет погрешности эксперимента [6, 8].

Обозначим узлы исходной таблицы данных через х_i, где 0<= i<= n - номер узла. Считаем известными значения экспериментальных данных в узловых точках f(x_i) = f_i. Введем непрерывную функцию ц(х) для аппроксимации дискретной зависимости f(x_i). В узлах функции ц(х) и f(x) будут отличаться на величину е_i = ц(х_i)-f(x_i). Отклонения еi могут принимать положительные и отрицательные значения. Чтобы не учитывать знаки, возведем каждое отклонение в квадрат и просуммируем квадраты отклонений по всем узлам:

Метод построения аппроксимирующей функции ц(х) из условия минимума величины Q называется методом наименьших квадратов (МНК).

Наиболее распространен способ выбора функции ц(х) в виде линейной комбинации:

ц(х) = c₀ ц₀(х) + c₁ ц₁(х) + … +c_m ц₀(х)

ц₀(х), ц₁(х),…, ц_m(х) - базисные функции; m <= n; c₀, c₁,...,c_m - коэффициенты, определяемые при минимизации величины Q.

Математически условия минимума суммы квадратов отклонений Q запишем, приравнивая нулю частные производные от Q по коэффициентам с_k, 0<=k<=n:

Из системы линейных алгебраических уравнений (3) определяются все коэффициенты с_k. Система (3) называется системой нормальных уравнений. Матрица этой системы имеет следующий вид:

и называется матрицей Грама. Элементы матрицы Грама являются скалярными произведениями базисных функций [5,8]

Расширенная матрица системы уравнений получится добавлением справа к матрице Грама столбца свободных членов

где скалярные произведения, являющиеся элементами столбца, определяются аналогично:

Отметим основные свойства матрицы Грама, полезные при программной реализации алгоритмов МНК:

1) матрица симметрична, т.е. а_i,j = а_j,i, что позволяет сократить объем вычислений при заполнении матрицы;

2) матрица является положительно определенной, следовательно, при решении системы нормальных уравнений методом исключения Гаусса можно отказаться от процедуры выбора главного элемента;

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

При обработке экспериментальных данных, определенных с погрешностью e в каждой узловой точке, обычно начинают с аппроксимации функцией ц(x), представимой одной-двумя базисными функциями. После определения коэффициентов ck вычисляют величину Q по формуле (1). Если получится, что корень(Q) > e, то необходимо расширить базис добавлением новых функций ц_k(х). Расширение базиса необходимо осуществлять до тех пор, пока не выполнится условие корень(Q)~=е [5, 8].

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

Выберем базисные функции ц_k(х) в виде последовательности степеней аргумента х, которые линейно независимы,

В этом случае мы будем аппроксимировать экспериментальную зависимость полиномом. Степень полинома m выбираем m<<n (при лагранжевой интерполяции m=n). Аппроксимирующая кривая в МНК не проходит через значения исходной функции в узлах, но проведена из условия наименьшего суммарного квадратичного отклонения. Экспериментальные данные "сглаживаются" с помощью функции ц(х). Если же выбрать m=n, то на основании единственности интерполяционного полинома получим функцию ц(х), совпадающую с каноническим интерполяционным полиномом степени n, аппроксимирующая кривая пройдет через все экспериментальные точки и величина Q будет равна нулю. Последнее обстоятельство используется для отладки и тестирования программ, реализующих алгоритмы МНК [6, 7].

Запишем расширенную матрицу системы нормальных уравнений для базиса (8):

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

Для решения систем уравнений с матрицей Грама разработаны методы сингулярного разложения. Если же m <= 4..5, то такие системы можно решать и более простым методом исключения Гаусса.

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

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

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

(1.5)

Составим вспомогательную таблицу

Таблица 1.1

n	x	y	x2	x3	x4	x5	x6
0	0,0000	0,8200	0,0000	0,0000	0,0000	0,0000	0,0000
1	0,1250	0,7408	0,0156	0,0020	0,0002	0,0000	0,0000
2	0,2500	0,6802	0,0625	0,0156	0,0039	0,0010	0,0002
3	0,3750	0,6410	0,1406	0,0527	0,0198	0,0074	0,0028
4	0,5000	0,6263	0,2500	0,1250	0,0625	0,0313	0,0156
5	0,6250	0,6388	0,3906	0,2441	0,1526	0,0954	0,0596
6	0,7500	0,6817	0,5625	0,4219	0,3164	0,2373	0,1780
7	0,8750	0,7578	0,7656	0,6699	0,5862	0,5129	0,4488
?	3,5000	5,5866	2,1875	1,5313	1,1416	0,8853	0,7050

Подставим данные из таблицы 1.1 в систему 1.5, и получим систему уравнений 1.6, решив систему уравнений 1.6, получим выражение с помощью которого можно найти коэффициенты аппроксимирующую структуры.



Вычислим коэффициенты аппроксимирующего полинома по выражению (1.7)

A₀ = 0.8211;

A₁ = -0.7204;

A₂ = 0.5604;

A₃ = 0.2048.

Поставив найденные коэффициенты аппроксимирующего полинома в формулу (1.8) получим уравнение кубической зависимости данной функции:

Y =



Рис. 1.13. График найденного полинома

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

1.3 Исследование спектральных методов цифровой обработки сигналов

Наряду представлением сигналов во временной области широко применяется представление сигналов в спектральной области.

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

К задачам, решаемым методами спектрального анализа, относятся, обнаружение, разрешение и оценивание параметров сигналов, сжатие данных, выделение информативных признаков, идентификация объектов (определение частотных, импульсных и других характеристик), распознавание образов (речи, изображений). Для случайных сигналов с помощью спектрального анализа решается общая задача выявления скрытых периодичностей и статистических (корреляционных) связей. Спектральный анализ детерминированных периодических (регулярных) сигналов и сигналов конечной длительности называют также гармоническим анализом [9, 11].

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

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

Прямое преобразование ДПФ осуществляется над конечным числом отсчетов (2 n ) входного сигнала, что позволяет преобразовать сигнал из временного представления в спектральное [13, 14].

Обратное ДПФ преобразование позволяет наоборот перевести сигнал из представления в спектральной области в его представление во временной области в том же количестве точек.

Спектральные коэффициенты F(k) получаются в результате свертки (попарного перемножения) входных отсчетов сигнала с базисными функциями Фурье - синусами и косинусами.

Так как базисные функции Sin и Cos равны нулю соответственно в точках 0 и 90 градусов фазовой плоскости, поэтому в разложении присутствуют как вещественная (Cos) так и мнимая (Sin) составляющие.

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

Где, a_k- Спектральные коэффициенты базиса Адамара, N - количество отсчетов, W(n) - матрица базиса Адамара, x(n) - дискретные отсчеты сигнала.

Формула обратного преобразования:

Количество спектральных коэффициентов меньше числа отсчетов сигнала из-за эффекта сходимости, многие коэффициенты равны или близки к нулю. В отличие от базиса Фурье здесь нет умножения на значения Sin или Cos, так как базисные функции имеют вид +1 или -1, т.е. отсчеты сигнала умножаются на единицу, фактически выполняется операция присвоения знака. [14, 15]

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

Требуется найти спектр сигнала на основе базиса функции Адамара.

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

Построим спектрограмму исходного сигнала:



Рисунок 1.14

Сделаем обратное преобразование и найдём исходный сигнал по формуле:

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

Постановка задачи.

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

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

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

Глава 2. Разработка спектральных методов согласование частоты дискретизации

2.1 Построение алгоритма полиномиального представление сигналов

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

Рис. 2.1. Блок-схема перевода сигнала в полиномиальную форму с помощью метода наименьших квадратов

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

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



где, m- степень полинома

С помощью полученного выражение (2.1) можно восстановить сигнал с нужным количеством отсчетов в соответствующем фрагменте сигнала.

В качестве основы для перевода сигнала из его обычного представления в виде последовательности отсчетов f(t) в спектр сигнала берётся одна из базисных функции преобразования Фурье. При этом блок-схема перевода сигнала в полиномиальную форму будит иметь структуру показанной на рис. 2.2.

Рис. 2.2. Блок-схема перевода сигнала в полиномиальную форму с помощью спектрального подхода.

Для представления фрагментов звуковых сигналов (при кусочно-полиномиальной аппроксимации) на практике достаточно выбрать m=1,2,3, что представляет выходные сигналы полиномами не выше третьей степени.

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

Сегодня спектральный анализ широко используется в разных областях, таких как радио и телевещание, мобильный связи, при обработке и передачи цифровых сигналов, изображений, а также видео изображений. Разработаны, многие базисной системы Фурье преобразований, таких как дискретное Фурье преобразование, дискретное косинусное преобразование, Вейвлет преобразование, преобразование Уоша-Адамара, преобразование Хаара, для каждого преобразование разработаны быстрые алгоритмы прямого преобразование. От поставленной задачи выберается та или иная система базиса. Существует множество других преобразований. Преобразование Виногpада Фурье и алгоритм первоначального множителя представляют собой оригинальные, но слишком сложные методы повышения скорости вычисления БПФ. Дискретное косинус преобразование используются для ускорения вычисления корреляции и свертки, а также в спектральном анализе, эти методы особенно подходит применяются еще и для сжатия данных при, например, передаче речи или видеосигналов, а также для записи медицинских сигналов, таких как сигналы ЭКГ или ЭЭГ. А также, они используются при распознавании шаблонов. При преобразовании Уолша сигнал раскладывается на прямоугольные импульсы, а не на синусоиды, и оно вычисляется быстрее, чем БПФ. Преобразование Адамара, построенное с помощью перестановки последовательности Уолша, вычисляется еще быстрее. Хотя для некоторых целей преобразования Уолша и Адамара позволяют получить определенные преимущества, они имеют ряд недостатков, которые ограничивают область их применения, преобразование Хаара особенно полезно для определения краев при обработке изображений и в подобных приложениях [13].

Система Фурье хорошо подходит для обработки гармонических сигналов колебательного характера. Для сигналов с другими свойствами -негармонических, пологих, имеющих характер трендов, с нарушениями стационарности, а также графических изображений (графики, фото) применяются дискретные преобразования с другими типами и формами базисных функций. К такому типу относятся базисы Адамара, Хаара, вейвлет - функции и многие другие аналогичные базисные системы. [14]

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

В качестве базисной системы выберем систему функций Уолша-Адамара. Преобразование Адамара не требует больших вычислительных ресурсов. В основе лежит понятие матрицы Адамара. [15]

В отличие от базиса Фурье здесь нет умножения на значения Sin или Cos, так как базисные функции имеют вид +1 или -1, т.е. отсчеты сигнала умножаются на единицу, фактически выполняется операция присвоения знака.

Матрица ядра этого преобразование содержит целочисленные коэффициенты из множества {-1; +1}. Очевидно, что при выполнении подобных преобразований существенно сокращается объем вычислений за счет исключения умножения в каждой базовой операции. Матрица ядра преобразования Уолша - Адамара для N =2^m может быть описана как результат кронекеровского произведения m матриц дискретных экспоненциальных функций (ДЭФ) E₂ размера 2х2:

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



аналогично можно получить и для N = 8:

Матрица ядра Адамара обладает следующими свойствами:

1) цикличностью a_N+k = a_k; a_N-k = a_-k

2) мультипликативностью a_k+l = a_k * a_l

3) симметричностью A_N = A_N^T

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

где c_k-спектральные коэффициенты, f(n)- отсчеты сигнала, H(n) матрица базисной системы.

...