Utworzenie szybkiej transformaty Fouriera (FFT) w C++ z wykorzystaniem programu z dyskretną transformatą Fouriera (DFT)

Witam,
Piszę do Was z prośbą o pomoc z utworzeniem szybkiej transformaty Fouriera (FFT) w Microsoft Visual Studio w języku C++. Mianowicie posiadam już poniższy kod z działającą DFT i pozostaje tylko kwestia zmian w kodzie (prawdopodobnie w części Obliczanie DFT), żeby uzyskać te same wyniki za pomocą FFT. Bardzo proszę nawet o jakieś drobne wskazówki.

#include <iostream>
#include<cstdlib>
#include<math.h>

using namespace std;

int main()



{
    double Fs = 8000.0;                // czestotliwosc probkowania 
    double Tpr = 1 / Fs;            // okres pomiaru jednej probki
    double N = 8.0;                    // ilosc pomierzonych probek 
    double rodzielczosc = Fs / N;    // rozdzielczosc wykonywania pomiarow
    double t = 0.0;                    // dyskretna dziedzina czasu, zwiększa się co czas pomiaru jednej probki Tpr
    double pi = 3.14159265359;        // stała Pi
    double pomiar[10];                // tablica do przechowywania wynikow obliczen



    double Re_DFT[10];
    double Im_DFT[10];



    double Re_suma_DFTx = 0;
    double Im_suma_DFTx = 0;



    double Re_suma_DFT[10];
    double Im_suma_DFT[10];



    double modul_DFT = 0;



    //...............................................
    // obliczanie próbek zadanego sygnału wejściowego



    cout << "---------------------------------------------" << endl;
    cout << "Pomierzone probki:" << endl << endl;



    for (int i = 0; i < N; i++)
    {
        pomiar[i] = (1 * sin(2 * pi * 1000 * t)) + (0.5 * sin(2 * pi * 2000 * t + (0.75 * pi)));
        cout << "probka " << i << "  " << pomiar[i] << endl;
        t = t + Tpr;
    }



    cout << endl;
    cout << "---------------------------------------------" << endl;
    cout << "wykonano pomiar " << N << " probek:" << endl << endl;
    cout << "Okres probkowania wynosi: " << Tpr << "[s]" << endl << endl;
    cout << "Rozdzielczosc probkowania wynosi : " << Fs / N << "[Hz]" << endl << endl;
    cout << "---------------------------------------------" << endl << endl;



    //.............................................
    // obliczanie DFT..............................



    for (int m = 0; m < N; m++)
    {
        cout << "+++++++++++++++++++++++++++++++++++++++++++++" << endl;
        cout << "Obliczenia DFT o czestotliwosci  " << m * rodzielczosc << " kHz" << endl;
        cout << "+++++++++++++++++++++++++++++++++++++++++++++" << endl;
        for (int n = 0; n < N; n++)
        {
            Re_DFT[n] = pomiar[n] * (cos(2 * pi * n * m / N));
            Im_DFT[n] = -1 * pomiar[n] * (sin(2 * pi * n * m / N));



            // zaokrąglam obliczone wartości Re_DFT[n] oraz Im_DFT[n] aby uniknąć błędu wynikającego z obliczeń cyfrowych
            if (Re_DFT[n] > -0.001 && Re_DFT[n] < 0.001)
            {
                Re_DFT[n] = 0.0;
            }
            if (Im_DFT[n] > -0.001 && Im_DFT[n] < 0.001)
            {
                Im_DFT[n] = 0.0;
            }



            cout << "Re_DFT" << "[" << m << "]" << "[" << n << "] = " << pomiar[n] << "*cos(2*" << pi << "*" << m << "*" << n << "/" << N << ") = " << Re_DFT[n] << endl;
            cout << "Im_DFT" << "[" << m << "]" << "[" << n << "] = " << pomiar[n] << "*sin(2*" << pi << "*" << m << "*" << n << "/" << N << ") = " << Im_DFT[n] << endl;

            Re_suma_DFTx = Re_suma_DFTx + Re_DFT[n];
            Im_suma_DFTx = Im_suma_DFTx + Im_DFT[n];



            if (n == N - 1)
            {
                Re_suma_DFT[m] = Re_suma_DFTx;
                Im_suma_DFT[m] = Im_suma_DFTx;



                modul_DFT = sqrt((Re_suma_DFT[m] * Re_suma_DFT[m]) + (Im_suma_DFT[m] * Im_suma_DFT[m]));



                if (Re_suma_DFT[m] > -0.001 && Re_suma_DFT[m] < 0.001)
                {
                    Re_suma_DFT[m] = 0.0;
                }
                if (Im_suma_DFT[m] > -0.001 && Im_suma_DFT[m] < 0.001)
                {
                    Im_suma_DFT[m] = 0.0;
                }
                if (modul_DFT > -0.001 && modul_DFT < 0.001)
                {
                    modul_DFT = 0.0;
                }



                cout << "---------------------------------------------" << endl;
                cout << " modul skladowej DFT o czestotliwosci " << m * rodzielczosc << " kHz wynosi: " << modul_DFT << endl;
                cout << " czesc rzeczywista DFT o czestotliwosci " << m * rodzielczosc << " kHz wynosi: " << Re_suma_DFT[m] << endl;
                cout << " czesc urojona DFT o czestotliwosci " << m * rodzielczosc << " kHz wynosi: " << Im_suma_DFT[m] << "j" << endl;
                cout << "---------------------------------------------" << endl << endl;



                Re_suma_DFTx = 0;
                Im_suma_DFTx = 0;
            }
        }
    }



    system("pause");
}

FFT i DFT to tak różne algorytmy, że ciężko przerobić jeden w drugi. Dawniej też stanąłem przed tym samym problemem i stwierdziłem, że lepiej znaleźć algorytm i go dostosować, niż przerabiać posiadany algorytm DFT.

Mam dwa projekty, w których robię FFT o wielkości będącej potęgą 2.
https://github.com/andrzejlisek/MorseCode - fftx.h i fftx.cpp - jest to klasa, która służy tylko do FFT.
https://github.com/andrzejlisek/AudioSpectrum - worker.js - wszystkie zmienne i funkcje dotyczące FFT zaczynają nazwę od FFT_.

Pierwszy z nich jest w C++ i tego projektu już od dawna nie rozwijam, a drugi jest w JavaScript, i czasem ten projekt aktualizuję, kod jest wzięty niemalże żywcem z projektu MorseCode.

W obu projektach (właściwie oba robią to samo, jeśli chodzi o FFT), obliczenia są rozdzielone na część, którą wystarczy przeliczyć jeden raz i część, którą przelicza się przy każdym generowaniu FFT.

Proponuje popatrzeć do projektu https://github.com/andrzejlisek/MorseCode i tam masz praktycznie to, czego potrzebujesz, możesz sobie wziąć żywcem.

@andrzejlisek: Dziękuję Bardzo. Na pewno skorzystam z porad :)

Największym problemem jest to, że piszesz wszystko w main.
Dla programików typu "hello world" to nie jest problem, ale czegoś takiego jak FFT jest to już nie do przyjęcia.

Możesz jeszcze temu się przyglądnąć
https://ng.techweb.com/authds/gateway?gateway=true&return=https%3A%2F%2Fwww.drdobbs.com%2Fcpp%2Fa-simple-and-efficient-fft-implementatio%2F199500857