В предыдущей главе мы рассматривали ситуацию, когда частота
колебания равнялась m/T, где m - целое.
Теперь рассмотрим ситуацию, когда это не так. Положим, что
m = n + q, где n - целое и 0 < q < 1.
Воспользуемся формулой (35). Поскольку первые несколько
условий для нецелого m не выполняются, то
остается последняя, самая сложная формула, помеченная словами
"Для остальных k":
Подставим в эту формулу n + q
вместо m и выполним упрощения,
воспользовавшись формулой (34) и введя обозначение ρ = 2πj.
Итого:
  ,  ρ = 2πj (37)
Теперь построим график функции, чтобы понять, как
она себя ведет. Ниже показана трехмерная поверхность.
По горизонтальной оси отложено k,
по вертикальной |Xk| и по
оси, уходящей вглубь плоскости, отложено q
от 0.01 до 0.99.
Рис. 1.
На рисунке видно два ярко выраженных ребра. Первое из них
всегда приходится на k = n и k = n + 1.
Второе ребро получается в результате зеркального эффекта.
Высота пика наименьшая в окрестности q = 0.5.
А наибольшая в окрестности q = 1 и q = 0
- то есть при целочисленном m.
К сожалению, пик не является единственным ненулевым коэффициентом
Фурье. Рядом с ним есть множество меньших, но не нулевых величин.
Если при целочисленном m можно наблюдать единственную
полоску, то при нецелом m = n + q эта полоска
размазывается:
Рис. 2.
На рисунке приведена практическая ситуация. Это - ДПФ для звука,
содержащегося в обычном WAV-файле. Высота синего штриха цвет соответствует
|Xk|. Исходный сигнал содержал
ноту "ля" второй октавы с частотой 440 гц и фазой в 90 градусов. ДПФ было
выполнено для N = 1000. Однако частота дискретизации
звука в WAV-файле составляла 44100 Гц, так что период дискретизации
был равен T = 1000/44100 секунд и из формулы m/T = 440
получим m = 440*(1000/44100) = 9.97,
то есть, не целое. В результате ярко выраженный пик окружают дополнительные
ненулевые значения.
На следующем рисунке:
Рис. 3.
показана "хорошая" ситуация, когда частота исходного звука составляла 441 Гц,
и m = 441*1000/44100 = 10, то есть целое. Вы видите только один
ненулевой отсчет.
Этот эффект будем называть эффектом размазывания. Вы видите, что он
определяет погрешность, с которой можно найти частоту
исходного колебания. Погрешность равна 1/T. При достаточно
большом отклонении m от целого эффект может быть очень заметен.
Например ниже вы видите ДПФ для сигнала, соответствующего ноте "ля-бемоль":
Рис. 4.
Точнее можно попытаться определить параметры m,
A и φ численными методами.
Для поиска φ следует учесть, что изменение
A не повлияет на комплексную фазу (аргумент)
коэффициентов Xk.
В самом деле, мы можем представить коэффициенты в виде:
Xk = (A/2)Z(m,&phi)k,
где Z(m,&phi)k - комплексное число, не зависящее
от действительного числа A, но зависящее
от m и φ:
Z(m,&phi)k =
Также не зависит от A отношение коэффициентов
Xk/Xl = Z(m,&phi)k/Z(m,&phi)l.
Это значит, что у нас есть две целевые функции, с помощью которых
мы можем найти частоту m/T и фазу φ.
Возьмем Xk, максимальное по модулю. Если соседние
отсчеты Xk-1 и Xk+1
равны нулю, то у нас нет эффекта размазывания и параметры
восстанавливаются так, как описано в предыдущей главе. На самом деле
нам придется сравнивать не с нулем, а с некоторым малым числом, поскольку
некоторая погрешность при вычислении ДПФ неизбежна.
Теперь, когда мы убедились в наличии эффекта размазывания, попробуем найти
m и φ после чего восстановим A
по формуле: A = |2Xk / Z(m,&phi)k|.
Сначала найдем два максимальных
отсчета Xk и Xk+1.
Теперь мы знаем, что искомое m лежит на
интервале (k, k+1).
Для нахождения m и φ нужно численно решить
систему уравнений:
|Z(m,&phi)k/Z(m,&phi)k+1| = |Xk/Xk+1|
Arg(Z(m,&phi)k) = Arg(Xk)
Здесь только две неизвестных - m и φ, неизвестная A исключена.
Систему можно решить итерационно. Для этого сначала фиксируем φ и подбираем m,
которое удовлетворяет первому уравнению системы. Потом - наоборот фиксируем найденное m и
подбираем φ, которое удовлетворяет второму уравнению системы. Потом снова возавращаемся
к первому уравнения - до тех пор, пока оба уравнения не окажутся сбалансированы с достаточной точностью.
Также существует метод вычисления частоты m, основанный на сопоставлении фаз. Стефан Бернси использует его
в своем алгоритме транспонирования.
Данный метод быстрее (не требует нескольких итераций для решения системы уравнений), однако требует двух
преобразований Фурье для одного фрагмента с небольшим сдвигом на время ΔT. Этот
алгоритм предполагает, что спектр не меняется за это короткое время - хотя бы в отношении частот.
Допустим, у нас есть фрагмент длиной T. Мы делаем преобразование Фурье не для всего этого
фрагмента, а для двух перекрывающихся фрагментов [0, T - ΔT] и [ΔT, T].
Затем рассматриваем каждую получившуюся гармонику и ее фазу в первом и втором преобразовании. По величине
сдвига фазы вносится поправка q и вычисляется точная частота (n + q)/T.
Наконец, эффект размазывания можно уменьшить, если применять оконные функции.
Хотя в результате применения
оконной функции, спектр может сильно исказиться, но все-таки форма пиков изменяется, сужается их "подошва".
Вот пример для T = 1 сек, N = 44100, m/T = 440.2 Гц, φ = 0.
Слева - спектр с "размазыванием", полученный обычным ДПФ. Справа - с применением оконной функции Хамминга.
Знаете ли Вы, что релятивистское объяснение феномену CMB (космическому микроволновому излучению) придумал человек выдающейся фантазии Иосиф Шкловский (помните книжку миллионного тиража "Вселенная, жизнь, разум"?). Он выдвинул совершенно абсурдную идею, заключавшуюся в том, что это есть "реликтовое" излучение, оставшееся после "Большого Взрыва", то есть от момента "рождения" Вселенной. Хотя из простой логики следует, что Вселенная есть всё, а значит, у нее нет ни начала, ни конца... Подробнее читайте в FAQ по эфирной физике.