Быстрое преобразование Фурье

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

МИНИСТЕРСТВО ОБРАЗОВАНИЯ И НАУКИ РОССИЙСКОЙ ФЕДЕРАЦИИ

ФЕДЕРАЛЬНОЕ ГОСУДАРСТВЕННОЕ БЮДЖЕТНОЕ

ОБРАЗОВАТЕЛЬНОЕ УЧРЕЖДЕНИЕ

ВЫСШЕГО ПРОФЕССИОНАЛЬНОГО ОБРАЗОВАНИЯ

«ВОРОНЕЖСКИЙ ГОСУДАРСТВЕННЫЙ ТЕХНИЧЕСКИЙ

УНИВЕРСИТЕТ» (ФГБОУ ВПО «ВГТУ», ВГТУ)

Кафедра высшей математики и физико-математического моделирования

КУРСОВАЯ РАБОТА

По дисциплине «Математика»

Тема: «Быстрое преобразование Фурье»

Выполнил студент

Группы РТ-172

Веселов Е.М.

Проверил: доцент Борщ Н.А.

Воронеж 2018



Содержание

1. Быстрое преобразование Фурье

2. Практическое задание

2.1 Спектральное представление и спектральные характеристики периодического сигнала

2.2 Спектральное представление и спектральные характеристики четной непериодической функции

2.3 Спектральное представление и спектральные характеристики произвольного непериодического сигнала



1. Быстрое преобразование Фурье

Быстрое преобразование Фурье (БПФ) - это алгоритм вычисления преобразования Фурье для дискретного случая. В отличие от простейшего алгоритма, который имеет сложность порядка O(N²), БПФ имеет сложность всего лишь O(Nlog₂N). Алгоритм БПФ был впервые опубликован в 1965 году в статье Кули (Cooly) и Тьюки (Tukey).

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

Если у вас нет времени и желания разбираться с теорией, то можете сразу скопировать текст программы на C++. Здесьнаходится заголовочный файл fft.h и исходник fft.cpp для быстрого преобразования Фурье для числа отсчетов, равного степени двойки. Вызывать надо функцию fft. А здесь находится заголовочный файл и исходник для произвольного (!) числа отсчетов. Он чуть медленнее, но скорость там тоже порядка Nlog₂N. Вызывать надо функцию universal_fft.

Определение 1.

Дана конечная последовательность x₀, x₁, x₂,...,x_N-1 (в общем случае комплексных). Дискретное преобразование Фурье (ДПФ) заключается в поиске другой последовательности X₀, X₁, X₂,...,X_N-1 элементы которой вычисляются по формуле:

(1).

Определение 2.

Дана конечная последовательность X₀, X₁, X₂,...,X_N-1 (в общем случае комплексных). Обратное дискретное преобразование Фурье (ДПФ) заключается в поиске другой последовательности x₀, x₁, x₂,...,x_N-1 элементы которой вычисляются по формуле:

фурье спектральный преобразование дискретный

(2).

Основным свойством этих преобразований (которое доказывается в соответствующих разделах математики) является тот факт, что из последовательности {x} получается (при прямом преобразовании) последовательность {X}, а если потом применить к {X} обратное преобразование, то снова получится исходная последовательность {x}.

Определение 3.

Величина

называется поворачивающим множителем.

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

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

Прямое преобразование Фурье можно выразить через поворачивающие множители. В результате формула (1) примет вид:

(3).

Эти коэффициенты действительно оправдывают свое название. Нарисуем на комплексной плоскости любое комплексное число, в виде вектора, исходящего из начала координат. Представим это комплексное число в показательной форме: re^jц, где r - модуль числа, а ц - аргумент. Модуль соответствует длине вектора, а аргумент - углу поворота:

Теперь возьмем какой-нибудь поворачивающий множитель . Его модуль равен единице, а фаза - 2р/N. Как известно, при умножении комплексных чисел, представленных в показательной форме, их модули перемножаются, а аргументы суммируются. Тогда умножение исходного числа на поворачивающий множитель не изменит длину вектора, но изменит его угол. То есть, произойдет поворот вектора на угол 2р/N (см. предыдущий рисунок).

