Загадки быстрого преобразования Фурье
• Метод фазово-амплитудной интерполяции (ФАИ)
• Точное определение частоты, амплитуды и фазы гармоник сигнала
Алгоритм быстрого преобразования Фурье (БПФ) - важный инструмент для анализа и обработки сигналов различной природы.
Он позволяет реконструировать амплитудный и фазовый спектры сигнала в частотной области представления по его амплитудным отсчётам во временной, при этом метод вычислительно оптимизированный при скромном расходе памяти.
Хотя в процессе преобразования никакая информация о сигнале не утрачивается (вычисления обратимы до округлений) алгоритму присущи некоторые особенности:
разрешающая способность, своего рода, "дискретная пикселизация" наподобие гистограммы
эффект спектральной утечки, обуславливающий характерное растекание информации по столбикам гистограммы для тех гармоник сигнала, частота которых не попадает весьма точно на деления дискретной частотной решётки
Эти факторы затрудняют высокоточный анализ и тонкую обработку результатов в дальнейшем.
В статье представлен действенный способ преодоления таких "неудобных" особенностей алгоритма.
Материалы снабжены иллюстрациями, однако для более глубокого погружения в тему рекомендуется скачать десктоп-приложение Solfeggio для высокоточной обработки и визуализации сигналов. С его помощью можно в деталях рассмотреть сигнал и "прощупать" его спектры на различных уровнях масштаба: от общей картинки вплоть до отдельных дискретных отсчётов.
ВАЖНО! Программа и алгоритмические наработки бесплатны для учебного и некоммерческого использования, однако есть ряд важных ограничений на области применения. Пожалуйста, перед внедрением представленных наработок в какие-либо проекты внимательно прочтите небольшое соглашение...
Именно в процессе создания Solfeggio были выявлены закономерности, которые в результате привели к разработке алгоритма фазово-амплитудной интерполяции для быстрого преобразования фурье (ФАИ для БПФ). При этом были сформулированы нерешённые задачи.
Теперь кратко пройдёмся по основным этапам и положениям теории обработки сигналов:
приёмник (микрофон или иной датчик) улавливает аналоговый сигнал и преобразует его в эквивалентные колебания напряжения электрического тока
аналогово-цифровой преобразователь (АЦП) с заданной частотой (то есть через равные интервалы времени) измеряет уровень напряжения входного сигнала и формирует массив (кадр) его временнЫх дискретных отсчётов, квантованных по уровню амплитуды
затем на основе сформированного кадра выполняется вычисление прямого БПФ, в результате которого на выходе получается массив спектральных данных, который даёт некоторое представление о гармониках исходного аналогового сигнала (амплитудах и фазах доминирующих частот), что с определённой погрешностью позволяет реконструировать свойства исходного аналогового сигнала (см. рис. 1)
Как упоминалось ранее, поскольку спектральное представление является дискретным, а не аналоговым, ему присуща разрешающая способность, своего рода, пикселизация наподобие столбиков гистограммы, а также может наблюдаться спектральная утечка, которая обуславливает характерное растекание информации о гармониках сигнала, частота которых не попадает довольно точно в деления дискретной частотной решётки. Эти явления на практике затрудняют последующую реконструкцию и анализ исходных сигналов.
Важно отметить, что применение обратного БПФ к полученному массиву спектральных данных позволяет однозначно восстановить исходный дискретный сигнал во временной области (с точностью до округлений), что указывает на то, что в процессе преобразования никакая информация из сигнала не утрачивается. Это означает, что вычисления обратимы.
В свою очередь согласно теореме Найквиста-Шеннона-Котельникова аналоговый сигнал может быть однозначно и без искажений восстановлен по своим дискретным временным отсчётам, когда частота дискретизации как минимум в два раза превышает частоту наиболее высокочастотной гармоники сигнала. Это наводит на мысль, что, интерполировав дискретный сигнал во временной области (обратно превратив в аналоговый) и повысив частоту дискретизации (повторно передискретизировав), можно достичь более высокого разрешения в фазо-частотной области представления за счёт увеличения вычислительной сложности и повышения потребления памяти. Этот подход применим, однако не всегда удобен на практике.
ЗАГАДКА №1
Но раз возможна интерполяция сигнала во временной области представления, может быть, существует какой-то метод интерполяции в частотной? Хотя на момент написания статьи математического обоснования этой гипотезе нет, экспериментальные данные указывают на то, что для малозашумленных сигналов подобный подход существует, и для многих распространённых практических ситуаций даёт он весьма точные и надёжные результаты.
ЗАГАДКА №2
Хотя точность метода высока, однако и ему присуща некоторая погрешность. Возможно, эту погрешность можно снизить, уточнив формулы. Вопрос также остаётся открытым.
Метод интерполяции интуитивен и на удивление вычислительно прост. Он был экспериментально обнаружен во время разработки программы-инструмента Solfeggio для высокоточного спектрального анализа.
Итак, сначала рассмотрим классический алгоритм быстрого преобразования Фурье на языке C#. Ниже представлены две модификации операции "бабочка" (с прореживанием по времени и по частоте), которые различаются порядком и структурой циклов, но в итоге ведут к одному результату
private static void DecimationInFrequency(ref Complex[] sample, bool direct)
{
if (sample.Length < 2) return;
var length = sample.Length / 2; // sample.Length >> 1
var sampleA = new Complex[length];
var sampleB = new Complex[length];
var abs = (Pi.Single / length).InvertSign(direct);
var rotorBase = new Complex(Cos(abs), Sin(abs));
var rotor = Complex.One; // rotor = rotorBase.Pow(0)
for (int i = 0, j = length; i < length; i++, j++)
{
var a = sample[i];
var b = sample[j];
sampleA[i] = a + b;
sampleB[i] = (a - b) * rotor;
rotor *= rotorBase; // rotor = rotorBase.Pow(i + 1)
}
DecimationInFrequency(ref sampleA, direct);
DecimationInFrequency(ref sampleB, direct);
for (int i = 0, j = 0; i < length; i++) // j += 2
{
sample[j++] = sampleA[i];
sample[j++] = sampleB[i];
}
}
private static void DecimationInTime(ref Complex[] sample, bool direct)
{
if (sample.Length < 2) return;
var length = sample.Length / 2; // sample.Length >> 1
var sampleA = new Complex[length];
var sampleB = new Complex[length];
var abs = (Pi.Single / length).InvertSign(direct);
var rotorBase = new Complex(Cos(abs), Sin(abs));
var rotor = Complex.One; // rotor = rotorBase.Pow(0)
for (int i = 0, j = 0; i < length; i++) // j+=2
{
sampleA[i] = sample[j++];
sampleB[i] = sample[j++];
}
DecimationInTime(ref sampleA, direct);
DecimationInTime(ref sampleB, direct);
for (int i = 0, j = length; i < length; i++, j++)
{
var a = sampleA[i];
var b = sampleB[i] * rotor;
sample[i] = a + b;
sample[j] = a - b;
rotor *= rotorBase; // rotor = rotorBase.Pow(i + 1)
}
}
Нормализация значений после операции "бабочка" выполняется следующим образом
public static partial class Butterfly
{
public static Complex[] Transform(this IEnumerable<Complex> sample, bool direct, bool inTime = false)
{
var workSample = sample.ToArray();
if (inTime) DecimationInTime(ref workSample, direct);
else DecimationInFrequency(ref workSample, direct);
double factor = direct ? workSample.Length / 2 : 2;
for (var i = 0; i < workSample.Length; i++) workSample[i] /= factor;
return workSample;
}
}
Метод фазово-амплитудной интерполяции (ФАИ)
Теперь рассмотрим метод фазово-амплитудной интерполяции (ФАИ).
public static IEnumerable<Bin> Interpolate(this IList<Bin> spectrum, List<int> resonances = default)
{
resonances?.Clear();
var count = spectrum.Count / 2 - 4;
for (var i = 0; i < count; i++)
{
/* Frequency (F); Magnitude (M); Phase (P); */
spectrum[i + 0].Deconstruct(out var aF, out var aM, out var aP);
spectrum[i + 1].Deconstruct(out var bF, out var bM, out var bP);
spectrum[i + 2].Deconstruct(out var cF, out var cM, out var cP);
spectrum[i + 3].Deconstruct(out var dF, out var dM, out var dP);
double GetPeakProbabilityByPhase() => Math.Abs(cP - bP) / Pi.Single;
double GetPeakProbabilityByMagnitude()
{
var bcM = (bM + cM) - (aM + dM) / Pi.Half;
return
(aM < bcM && bcM > dM)
&&
(bM * 0.95 < bcM && bcM > cM * 0.95)
? 0.95
: 0.05;
}
var peakProbabilityByPhase = GetPeakProbabilityByPhase();
var peakProbabilityByMagnitude = GetPeakProbabilityByMagnitude();
var peakProbability = peakProbabilityByPhase * peakProbabilityByMagnitude;
if (peakProbabilityByMagnitude > 0.5 && peakProbabilityByPhase < 0.5)
resonances?.Add(i);
if (peakProbability > 0.5)
{
/*
var halfStepF = (cF - bF) / 2;
var bcMiddleF = (bF + cF) / 2;
var bcOffsetScale = (cM - bM) / (cM + bM);
var bcF = bcMiddleF + bcOffsetScale * halfStepF;
*/
var bcF = (bF * bM + cF * cM) / (bM + cM);
var bcM = (bM + cM) - (aM + dM) / Pi.Half;
var bcP = bP + (bcF - bF) * (cP - bP) / (cF - bF);
/* y(x) = y0 + ( x - x0) * (y1 - y0) / (x1 - x0) */
var abF = aF + (bcF - bF);
var abM = aM;
var abP = aP;
var dcF = dF + (bcF - cF);
var dcM = dM;
var dcP = dP;
yield return new(in abF, in abM, in abP);
yield return new(in bcF, in bcM, in bcP);
yield return new(in dcF, in dcM, in dcP);
i += 3;
}
else
{
yield return new(in aF, in aM, in aP);
}
}
}
Точное определение частоты, амплитуды и фазы гармоник сигнала
Ключевыми являются три строчки
var bcF = (bF * bM + cF * cM) / (bM + cM);
var bcM = (bM + cM) - (aM + dM) / Pi.Half;
var bcP = bP + (bcF - bF) * (cP - bP) / (cF - bF);
/* y(x) = y0 + ( x - x0) * (y1 - y0) / (x1 - x0) */
Именно по этим формулам с довольно высокой точностью вычисляется частота, амплитуда и фаза размытой по спектру гармоники. Первая и треться формулы для частоты и фазы выведены интуитивно, вторая для амплитуды подобрана экспериментальным путём. На момент написания статьи строго математического обоснования этим формулам не найдено.
Примечательно, что представленные формулы наилучшим образом работают для прямоугольного окна, то есть не требуют дополнительной предобработки дискретного кадра. Применение каких-либо других оконных функций лишь снижает их точность.
Наличие нескольких близлежащих гармоник в спектре может вносить дополнительные искажения в точность интерполяции, но для отдалённых друг от друга гармоник погрешности снижаются к минимуму.
Из каких наблюдений получены формулы?
Рассмотрим амплитудный и фазовый спектры единичной гармоники при попадании частоты на деление дискретной решётки, а также в разные части интервала между двумя близлежащими делениями.
Частоту i-того деления решётки можно рассчитать по следующей формуле
var binToFrequency = sampleRate / frameSize;
var binFrequency = i * binToFrequency;При попадании гармоники в деление решётки на амплитудном спектре наблюдается чёткий пик, а эффект спектрального размытия отсутствует, значёния амплитуд на остальных делениях близки к нулевым. Значение фазы также определяются с высокой точностью, однако на фазовом спектре начинают наблюдаться резкие фазовые скачки.
При лёгком сдвиге частоты в интервал между двумя делениями решётки на амплитудном спектре начинает проявляться эффект утечки и характерный срез пика у гармоники, который достигает своего максимума на средине интервала. При этом интересная картина наблюдается на фазовом спектре: значение фазы гармоники теперь лежит на прямой, которая соединяет значения на двух близлежащих делениях решётки, а остальная часть фазового спектра образует гладкие кривые.
Внимательное рассмотрение процесса сдвига раскрывает следующую закономерность
var halfStepF = (cF - bF) / 2;
var bcMiddleF = (bF + cF) / 2;
var bcOffsetScale = (cM - bM) / (cM + bM);
var bcF = bcMiddleF + bcOffsetScale * halfStepF;
Выполнение подстановок и преобразований с упрощениями приводит к следующим выражениям для расчёта точной частоты гармоники и её фазы
var bcF = (bF * bM + cF * cM) / (bM + cM);
var bcP = bP + (bcF - bF) * (cP - bP) / (cF - bF);
/* y(x) = y0 + ( x - x0) * (y1 - y0) / (x1 - x0) */Значение амплитуды угадать немного сложнее, но с достаточно высокой точность оно описывается выражением
var bcM = (bM + cM) - (aM + dM) / Pi.Half;Необходимо упомянуть методы выявления размытых пиков в сигнале для их дальнейшей интерполяции. Существует два характерных признака по которым можно определить кандидата: срезанный пик на амплитудном спектре и скачок на фазовом. Но более надёжным является комбинированный вероятностный подход.
double GetPeakProbabilityByPhase() => Math.Abs(cP - bP) / Pi.Single;
double GetPeakProbabilityByMagnitude()
{
var bcM = (bM + cM) - (aM + dM) / Pi.Half;
return
(aM < bcM && bcM > dM)
&&
(bM * 0.95 < bcM && bcM > cM * 0.95)
? 0.95
: 0.05;
}
var peakProbabilityByPhase = GetPeakProbabilityByPhase();
var peakProbabilityByMagnitude = GetPeakProbabilityByMagnitude();
var peakProbability = peakProbabilityByPhase * peakProbabilityByMagnitude;
if (peakProbabilityByMagnitude > 0.5 && peakProbabilityByPhase < 0.5)
resonances?.Add(i);
if (peakProbability > 0.5)Выявление резонансов
Следует отметить, что характерные амплитудные срезы могут также возникать при резонансах, когда две или более значимые гармоники лежат в одном или близлежащих интервалах дискретной решётки, но в таком случае на фазовом спектре начинает возникать эффект слияния двух скачков, что позволяет отличить резонанс от явления спектральной утечки.
В общем случае резонирующими являются две любые гармоники, просто характерный эффект биений начинает становится наиболее заметным, когда их частоты достаточно близки.
ЗАГАДКА №3
Интересное наблюдение: резонанс очень напоминает амплитудную модуляцию низкочастотного сигнала высокочастотным, например, резонанс 440 и 444 герц похож на 4 герцевый сигнал (или частоты такого порядка), модулированной частотой примерно в 442 герца. Данное предположение требует исследования, уточнений и дальнейшей проверки.
На заключительном этапе, после выполнения ФАИ, из спектральных данных легко выделить чистые гармонические составляющие, которые уже имеют не попиксельное, а растровое представление (частота/амплитуда/фаза).
public static IEnumerable<Bin> EnumeratePeaks(this IList<Bin> spectrum, double silenceThreshold = 0.01)
{
var count = spectrum.Count - 3;
for (var i = 0; i < count; i++)
{
spectrum[i + 0].Deconstruct(out var aF, out var aM, out var aP);
spectrum[i + 1].Deconstruct(out var bF, out var bM, out var bP);
spectrum[i + 2].Deconstruct(out var cF, out var cM, out var cP);
if ((aM + cM) * 0.25d < bM && aM < bM && bM > cM && bM > silenceThreshold)
yield return spectrum[i + 1];
}
}Метод позволяет выявить набор наиболее выраженных гармоник исходного сигнала, после чего их легко представить в виде комплексных чисел либо косинусоид для дальнейшей высокоточной обработки с помощью базовых математических операций.
В некоторых случаях это обеспечивает ощутимое информационное сжатие сигналов за счёт потери малозначимых и шумовых гармоник, благодаря чему можно производить дальнейшие вычисления и преобразования более компактно и оптимально, чем при работе с исходными растровыми массивами.
ИТОГИ
Метод фазово-амплитудной интерполяции открывает новые возможности в области обработки и сжатия сигналов, поскольку позволяет из дискретно-растрового представления сигнала выделить с высокой точностью его аналого-векторную информационную составляющую.
ЗЕРКАЛА СТАТЬИ
EN: makeloft, gitlab, gihub, habr
RU: makeloft, gitlab, gihub, habr
ПЕСОЧНИЦА