Ломоносова, Институт математики, информационных и космических технологий
 Скачать 304.76 Kb.
|
|
Алгоритмы преобразования Фурье. Применение в цифровой обработке сигналов Северный (Арктический) Федеральный Университет им. М.В. Ломоносова, Институт математики, информационных и космических технологий Тюльпин А.А. alekseytyulpin@gmail.com Архангельск – 2013 Часть I. Математический аппарат Тригонометрический ряд Фурье Тригонометрическим рядом Фурье периодической функции f (x), имеющей период T называется ряд a 0 2 + ∞ X k=1 (a k cos kx T + b k sin kx T ), где a 0 = 2 T T /2 Z −T /2 f (x)dx, a n = 2 T T /2 Z −T /2 f (x) cos(kx)dx, b n = 2 T T /2 Z −T /2 f (x) sin(kx)dx Разложение в комплексный ряд Фурье В более общем виде рядом Фурье элемента гильбертова пространства называется разложение этого элемента по ортогональному базису. Существует множество систем ортогональных функций: Уолша, Лагера, Котельникова и др. В пространстве L 2 [a, b] одной из таких систем является ситема функций e −ikx k ∈ Z . Разложение функции f (x) в ряд Фурье по данной системе будет иметь следующий вид: ∞ X k=−∞ C k e −ikx где C k = 1 T T /2 Z −T /2 f (x)e −i2πkx T dx — комплексный интеграл Фурье Преобразование Фурье Преобразованием Фурье будем называть операцию, сопоставляющуюю функции f (x) функцию ˜ f (ω), описывающую коэффициенты разложения f (x) в ряд Фурье: ˜ f (ω) = ∞ Z −∞ f (x)e −iωx dx Обратным преобразованием Фурье будем называть операцию, обращающую преобразование Фурье: f (x) = 1 √ 2π ∞ Z −∞ ˜ f (ω)e iωx dω Часть II. Применение преобразования Фурье в цифровой обработке сигналов Разложение сигнала на гармонические колебания t = 1s, df = 1000Hz, f = 10Hz x(t) = cos(2πf t)+cos(2π(3f )t)+1.2 cos(2π(4.8f )t)+0.3 cos(2π(20f )t) 0.0 0.2 0.4 0.6 0.8 1.0 Time [s] 4 3 2 1 0 1 2 3 4 Amplitude 0 200 400 600 800 1000 Frequency [Hz] 0.0 0.2 0.4 0.6 0.8 1.0 1.2 1.4 Amplitude Свертка функций Свертка – операция в функциональном анализе, показывающая схожесть одной функции с отраженной и сдвинутой копией другой. На рисунке представлен пример сверки двух прямоугольных импульсов: Вычисление свертки Свертка сигналов во временной области: (x ∗ h)(t) = ∞ Z −∞ x(τ )h(t − τ )dτ Пусть X(ω) и H(ω) – образы сигналов x(t) и h(t) соответственно. Тогда преобразование Фурье свертки данных сигналов можно вычислить по формуле: ∞ Z −∞ (x ∗ h)(t)e −iωt dt = X(ω) · H(ω) На основании приведенных выше формул следует, что: (x ∗ h)(t) = 1 √ 2π ∞ Z −∞ X(ω) · H(ω)e iωx dω Цифровая фильтрация 0 200 400 600 800 1000 0.5 0.0 0.5 1.0 1.5 Spectrum of filter 0 10 20 30 40 50 1.5 1.0 0.5 0.0 0.5 1.0 1.5 2.0 Filter operator 0 200 400 600 800 1000 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 Gibbs phenomen Часть III. Алгоритмы вычисления преобразования Фурье. Быстрая свертка Дискретные сигналы Дискретизация сигнала по времени – процедура, состоящая в замене несчетного множества его значений их счетным (дискретным) множеством, которое содержит информацию о значениях непрерывного сигнала в определенные моменты времени. Выбор частоты дискретизации Теорема Котельникова-Найквиста-Шеннона Непрерывный сигнал x(t) с ограниченным спектром можно точно восстановить (интерполировать) по его отсчетам x(n·dT ), взятым через интервалы dT = 1 2F , где F – верхняя частота спектра сигнала. Дискретное преобразование Фурье (ДПФ) Обозначения: x n – значение сигнала в момент времени t n X(ω k ) = X k – значение спектра сигнала в точке ω k N – количество отсчетов ω k = 2πk f k – k-я частота спектра j – мнимая единица X k = N −1 X n=0 x n e −jωkn N x n = 1 N N −1 X k=0 X k e jωkn N Вычислительная сложность ДПФ Пусть {k ∈ 0, N − 1} – отсчеты спектра сигнала x(t). Тогда для вычисления как прямого, так и обратного преобразования Фурье потребуется O(n 2 ) операций: def DFT(x, df): X = zeros(df,dtype = np.complex) N = len(x) for k in xrange(df): for n in xrange(N): X[k] += x[n] * np.exp(- 1 j * 2 * np.pi * k * n / N) return X def IDFT(X): N = len(X) x = zeros(N,dtype = np.complex) for n in xrange(N): for k in xrange(N): x[n] += X[k] * np.exp( 1 j * 2 * np.pi * k * n / N) x[n] /= N return x Введение в быстрое преобразование Фурье (БПФ) Виды БПФ: БПФ с прореживанием по частоте (DIF FFT) БПФ с прореживанием по времени (DIT FFT) Первая программная реализация алгоритма БПФ была осуществлена в начале 60-х годов XX века Джоном Кули в вычислительном центре IBM под руководством Джона Тьюки, а в 1965 году ими же была опубликована статья, посвященная алгоритму быстрого преобразования Фурье (Cooley–Tukey FFT). БПФ. Основная идея Основная идея БПФ состоит в разбиении спектра на 2 части, и вычислении БПФ уже для каждой из частей размерности N/2. БПФ с прореживанием по времени Рассмотрим один из алгоритмов вычиcления преобразования Фурье последовательностей отсчетов длины N = 2 l , l ∈ R за O(N log 2 N ) операций – алгоритм прореживания по времени. Обозначим W nk N = e −jn2πk N . Разобьем последовательность отсчетов сигнала на четные и нечетные: X k = N −1 X n=0 x n · W nk N = N/2−1 X n=0 x 2n · W 2nk N + + N/2−1 X n=0 x 2n+1 · W (2n+1)k N = N/2−1 X n=0 x 2n · W nk N/2 | {z } X 0 (ω k ) + + W k N N/2−1 X n=0 x 2n+1 · W nk N/2 | {z } X 1 (ω k ) = X 0 (ω k ) + W k N X 1 (ω k ) БПФ с прореживанием по времени Рассмотрим разбиение спектра на две части: X = {X k |k = 0, N − 1} = = {X m |m = 0, N/2 − 1} [ {X r |r = N/2, N − 1} = = {X m |m = 0, N/2 − 1} [ {X r+N/2 |r = 0, N/2 − 1} Вычислим вторую половину спектра: X k+N/2 = N/2−1 X n=0 x 2n · W n(k+N/2) N/2 + + W k+N/2 N N/2−1 X n=0 x 2n+1 · W n(k+N/2) N/2 БПФ с прореживанием по времени Рассмотрим некоторые особенности коэффициентов W nk N W n(k+N/2) N/2 = W nN/2 N/2 · W nk N/2 W nN/2 N/2 = e −j2πnN/2 N/2 = e −j2πn = cos(−2πn) − j sin(−2πn) = 1 Получаем: W n(k+N/2) N/2 = W nk N/2 Также, стоит отметить, что W k+N/2 N = e −j2πN 2N W k N = −W k N БПФ с прореживанием по времени Для первой половины спектра: X k = X 0 (ω k ) + W k N X 1 (ω k ) Для второй половины спектра: X k+N/2 = X 0 (ω k ) − W k N X 1 (ω k ) def W(n, N): return np.exp(- 1 j * 2 * np.pi * n/N) def fft(x): N = len(x) if N <= 1 : return x even = fft(x[ 0 :: 2 ]) odd = fft(x[ 1 :: 2 ]) return [even[n] + W(n, N) * odd[n] for n in xrange(N/ 2 )] + \ [even[n] - W(n, N) * odd[n] for n in xrange(N/ 2 )] Литература A. Oppenheim, R. Schafer. Digital Signal Processing L. Franks. Signal theory H. Nussbaumer. Fast Fourier Transform and Convolution Algorithms E. Titchmarsh. Introduction to the Theory of Fourier Integrals E. Chu, A. George. Inside the FFT Black Box: Serial and Parallel Fast Fourier Transform Algorithms А. Н. Колмогоров, С.В. Фомин. Элементы теории функций и функционального анализа Сайт профессора Давыдова – http://prodav.narod.ru/ О. В. Бесов Тригонометрические ряды Фурье Сайт http://dsplib.ru/ Спасибо за внимание! Алгоритмы преобразования Фурье. Применение в цифровой обработке сигналов Северный (Арктический) Федеральный Университет им. М.В. Ломоносова, Институт математики, информационных и космических технологий Тюльпин А.А. alekseytyulpin@gmail.com Архангельск – 2013 |