Преобразование Фурье в цифровой информации

Декомпозиция сигнала в частотную область
Все началось в 1807 году, когда французский математик и физик Жозеф Фурье представил тригонометрическое разложение рядов (в настоящее время известное как метод рядов Фурье) для решения уравнения теплопроводности в частных производных в металлической пластине. Идея Фурье состояла в том, чтобы разложить сложную периодическую функцию на сумму простейших осциллирующих функций - синусов и косинусов.

Джозеф Фурье

Если функция f(x) периодична с периодом T и интерграфируема (ее интеграл конечен) на интервале [x 0, x 0 + T], то ее можно преобразовать в ряд Дискретное преобразование Фурье

Определение
Ряд Фурье можно назвать прародителем преобразования Фурье, которое в случае цифровых сигналов (дискретное преобразование Фурье) описывается формулой:

Преобразование Фурье является обратимым, и мы можем вернуться к временной области путем вычисления:

В некоторых обозначениях мы можем наблюдать, что деление на N переводится в обратное вычисление - это не нарушает вычисление, если мы не применяем деление на N как в прямом, так и в обратном преобразовании Фурье.
Реализации прямого и обратного дискретного преобразования Фурье показаны ниже.
public >возвращает<
Преобразование Фурье сигнала x///>
name="x" param< Signal x samples values
///>
summary summary<
//Complex[] DFT(Double[] x)
{
int N = x.Length; // Количество отсчетов
Complex[] X = новый комплекс[N];
for (int k = 0; k < N; k++)
{
X[k] = 0;
for (int n = 0; n < N; n++)
{
X[k] += x[n] * Complex.Exp(-Комплекс.ImaginaryOne * 2 * Math.PI * (k * n) / Convert.ToDouble(N));
}
X[k] = X[k] / N;
}
возврат X;
}

///

/// Вычисляет обратное дискретное преобразование Фурье заданного спектра X //

///
Комплексные значения спектра

// Выборки сигнала во временной области

public Double[] iDFT (Complex[] X)
{
int N = X.Длина; // Число элементов спектра
Double[] x = новый Double[N];

for (int n = 0; n < N; n++)

{
Комплексная сумма = 0;

for (int k = 0; k < N; k++)

{
sum += X[k] * Complex.Exp(Complex.ImaginaryOne * 2 * Math.PI * (k * n) / Convert.ToDouble(N));

}

x[n] = sum.Real; //В результате мы ожидаем только реальные значения (если наши расчеты верны, мнимые значения должны быть равны или близки к нулю)

}

return x;

}
Ошибка, часто совершаемая в кодовой реализации обратного преобразования Фурье, заключается в переводе комплексного значения в вещественное с помощью значения его величины - в таком случае мы получим | x (n)| вместо x (n).

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

Для проверки алгоритма создадим сигнал, представляющий собой сумму двух синусоидальных волн:

Обозначим через x 1 (n) исходный сигнал и через x 2 (n) сигнал нарушения (шум). Сумма этих сигналов x (n) = x 1 (n) + x 2 (n) в этом случае является нарушенным сигналом.

Рис.1. Исходные и прерывающие сигналы


Рис.2. Сумма сигналов (нарушенный сигнал)

Поскольку мы создали наш сигнал из суммы двух синусоидальных волн, то согласно теореме Фурье мы должны получить его частотное изображение, сосредоточенное вокруг двух частот f1 и f2, а также его противоположностей -f1 и -f2.

Рис.4. Спектр с удаленной частотой прерывания и его обратный результат DFT
Рис.3. Амплитудный спектр нарушенного сингала.

Если теорема Фурье верна, то, удалив из спектра полосы, которые пришли от сигнала нарушения, мы должны получить исходный сигнал.

Рис.4. Спектр с удаленной частотой прерывания и его обратный результат DFT

Как видно из рисунка 4, наш выходной сигнал близок к исходному, а искажения вызваны дополнительными полосами, возникшими в результате численных операций. Давайте уберем их все и оставим только две основные полосы.

Рис.5. Спектр только с основной частотой и его обратный результат DFT.

На рисунке 3 мы не видим отрицательных значений, так почему же мы сказали, что вычисляемый спектр сосредоточен вокруг двух частот f1 и f2, а также его противоположностей - f1 и -f2?

Это потому, что мы используем для анализа сигналов, которые определены в диапазоне, симметричном началу координат. Если мы обозначим наш нарушенный сигнал x (n) как:



Рис.7. Сдвинутый сигнал

и реализовать незначительное изменение в нашем алгоритме:

public >returns<
Преобразование Фурье сигнала x///>
name="SamplesNumbers" param< SamplesNumbers vector
///>name="x"paramvalues
///>
summary summary<
/// Complex[] DFT(Double[] x, int[] SamplesNumbers)
{
int N = x.Length; // Количество отсчетов
Complex[] X = новый комплекс[N];

for (int k = 0; k < N; k++)

{
X[k] = 0;

Double s1 = SamplesNumbers[k]; // Получить индекс выборки

для (int n = 0; n < N; n++)
{
Double s2 = SamplesNumbers[n]; // Получить индекс выборки
X[k] += x[n] * Complex.Exp(-Комплекс.ImaginaryOne * 2 * Math.PI * (s1 * s2) / Convert.ToDouble(N));
}

X[k] = X[k] / N;

}
return X;
}
мы получим амплитудный спектр, показанный ниже.


