Обратный БПФ в C#

Question

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

Я хочу сделать это с помощью быстрого преобразования Фурье (FFT). Это мой рекурсивный C# FFT:

void ft(float n, ref Complex[] f) 
{ 
    if (n > 1) 
    { 
     Complex[] g = new Complex[(int) n/2]; 
     Complex[] u = new Complex[(int) n/2]; 

     for (int i = 0; i < n/2; i++) 
     { 
      g[i] = f[i * 2]; 
      u[i] = f[i * 2 + 1]; 
     } 

     ft(n/2, ref g); 
     ft(n/2, ref u); 

     for (int i = 0; i < n/2; i++) 
     { 
      float a = i; 
      a = -2.0f * Mathf.PI * a/n; 
      float cos = Mathf.Cos(a); 
      float sin = Mathf.Sin(a); 
      Complex c1 = new Complex(cos, sin); 
      c1 = Complex.Multiply(u[i], c1); 
      f[i] = Complex.Add(g[i], c1); 

      f[i + (int) n/2] = Complex.Subtract(g[i], c1); 
     } 
    } 
} 

Вдохновляющий пример был

Я сравнивались мои результаты с теми из wolframalpha для того же входа 0.6,0.7,0.8,0.9 но результаты не должны быть одинаковыми. Мои результаты в два раза больше, чем Вольфрам, а мнимая часть - это 2 раза от Вольфрама.

Кроме того, вики показывает, что обратное БПФ может быть вычислена с

Но я сравниваю входы и выходы, и они разные.

Есть ли у кого-нибудь идеи, что не так?

Answer 1

Различные реализации часто используют разные определения дискретного преобразования Фурье (ДПФ) с соответствующими результатами. Соответствие между реализациями обычно довольно тривиально (например, коэффициент масштабирования).

Более конкретно, ваша реализация основана на следующем определении ДПФ:

С другой стороны, Вольфрам альфа по умолчанию использует a definition, что после корректировки на 0 на основе индексации выглядит :

Соответственно, можно преобразовать результат вашей реализации, чтобы соответствовать Wolfram Alpha с:

void toWolframAlphaDefinition(ref Complex[] f) 
{ 
    float scaling = (float)(1.0/Math.Sqrt(f.Length)); 
    for (int i = 0; i < f.Length; i++) 
    { 
    f[i] = scaling * Complex.Conjugate(f[i]); 
    } 
} 

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

вы предоставили бы:

void inverseFt(ref Complex[] f) 
{ 
    for (int i = 0; i < f.Length; i++) 
    { 
    f[i] = Complex.Conjugate(f[i]); 
    } 
    ft(f.Length, ref f); 
    float scaling = (float)(1.0/f.Length); 
    for (int i = 0; i < f.Length; i++) 
    { 
    f[i] = scaling * Complex.Conjugate(f[i]); 
    } 
} 

Вызов ft на исходная последовательность 0.6, 0.7, 0.8, 0.9 должна, таким образом, получить преобразованную последовательность 3, -0.2+0.2j, -0.2, -0.2-0.2j.

Дополнительный вызов inverseFt на этой последовательности преобразований должен привести вас к исходной последовательности 0.6, 0.7, 0.8, 0.9 (в пределах некоторой разумной ошибки с плавающей запятой), как показано в this live demo.