Что не так с этой реализацией преобразования Фурье

Я пытаюсь реализовать дискретное преобразование Фурье, но он не работает. Я, наверное, где-то написал ошибку, но еще не нашел ее.

Исходя из следующей формулы:

terere

Эта функция выполняет первый цикл, перейдя через X0 - Xn-1... enter image description here

    public Complex[] Transform(Complex[] data, bool reverse)
    {
        var transformed = new Complex[data.Length];
        for(var i = 0; i < data.Length; i++)
        {
            //I create a method to calculate a single value
            transformed[i] = TransformSingle(i, data, reverse);
        }
        return transformed;
    }

И фактическое вычисление, вероятно, там, где ошибка.

    private Complex TransformSingle(int k, Complex[] data, bool reverse)
    {
        var sign = reverse ? 1.0: -1.0;
        var transformed = Complex.Zero;
        var argument = sign*2.0*Math.PI*k/data.Length;
        for(var i = 0; i < data.Length; i++)
        {
            transformed += data[i]*Complex.FromPolarCoordinates(1, argument*i);
        }
        return transformed;
    }

Далее объяснение остальной части кода:

var sign = reverse ? 1.0: -1.0; Обратный DFT не будет иметь -1 в аргументе, тогда как обычный DFT имеет аргумент -1.

enter image description here

var argument = sign*2.0*Math.PI*k/data.Length; является аргументом алгоритма. Эта часть:

enter image description here

то последняя часть

transformed += data[i]*Complex.FromPolarCoordinates(1, argument*i);

Я думаю, что я тщательно скопировал алгоритм, поэтому не вижу, где я ошибся...

Дополнительная информация

Как показал Адам Гритт в своем ответе, есть хорошая реализация этого алгоритма AForge.net. Возможно, я смогу решить эту проблему за 30 секунд, просто скопировав их код. Тем не менее, я до сих пор не знаю, что я сделал неправильно в моей реализации.

Мне действительно любопытно, где мой недостаток, и что я неправильно интерпретировал.

Ответы

Ответ 1

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

transformed += data[i]*Complex.FromPolarCoordinates(1, argument*i);

когда это должно быть больше похоже:

transformed += data[i]*Math.Pow(Math.E, Complex.FromPolarCoordinates(1, argument*i));

Если вы не включили этот метод в метод FromPolarCoordinates()

UPDATE: Я нашел следующий бит кода в библиотеке AForge.NET Framework, и он показывает дополнительные операции Cos/Sin, которые выполняются не в вашем код. Этот код можно найти в полном контексте в методе Sources\Math\FourierTransform.cs: DFT.

for ( int i = 0; i < n; i++ )
{
    dst[i] = Complex.Zero;

    arg = - (int) direction * 2.0 * System.Math.PI * (double) i / (double) n;

    // sum source elements
    for ( int j = 0; j < n; j++ )
    {
        cos = System.Math.Cos( j * arg );
        sin = System.Math.Sin( j * arg );

        dst[i].Re += ( data[j].Re * cos - data[j].Im * sin );
        dst[i].Im += ( data[j].Re * sin + data[j].Im * cos );
    }
}

Используется специальный класс Complex (как и до 4.0). Большая часть математики похожа на то, что вы реализовали, но внутренняя итерация выполняет дополнительные математические операции над реальными и мнимыми частями.

ДАЛЬНЕЙШЕЕ ОБНОВЛЕНИЕ: После некоторой реализации и тестирования я обнаружил, что приведенный выше код и код, предоставленные в вопросе, дают одинаковые результаты. Я также нашел, основываясь на комментариях, какая разница между тем, что генерируется из этого кода, и тем, что производится WolframAlpha. Разница в результатах заключается в том, что, по-видимому, Вольфрам применяет нормировку 1/sqrt (N) к результатам. В Wolfram Link предоставляется, если каждое значение умножается на Sqrt (2), то значения те же, что и те, которые генерируются вышеуказанным кодом (в сторону округления ошибок). Я проверил это, передав 3, 4 и 5 значениям в Wolfram и обнаружил, что мои результаты были разными по Sqrt (3), Sqrt (4) и Sqrt (5). На основе Discrete Fourier Transform информация, предоставленная wikipedia, упоминает нормализацию, чтобы сделать преобразования для DFT и IDFT унитарными. Это может быть проспект, что вам нужно посмотреть, чтобы изменить код или понять, что может делать Вольфрам.

Ответ 2

Ваш код фактически почти прав (вам не хватает 1/N в обратном преобразовании). Дело в том, что используемая вами формула обычно используется для вычислений, потому что она легче, но в чисто теоретических средах (и во Вольфраме) вы должны использовать нормализацию на 1/sqrt (N), чтобы сделать преобразования унитарными.

то есть. ваши формулы:

Xk = 1/sqrt(N) * sum(x[n] * exp(-i*2*pi/N * k*n))

x[n] = 1/sqrt(N) * sum(Xk * exp(i*2*pi/N * k*n))

Это просто вопрос условности в нормализации, меняются только амплитуды, поэтому ваши результаты не так уж плохи (если вы не забыли 1/N в обратном преобразовании).

Приветствия