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

Вернемся к  предыдущей ситуации.

Дана функция f(t) на отрезке [0, T].

Выполнена ее дискретизация, для чего отрезок разбит на N равных частей в точках t_n = Tn/N и вычислены значения функции в этих точках: {x} : x_n = f(t_n) = f(Tn/N).

Пусть выполнено  прямое дискретное преобразование Фурье (далее - ДПФ) {X} : X_k = NA_ke ^jφ_k, и функция разложена на сумму из N гармоник:

G_k(t) = A_k cos(2πtk / T + φ_k)

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

f(t) = A cos(2πtm / T + φ).

Получится ли в  результате ее преобразования последовательность {X}, в которой все элементы равны нулю, кроме элемента X_m = NA_me ^jφ_m, который дает как раз эту гармонику?

G_m(t) = A_m cos(2πtm / T + φ_m) = f(t), A_m = A, φ_m = φ

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

G_m(t) = (A/2) cos(2πtm / T + φ) = f(t) / 2 = f'(t)

и

G_N-m(t) = (A/2) cos(2πt(N - m) / T - φ) = f''(t)

Как видите у  них половинные амплитуды, противоположные  фазы, а частоты зеркально симметрично  расположены на отрезке [0, N]. Это - тот самый зеркальный эффект.

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

Преобразуем сумму  этих гармоник по формуле суммы косинусов:

Итого:

f'(t) + f''(t) = A cos(πtN / T) cos(2πtm / T - πtN / T + φ)    (30)

А нам требовалось:

f(t) = A cos(2πtm / T + φ)    (31)

Однако, формулы (30) и (31) дают один и тот же результат в точках t_n = Tn / N. В самом деле, подставим Tn / N вместо t сначала в (30):

f'(t) + f''(t) =  
= A cos(πTnN / TN) cos(2πTnm / TN - πTnN / TN + φ) =  
= A cos(πn) cos(2πnm / N - πn + φ) = ...

Второй множитель  разложим по формуле косинуса разности, отделив πn:

... = A cos(πn) [cos(2πnm / N + φ) cos(πn) +  
+ sin(2πnm / N + φ) sin(πn)] = ...

Учитывая, что  для целого n выполняется sin(πn) = 0 и cos²(πn) = 1, получаем:

... = A cos(πn) [cos(2πnm / N + φ) cos(πn)] =  
= A cos²(πn) cos(2πnm / N + φ) = A cos(2πnm / N + φ)     (32)

Теперь подставим Tn / N вместо t в (31):

f(t) = A cos(2πtm / T + φ) = A cos(2πTnm / TN + φ) =  
= A cos(2πnm / N + φ)     (33)

Формулы (32) и (33) совпадают, что и требовалось  доказать.

Из этого примера следует важный вывод. Заданная дискретная последовательность {x} может быть разложена в общем случае на разные суммы гармоник G_k(t). Даже в элементарном случае, когда исходная функция представляла собой одну гармонику, в результате можно получить две. То есть, разложение дискретной последовательности на гармоники неоднозначно.

Этим эффектом мы обязаны именно дискретизации. Дело в том, что если вместо ДПФ использовать его непрерывный аналог - разложение в ряд Фурье непрерывной функции  или непрерывное преобразование Фурье f(t), то мы получим единственую правильную гармонику G_m(t) = A cos(2πtm / T + φ) = f(t). Если же мы применяем ДПФ, то получим сумму гармоник, которая только в точках дискретизации совпадает с исходной функцией:

На этом графике для N = 8 и m = 2 синим цветом показана исходная гармоника f(t) и две гармоники, которые получаются в результате преобразвания Фурье: f'(t) зеленым цветом и f''(t) красным. В точках дискретизации, отмеченных вертикальными штрихами, сумма гармоник f'(t) и f''(t) совпадает с гармоникой f(t).

Заметим также, что тот же результат преобразования получился бы, если бы мы в качестве исходной функции f(t) взяли 2f''(t) или f'(t) + f''(t). Это следует из того, что  в результате дискретизации была бы получена та же последовательность {x} и результаты ДПФ, естественно, дали бы то же самое.

Итак, мы имеем  правило:

Разложение на гармоники, когда исходные данные представлены дискретным набором точек {x} является принципиально неоднозначным. Функции  
f(t) = A cos(2πtm / T + φ),  
2f''(t) = A cos(2πt(N-m) / T - φ) и  
f'(t) + f''(t) = (A/2) cos(2πtm / T + φ) + (A/2) cos(2πt(N-m) / T - φ)  
дают после дискретизации одни и те же исходные данные и те же результаты ДПФ.