Если теперь посмотреть на формулу (3), то станет ясен геометрический смысл преобразования Фурье: он состоит в том, чтобы представить N комплексных чисел-векторов из набора {x}, каждое в виде суммы векторов из набора {X}, повернутых на углы, кратные 2р/N.

Теорема 0.

Если комплексное число представлено в виде e ^j2рN, где N - целое, то это число e ^j2рN = 1.

Доказательство:

По формуле Эйлера, и ввиду периодичности синуса и косинуса:

e ^j2рN = cos(2рN) + j sin(2рN) = cos 0 + j sin 0 = 1 + j0 = 1

Теорема 1.

Величина периодична по k и по n с периодом N. То есть, для любых целых l и m выполняется равенство:

(4).

Доказательство:

(5)

Величина -h = -(nl+mk+mlN) - целая, так как все множители целые, и все слагаемые целые. Значит, мы можем применить Теорему 0:



Что и требовалось доказать по (4).

Теорема 2.

Для величины справедлива формула:

Доказательство:

Алгоритм быстрого преобразования Фурье (БПФ) - это оптимизированный по скорости способ вычисления ДПФ. Основная идея заключается в двух пунктах.

1. Необходимо разделить сумму (1) из N слагаемых на две суммы по N/2 слагаемых, и вычислить их по отдельности. Для вычисления каждой из подсумм, надо их тоже разделить на две и т.д.

2. Необходимо повторно использовать уже вычисленные слагаемые.

Применяют либо "прореживание по времени" (когда в первую сумму попадают слагаемые с четными номерами, а во вторую - с нечетными), либо "прореживание по частоте" (когда в первую сумму попадают первые N/2 слагаемых, а во вторую - остальные). Оба варианта равноценны. В силу специфики алгоритма приходится применять только N, являющиеся степенями 2. Рассмотрим случай прореживания по времени.

Теорема 3.

Определим еще две последовательности: {x_[even]} и {x_[odd]} через последовательность {x} следующим образом:

x_[even]n = x_2n,

x_[odd]n = x_2n+1, (6)

n = 0, 1,..., N/2-1

Пусть к этим последовательностям применены ДПФ и получены результаты в виде двух новых последовательностей {X_[even]} и {X_[odd]} по N/2 элементов в каждой.

Утверждается, что элементы последовательности {X} можно выразить через элементы последовательностей {X_[even]} и {X_[odd]} по формуле:

(7).

Доказательство:

Начинаем от формулы (2), в которую подставляем равенства из (6):

(8)

Теперь обратим внимание на то, что:

(9)

Подставляя (9) в (8) получаем:

(10)

Сравним с формулами для X_[even]k и X_[odd]k при k = 0,1,…,N/2-1:

(11)

Подставляя (11) в (10) получим первую часть формулы (7):

Это будет верно при k = 0,1,…,N/2-1.

Согласно теореме 1:

(12)

Подставим (12) в (10), и заменим по теореме 2:



(13)

Для k = N/2,…,N-1 по формуле (2):

(14)

Подставляем (14) в (13):

Эта формула верна для k = N/2,…,N-1 и, соответственно, (k - N/2) = 0,1,…,N/2-1 и представляет собой вторую и последнюю часть утверждения (7), которое надо было доказать.

Формула (7) позволяет сократить число умножений вдвое (не считая умножений при вычислении X_[even]k и X _[odd]k), если вычислять X_k не последовательно от 0 до N - 1, а попарно: X₀ и X_N/2, X₁ и X_N/2+1,..., X_N/2-1 и X_N. Пары образуются по принципу: X_k и X_N/2+k.

Теорема 4.

ДПФ можно вычислить также по формуле:

(15)

Доказательство:

Согласно второй части формулы (7), получим:

Это доказывает второе равенство в утверждении теоремы, а первое уже доказано в теореме 3.

Также по этой теореме видно, что отпадает необходимость хранить вычисленные X_[even]k и X_[odd]k после использования при вычислении очередной пары и одно вычисление можно использовать для вычисления двух элементов последовательности {X}.

