c# - rapida - Qué está mal con esta implementación de la transformada de Fourier

transformada rapida de fourier en c

Question

Estoy tratando de implementar una transformada discreta de Fourier, pero no está funcionando. Probablemente he escrito un error en alguna parte, pero aún no lo he encontrado.

De acuerdo con la siguiente fórmula:

Esta función realiza el primer ciclo, girando sobre X0 - Xn-1 ...

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

Y el cálculo real, es probablemente donde está el error.

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

A continuación la explicación del resto del código:

var sign = reverse ? 1.0: -1.0; El DFT invertido no tendrá -1 en el argumento, mientras que un DFT regular tiene un -1 en el argumento.

var argument = sign*2.0*Math.PI*k/data.Length; es el argumento del algoritmo. Esta parte:

entonces la última parte

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

Creo que copié cuidadosamente el algoritmo, así que no veo dónde cometí el error ...

Información Adicional

Como ha demostrado Adam Gritt en su respuesta, AForge.net ofrece una buena implementación de este algoritmo. Probablemente pueda resolver este problema en 30 segundos simplemente copiando su código. Sin embargo, todavía no sé lo que hice mal en mi implementación.

Tengo mucha curiosidad sobre dónde está mi error y qué he interpretado mal.

Answer 1

Mis días de hacer matemáticas complejas están un poco por detrás de mí en este momento, por lo que puede que también me esté perdiendo algo. Sin embargo, me parece que estás haciendo la siguiente línea:

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

cuando debería ser más como:

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

A menos que tenga esto envuelto en el método FromPolarCoordinates()

ACTUALIZACIÓN: Encontré el siguiente fragmento de código en la biblioteca de AForge.NET Framework y muestra operaciones Cos / Sin adicionales que no se están manejando en su código. Este código se puede encontrar en contexto completo en el método 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 );
    }
}

Está utilizando una clase compleja personalizada (como era antes de 4.0). La mayoría de las matemáticas son similares a las que ha implementado, pero la iteración interna realiza operaciones matemáticas adicionales en las porciones Real e Imaginario.

NUEVA ACTUALIZACIÓN: Después de algunas implementaciones y pruebas, descubrí que el código anterior y el código proporcionado en la pregunta producen los mismos resultados. También encontré, basado en los comentarios, la diferencia entre lo que se genera a partir de este código y lo que produce WolframAlpha. La diferencia en los resultados es que parecería que Wolfram está aplicando una normalización de 1 / sqrt (N) a los resultados. En el enlace de Wolfram proporcionado si cada valor se multiplica por Sqrt (2), los valores son los mismos que los generados por el código anterior (aparte de los errores de redondeo). Probé esto pasando 3, 4 y 5 valores a Wolfram y encontré que mis resultados eran diferentes por Sqrt (3), Sqrt (4) y Sqrt (5) respetuosamente. Basado en la información de la Transformada de Fourier Discreta proporcionada por wikipedia, menciona una normalización para hacer que las transformaciones para DFT e IDFT sean unitarias. Esta podría ser la vía hacia la que debe mirar hacia abajo para modificar su código o comprender qué puede estar haciendo Wolfram.

Answer 2

Su código es realmente casi correcto (le falta un 1 / N en la transformación inversa). La cuestión es que la fórmula que usaste normalmente se usa para cálculos porque es más ligera, pero en entornos puramente teóricos (y en Wolfram), usarías una normalización de 1 / sqrt (N) para hacer que las transformadas sean unitarias.

es decir, tus fórmulas serían:

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

Es solo una cuestión de convención en la normalización, solo cambian las amplitudes por lo que sus resultados no fueron tan malos (si no hubiera olvidado el 1 / N en la transformación inversa).

Aclamaciones