Доказательство  зеркального эффекта.

В начале главы  упоминалось о том, что в результате ДПФ гармонической функции на практике получаются две гармоники. Однако этот эмпирический факт не доказывался. Докажем теперь строго, какие гармоники дает произвольная гармоническая функция f(t) = A cos(2πtm / T + φ) при целочисленном m [0, N[.

Напомним формулу  прямого ДПФ:

В данном случае

x_n = f(t_n) = f(Tn / N) = A cos(2πTnm / NT + φ) = A cos(2πnm / N + φ)

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

w_n = 2πn / N

Z_k,n = (f(t_n) / A) e^{-j2πkn / N} = cos(w_nm + φ) e^-jw_n^k

В результате формула  прямого ДПФ упрощается до:

X_k = A Z_k,n

Теперь преобразуем Z_k,n:

Z_k,n = cos(w_nm + φ) e^-jw_n^k =  
применяем формулу Эйлера: 
= cos(w_nm + φ) (cos(-w_nk) + j sin(-w_nk)) =  
= cos(w_nm + φ) (cos(w_nk) - j sin(w_nk)) =  
раскрываем скобки: 
= cos(w_nm + φ) cos(w_nk) - j cos(w_nm + φ) sin(w_nk) =  
применяем формулы произведения косинусов и синуса на косинус: 
= (1/2)[cos((w_nm + φ) - w_nk) + cos((w_nm + φ) + w_nk)] -  
= - (j/2)[sin(w_nk - (w_nm + φ)) + sin(w_nk + (w_nm + φ))] =  
перегруппировываем слагаемые: 
= (1/2)[cos((w_nm + φ) - w_nk) + cos((w_nm + φ) + w_nk) -  
- j sin(w_nk - (w_nm + φ)) - j sin(w_nk + (w_nm + φ))] =  
= (1/2)[cos((w_nm + φ) - w_nk) + cos((w_nm + φ) + w_nk) +  
+ j sin((w_nm + φ) - w_nk) - j sin((w_nm + φ) + w_nk)] =  
= (1/2)[cos(w_nm + φ - w_nk) + j sin(w_nm + φ - w_nk) +  
+ cos(w_nm + φ + w_nk) - j sin(w_nm + φ + w_nk)] =  
применяем формулу Эйлера (только наоборот): 
= (1/2)[e ^j(w_n^{m + φ - w}_n^k) + e ^-j(w_n^{m + φ + w}_n^k)] =  
и выносим за скобки все, что можно: 
= (1/2)[e ^jφe ^jw_n^{(m - k)} + e ^-jφe ^-jw_n^{(m + k)}]

Теперь подставляем  полученные величины в сумму ДПФ и преобразуем:

X_k = A Z_k,n =  
= A (1/2)[e ^jφe ^jw_n^{(m - k)} + e ^-jφe ^-jw_n^{(m + k)}] =  
= (A/2)e ^jφ e ^jw_n^{(m - k)} + (A/2)e ^-jφ e ^-jw_n^{(m + k)} =  
подставляем w_n: 
= (A/2)e ^jφ e ^{j2πn(m - k) / N} + (A/2)e ^-jφ e ^{-j2πn(m + k) / N} =  
вводим обозначения для сумм: 
= (A/2)e ^jφS₁ + (A/2)e ^-jφS₂

Легко видеть, что  суммы S₁ и S₂ являются геометрическими прогрессиями, а формула суммы геометрической прогрессии нам известна:

S_N = a₀(q^N - 1) / (q - 1), q ≠ 1    (34)

Первый элемент  прогрессии в обоих случаях равен a₀ = 1.

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

q₁ = e ^{j2π(m - k) / N} для S₁

и

q₂ = e ^{-j2π(m + k) / N} для S₂.

Условие q ≠ 1 вынуждает  нас решить уравнения:

e ^{j2π(m - k) / N} = 1,

и

e ^{-j2π(m + k) / N} = 1,

Учитывая Теорему 0, получим, что условие q ≠ 1 не выполняется  при k = m для S₁ и при k = (N - m) для S₂.

В случае, когда  выполняются оба условия: k = m и k = (N - m), то есть k = m = N /2 обе суммы нельзя считать по формуле геометрической прогресии.

В случае k = m для S₁ придется выполнить небольшие дополнительные преобразования:

S₁ = e ^{j2πn(m - m) / N} = 1 = N

Аналогично в  случае k = N - m для S₂:

S₂ = e ^{-j2πn(m + N - m) / N} = e ^-j2πn = 1 = N

В случае k = m = N / 2 имеем:

X_k = (A/2)Ne ^jφ + (A/2)Ne ^-jφ =  
= (A/2)N(e ^jφ + e ^-jφ) =  
= (A/2)N(cos φ + j sin φ + cos φ - j sin φ) =  
= (A/2)N(2 cos φ) =  
= ANcos φ

В случае k = m = 0 имеем:

X_k = (A/2)e ^jφN + (A/2)e ^-jφN = (A/2)N(e ^jφ + e ^-jφ) =

= (A/2)N(cos φ + jsin φ + cos φ - jsin φ) = ANcos φ

Наконец, получаем формулу для X_k:

Для k = m = N / 2 или k = m = 0: 
X_k = ANcos φ  
Для k = m ≠ N / 2: 
X_k = (A/2)Ne ^jφ + (A/2)e ^-jφ(e ^-j4πm - 1) / (e ^{-j4πm / N} - 1)  
Для k = (N - m) ≠ N / 2: 
X_k = (A/2)e ^jφ(e ^j4πm - 1) / (e ^{j4πm / N} - 1) + (A/2)Ne ^-jφ  
Для остальных k: 
X_k = (A/2)e ^jφ(e ^{j2π(m - k)} - 1) / (e ^{j2π(m - k) / N} - 1) +  
+ (A/2)e ^-jφ(e ^{-j2π(m + k)} - 1) / (e ^{-j2π(m + k) / N} - 1)    (35)

Заметим, что  эта формула получена без использования  факта целочисленности m и k.

Теперь учтем  целочисленность. Для этого применим Теорему 0 и заменим в формуле (35) экспоненты на 1 везде, где выполняется это условие:

Для k = m = N / 2 или k = m = 0: 
X_k = ANcos φ  
Для k = m ≠ N / 2: 
X_k = (A/2)Ne ^jφ + (A/2)e ^-jφ(1 - 1) / (e ^{-j4πm / N} - 1)  
Для k = (N - m) ≠ N / 2: 
X_k = (A/2)e ^jφ(1 - 1) / (e ^{j4πm / N} - 1) + (A/2)Ne ^-jφ  
Для остальных k: 
X_k = (A/2)e ^jφ(1 - 1) / (e ^{j2π(m - k) / N} - 1) +  
+ (A/2)e ^-jφ(1 - 1) / (e ^{-j2π(m + k) / N} - 1)

Сокращаем везде, где получаются нули, и приходим к формулам:  
 
Для k = m = N / 2 или k = m = 0: 
X_k = ANcos φ  
Для k = m ≠ N / 2: 
X_k = (A/2)Ne ^jφ  
Для k = (N - m) ≠ N / 2: 
X_k = (A/2)Ne ^-jφ  
Для остальных k: 
X_k = 0    (36)

Вывод:

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

f(t) = A cos(2πtm / T + φ),

2f''(t) = A cos(2πt(N-m) / T - φ) и

f'(t) + f''(t) = (A/2) cos(2πtm / T + φ) + (A/2) cos(2πt(N-m) / T - φ)

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

f'(t) + f''(t).

При этом все коэффициенты ДПФ равны нулю за исключением 

X_m = (A/2)Ne ^jφ

и

X_{N - m} = (A/2)Ne ^-jφ

кроме частных случаев m = N / 2 и m = 0, в которых единственный ненулевой коэффициент:

X_m = ANcos φ

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

Исправление зеркального эффекта.

Таким образом  зеркальный эффект в подавляющем  большинстве случаев искажает исходную картину, поскольку в действительности очень редко на вход подается сумма  двух гармонических сигналов f'(t) + f''(t) именно с таким соотношением частот: mν и (N - m)ν. В результате исходный спектр искажается, словно отражаясь в зеркале:

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

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

X_m = (A/2)Ne ^jφ  
f(t) = A cos(2πtm / T + φ).

Частота восстанавливается  проще всего: ν = m/T, где m - индекс найденного ненулевого элемента X_m. Амплитуда и фаза восстанавливаются по формуле (29):

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