На этом шаге будет выполнено N/2 умножений комплексных чисел. Если мы применим ту же схему для вычисления последовательностей {X_[even]} и {X_[odd]}, то каждая из них потребует N/4 умножений, итого еще N/2. Продолжая далее в том же духе log₂N раз, дойдем до сумм, состоящих всего из одного слагаемого, так что общее количество умножений окажется равно (N/2)log₂N, что явно лучше, чем N² умножений по формуле (2).

Рассмотрим БПФ для разных N. Для ясности добавим еще один нижний индекс, который будет указывать общее количество элементов последовательности, к которой этот элемент принадлежит. То есть X_{R}k - это k-й элемент последовательности {X_{R}} из R элементов. X_{R}[even]k - это k-й элемент последовательности {X_{R}[even]} из R элементов, вычисленный по четным элементам последовательности {X_{2R}}. X_{R}[odd]k - это k-й элемент последовательности {X_{R}[odd]}, вычисленный по нечетным элементам последовательности {X_{2R}}.

В вырожденном случае, когда слагаемое всего одно (N = 1) формула (1) упрощается до:



,

Поскольку в данном случае k может быть равно только 0, то X_{1}0 = x_{1}0, то есть, ДПФ над одним числом дает это же самое число.

Для N = 2 по теореме 4 получим:

Для N = 4 по теореме 4 получим:

Отсюда видно, что если элементы исходной последовательности были действительными, то при увеличении N элементы ДПФ становятся комплексными.

Для N = 8 по теореме 4 получим:



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

Тогда формулы для N=4 будут использовать те же W-множители, что и формулы для N=8:

Подведем итог:

В основе алгоритма БПФ лежат следующие формулы:



2. Практические задания

2.1 Спектральное представление и спектральные характеристики периодического сигнала

2.1.1. Дана кусочно-линейная функция

(1)

являющаяся математической моделью некоторого сигнала.

2.1.2 Продолжим эту функцию периодически на всю числовую ось, получим периодическую функцию f(x) c периодом Т = 2, график которой изображен на рис.1.

Рис. 1 График периодически продолженной функции

Эта функция удовлетворяет условиям Дирихле:

1) f(x) непрерывна на главном периоде [-1;1], за исключением конечного числа точек разрыва первого рода;

2) производная имеет на отрезке [-1;1] конечное число точек разрыва первого рода.

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

[-1;1] кусочно-гладкой. Ряд Фурье кусочно-гладкой функции сходится к значению f(x) в каждой точке непрерывности функции и к значению в каждой точке разрыва.

2.1.3 График суммы ряда Фурье приведен на рис.2. В точках непрерывности функции график суммы ряда полностью совпадает с графиком функции. В точках разрыва первого рода значения суммы ряда отличаются от значений заданной функции и равны среднему арифметическому левого и правого пределов функции в этой точке.

Рис. 2 График суммы ряда Фурье периодической функции

2.1.4. Для функции f(x), периодической с периодом Т=2 и удовлетворяющей на отрезке [-1;1] условиям Дирихле, можно записать ряд Фурье:

(2)

Здесь

(3)

(4)

Соотношение (2) означает, что функция f(x) соответствует ряд Фурье, записанный справа. Согласно теореме Дирихле, равенство левой и правой частей в (2) выполняется в отдельных точках. Из этого, однако, еще не следует равенство функций, т.е. что ряд в правой части (2) сходится, причем именно к f(x). Дейтсвительно, если значения двух функций отличаются только в конечном числе точек, то интегралы (3), (4), определяющие их коэффициенты Фурье, будут одинаковыми. Такие функции имеют один и тот же ряд Фурье.

Вычислим по формулам (3), (4) коэффициенты ряда Фурье заданной функции:

(5)

= (6)

(7)

По формуле (2) запишем искомый ряд Фурье

(8)

Равенство в (8) имеет место во всех точках непрерывности функции f(x).

В точках разрыва значение суммы ряда равно:



(9)

Проверим выполнение условия (9) в точке x=0. Подставляя x=0 в правую часть (10) и учитывая, что cos0=1, sin0=0, получим

=1,625+0,125=1,75.

