Быстрое преобразование Фурье. Зеркальный эффект в БПФ

Предположим, что исходный сигнал состоял из суммы гармоник. f_s(t) = A_s cos(2πtm_s / T + φ_s). Пусть мы этот сигнал подвергли дискретизации, выполнили над ним прямое и обратное преобразование Фурье. Представили в виде суммы гармоник G_k(t) = A_k cos(2πtk / T + φ_k), как это описано в предыдущей главе. Спрашивается, эти гармоники G_k - те же самые, что и исходные гармоники f_s или нет? Оказывается, нет, не те. Но кое-какую информацию об исходных гармониках все же можно попытаться восстановить.

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

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

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

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

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

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

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

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

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

Однако, формулы (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)

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[.

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

Z_k,n = cos(w_nm + φ) e^-jw_nk =
применяем формулу Эйлера:
= 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_nm + φ - w_nk)} + e^{-j(w_nm + φ + w_nk})] =
и выносим за скобки все, что можно:
= (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₂ являются геометрическими прогрессиями, а формула суммы геометрической прогрессии нам известна:

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

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

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

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

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

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 = 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)

Вывод

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

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

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

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

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

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

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

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

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

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

Знаете ли Вы, что только в 1990-х доплеровские измерения радиотелескопами показали скорость Маринова для CMB (космического микроволнового излучения), которую он открыл в 1974. Естественно, о Маринове никто не хотел вспоминать. Подробнее читайте в FAQ по эфирной физике.