Рис.8. Амплитудный спектр сдвинутого сигнала

То же самое изменение должно быть реализовано в обратном методе DFT.

public >возвращает<
Выборки сигнала во временной области///>
name="SamplesNumbers" param ///>
name="X" param< Spectrum complex values
///>
summary summary<
/// Double[] iDFT(Complex[] X, int[] SamplesNumbers)
{
int N = X.Длина; // Число элементов спектра
Double[] x = новый Double[N];

for (int n = 0; n < N; n++)

{
Комплексная сумма = 0;
Double s1 = SamplesNumbers[n]; // Получить индекс выборки

для (int k = 0; k < N; k++)

{
Double s2 = SamplesNumbers[k]; // Получить индекс выборки
sum += X[k] * Complex.Exp(Комплекс.ImaginaryOne * 2 * Math.PI * (s1 * s2) / Convert.ToDouble(N));
}

x[n] = sum.Real; // В результате мы ожидаем только действительные значения (если наши вычисления верны, мнимые значения должны быть равны или близки к нулю)

}

return x;

}

Обратите внимание, что, сдвигая сигнал в симметричное положение к оригиналу, мы вычисляем сдвинутый спектр:

X(k)→X(k−N/2)
X(k)→X(k−N/2)

Тот же эффект может быть достигнут, если мы изменим исходный сигнал, используя формулу

x(n)→(-1)n⋅x(n)
x(n)→(−1)n⋅x(n)
В этом случае нет необходимости хранить фактические числа выборок. Для этого используйте метод, показанный ниже.
public >returns<
Shift signal///>
name="x" param< Original signal
///>
summary summary<
/// Double[] Shift(Double[] x)
{
int N = x.Length;

for (int i = 0; i < N; i++)

{
x[i] = (int)Math.Pow(-1, i) * x[i];
}

return x;

}

Дискретное преобразование Фурье с уменьшенным числом вычислений

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

Обозначим вспомогательной переменной W N:

Тогда наша формула принимает вид

Также мы можем заметить, что если k = 0 или n = 0, мы получаем:
Соответственно, обозначение DTF может быть записано в матричном виде:

Отсюда следует, что для определения дикретного преобразования Фурье сигнала N-выборок мы должны выполнить N 2 операций умножения и определить (N-1) 2 коэффициентов в матрице W N. Напомним, что операция умножения является одной из самых трудоемких операций для процессоров.
Вернемся к фактору W N и рассмотрим N=4 отсчетов сигнала x(n). Для вычисления X(k) нам нужно найти значения коэффициентов матрицы:


Это наблюдение приводит нас к выводу, что при определении DFT мы вычисляем только N / 2 - 1 спектральных коэффициентов, а другой мы можем получить из предыдущих - это свойство симметрии спектра. Благодаря этому наблюдению в четыре раза сократилось количество операций внутри алгоритма. Это наблюдение верно только в том случае, если мы имеем дело с физическими сигналами, то есть сигналами с реальными значениями.
public >returns<
Преобразование Фурье сигнала x///>
name="x" param< Signal x samples values
///>
summary summary<
/// Complex[] DFT2(Double[] x)
{
int N = x.Length; // Количество отсчетов
Complex[] X = новый комплекс[N];

// Вычисление элементов с нулевым индексом

X[0] = x.Sum() / N;

for (int k = 1; k < N/2 + 1; k++)

{
X[k] = 0;

for (int n = 0; n < N; n++)

{
X[k] += x[n] * Complex.Exp(-Комплекс.ImaginaryOne * 2 * Math.PI * (k * n) / Convert.ToDouble(N));
}

X[k] = X[k] / N;

// Пропустить дополнительный расчет для центрального элемента
, если (k != N / 2)
{
X[N - k] = новый комплекс(X[k]. Вещественный, -X[k].Мнимое);
}
}

возврат X;

}

///

/// Вычисляет прямое дискретное преобразование Фурье заданного цифрового сигнала x по номерам выборок и приведенному числу умножений ///

///
Значения выборок сигнала x

//
Samples numbers vector

/// Преобразование Фурье сигнала x
public Complex[] DFT2(Double[] x, int[] SamplesNumbers)
{
int N = x.Length; // Количество отсчетов
Complex[] X = новый комплекс[N];

for (int k = 0; k < N / 2 + 1; k++)

{
X[k] = 0;
Double s1 = SamplesNumbers[k];
for (int n = 0; n < N; n++)
{
Double s2 = SamplesNumbers[n];
X[k] += x[n] * Complex.Exp(-Комплекс.ImaginaryOne * 2 * Math.PI * (s1 * s2) / Convert.ToDouble(N));
}
X[k] = X[k] / N;
// Пропустить дополнительный расчет для центрального элемента
, если (k != N/2)
{
X[N - k - 1] = новый комплекс(X[k]. Вещественный, -X[k].Мнимое);
}
}
return X;
}