Ряд Фурье (2) в общем случае можно записать в комплексной форме с комплексными коэффициентами. Воспользуемся формулами Эйлера:

= (12)

Если ввести обозначения

(13)

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



(14)

Запишем теперь в комплексной форме ряд Фурье заданной функции:

где ,

(15)

2.1.5. Запишем частичные суммы ряда (8):

;

(16)

Графики частичных сумм приведены на рис. 3.

Рис. 3 Графики частичных сумм ряда Фурье

Вывод: с ростом n графики частичных сумм в точках непрерывности

приближаются к графику функции f(x). В точках разрыва значения частичных сумм приближается к

2.1.6. Выражение в ряде Фурье называется n-й гармоникой.

Известно, что

(17)

где (18)

или с учетом четверти.

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

Так как

(19)

Продолжим эти вычисления для n = 4,5,…,8 и занесем данные в таблицу. Откладывая на графиках вертикальные отрезки соответствующей длины, получим амплитудную и фазовую диаграммы данной функции. Диаграммы приведены на рис. 4 и 5.

n
1	0.101	-0.477	0.488	-4.503
2	0	-0.557	0.557	1.57
3	0.011	-0.159	0.16	-4.642
4	0	-0.279	0.279	1.57
5	0.004053	-0.095	0.096	-4.67
6	0	-0.186	0.186	1.57
7	0.002068	-0.068	0.068	-4.682
8	0	-0.139	0.139	1.57

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

Вывод: с ростом n амплитудный спектр убывает (не монотонно).

Рис. 5 Фазовая диаграмма периодического сигнала

Вывод: с ростом n фазовый спектр

2.1.7. При рассмотрении энергетических характеристик периодического сигнала основной интерес представляет средняя мощность, которая совпадает с мощностью, средней за один период T=2l:

(20)

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

(21)

откуда (22)

Проверим выполнение равенства Парсеваля для заданной функции (1). Вычислим интеграл в левой части (21):

(23)

Вычислим правую часть равенства (21). Для этого подсчитаем квадраты коэффициентов:

, . (24)

Найдем сумму ряда в правой части (21):

(25)

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

Окончательно получим

(26)

Для ряда (8) равенство Парсеваля выполняется.

Согласно (22) и (23), средняя мощность сигнала

(27)

Запишем математический смысл равенства Парсеваля.

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

(28)

Достигает минимума, когда коэффициенты ряда вычисляется по формулам (3), (4).

С математической точки зрения выполнение равенства Парсеваля означает, что ряд Фурье сходится в среднем к функции f(x), т.е. что среднеквадратическая ошибка стремится к нулю и выполняется соотношение



. (29)

Из этого, однако, еще не следует, что ряд Фурье сходится к функции f(x) равномерно, т.е. что при достаточно больших N при всех значениях x их отрезка [-1;1] модуль разности можно сделать сколь угодно малым:

(30)

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

2.1.8. Продолжим функцию периодически на всю числовую ось четным образом. Получим периодическую функцию с периодом T=2l=2, график которой приведен на рис. 6.

Рис. 6 График четной периодической функции и суммы ряда Фурье

Вычислим коэффициенты ряда Фурье этой четной функции. Все коэффициенты =0 в силу нечетности подынтегральных функций .

=2,5

= (31)

Действительная форма ряда Фурье четной функции имеет вид

(32)

Ряд Фурье четной функции содержит только косинусоидальные гармоники.

Запишем комплексную форму ряда Фурье. Для четной функции все следовательно, все коэффициенты комплексного ряда Фурье действительные числа:

(33)

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

(34)

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

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

Рис. 7 График амплитудного спектра четной периодической функции

2.1.9. Продолжим функцию периодически на всю числовую ось нечетным образом. Получим периодическую функцию с периодом Т=2, график которой приведен на рис. 8. Нечетная периодическая функция удовлетворяет условиям Дирихле, следовательно, для нее можно записать ряд Фурье. Вычислим коэффициенты этого ряда:

. (35)



Рис. 8 График нечетной периодической функции

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

(36)

Ряд Фурье нечетной функции содержит только синусоидальные гармоники.

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

(37)

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

(38)

