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

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

12 голосов
/ 19 апреля 2011

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

На основе следующей формулы:

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; Обратное ДПФ не будет иметь -1 в аргументе, в то время как обычное ДПФ действительно имеет -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 секунд, просто скопировав их код.Однако я до сих пор не знаю, что я сделал неправильно в своей реализации.

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

Ответы [ 2 ]

5 голосов
/ 19 апреля 2011

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

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 );
    }
}

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

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

1 голос
/ 21 апреля 2011

Ваш код на самом деле почти правильный (вам не хватает 1 / N в обратном преобразовании). Дело в том, что формула, которую вы использовали, обычно используется для вычислений, потому что она легче, но в чисто теоретических средах (и в Wolfram) вы должны использовать нормализацию на 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 в обратном преобразовании).

Приветствия

...