Для нечетной функции принято считать = (так как ) и рассматривать амплитудный спектр коэффициентов , которые откладываются на диаграмме с соответствующим знаком (рис. 9).

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

Рис. 9 График амплитудного спектра нечетной периодической функции

2.1.10. Из формул (8), (32) и (36), дающих разложение в ряд Фурье на отрезке [0;1] одной и той же функции , можно сделать следующий вывод:

Вид ряда Фурье зависит от того, как была периодически продолжена функция на всю числовую ось. При четном продолжении ряд содержит только косинусы, при нечетном - синусы, произвольное продолжение содержит оба вида гармоник. Все три ряда, имея различные коэффициенты, в точках x(0;1) сходятся к одним и тем же значениям. Однако четное продолжение более предпочтительно, т.к. ряд Фурье в этом случае сходится равномерно.



2.2 Спектральное представление и спектральные характеристики четной непериодической функции

2.2.1. Пусть задан прямоугольный импульс:

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

Рис. 10 График четного периодического импульса

Периодическая функция g2(t) удовлетворяет условиям Дирихле и может быть разложена в ряд Фурье. Функция g2(t) четная, следовательно, bn=0 и все коэффициенты комплексного ряда Фурье являются действительными числами:

Найдем эти коэффициенты:

Ряд Фурье в комплексной форме имеет вид:

Спектральная диаграмма действительных коэффициентов Cn приведена на рис. 11. Для четной функции фазовый спектр ?n = 0, коэффициенты Cn, n = 0, ±1, ±2,... откладываются на диаграмме со знаком в виде вертикальных отрезков соответствующей длины. Пунктирной линией показана огибающая дискретного спектра.

Рис. 11 График спектральной плотности четного импульса

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

Запишем спектральную плотность для данного прямоугольного импульса:

Здесь используется функция:

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

Рис. 12 График спектральной плотности четного импульса

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

2.2.3. Рассмотрим спектральную плотность импульса g2(t), полученную из пункта 2.1. с увеличением амплитуды в |a| = 2 раза:

Воспользуемся свойством линейности спектральной плотности: если .

Таким образом, спектральная плотность сигнала g2(t) равна:

Вывод: при увеличении амплитуды импульса в a раз, амплитудный спектр сигнала увеличивается во столько же раз, а фазовый спектр сигнала не изменяется.

2.2.4. Рассмотрим спектральную плотность импульса g3(t), полученного из g2(t) с увеличением длительности в |b| = 3 раза:

Воспользуемся свойством изменения масштаба времени:

если

Так как , , то спектральная плотность более длительного импульса будет равна:

.

График спектральной плотности для сигнала представлен на рис.2.7.

Вывод: при увеличении длительности импульса в |b| раз амплитуда спектральной плотности увеличивается во столько же раз, и во столько же раз график спектральной плотности сжимается по оси частот.

2.2.5. Рассмотрим сигнал (рис. 13), поученный из сигнала сдвигом на секунды по оси времени.

Рис. 13 График импульса, сдвинутого по оси времени

Для определения спектральной плотности воспользуемся свойством сдвига сигнала во времени:

Так как , то новая спектральная плотность будет равна:

.

Данная функция является уже комплекснозначной. Однако

,

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

.

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

2.3 Спектральное представление и спектральные характеристики произвольного непериодического сигнала.

Найдем спектральную плотность сигнала

График этой функции приведен на рис. 14.



Рисунок 14

Функция q(t) является абсолютно интегрируемой:

и имеет только один максимум в точке . Следовательно, спектральная плотность функции q(t) существует и определяется:

Амплитудный спектр сигнала равен модулю спектральной плотности:

График амплитудного спектра сигнала для варианта 2 приведен на рис. 15.

Рисунок 15

Найдем фазовый спектр сигнала, он равен аргументу спектральной плотности:

График фазового спектра сигнала изображен на рис.16.

Рисунок 16

Проверим выполнение равенства Парсеваля для непериодического сигнала q(t):

Сначала вычислим интеграл в левой части этого неравенства:

Затем вычислим интеграл в правой части неравенства Парсеваля:

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

Размещено на Allbest.ru