Biblioteka DFT - dyskretna transformata furiera (ang. FFT)

Witam,
muszę napisać program, wykorzystujący bibliotekę dyskretnej transformaty furiera (ang.FFT)
Ma to wyglądać następująco(kod napisany w C#):

using System;
using System.Collections.Generic;
using System.Linq;
using System.Numerics;
using System.Text;
using System.Threading.Tasks;
using System.Windows;

namespace Adaline
{
    class DFT
    {
        public static double[] Amplituda { get; set; }
        public static Complex[] ciagPoDFT;


        public static double[] przeksztalc(double[] dane)
        {
            int N = dane.Count();

            Complex[] x = new Complex[N];
            for (int i = 0; i < N; i++)
            {
                x[i] = new Complex(dane[i], 0);
            }

            MathNet.Numerics.IntegralTransforms.Fourier.Forward(x);
         
            Amplituda = new double[N];

            for (int i = 0; i < N; i++)
            {
              Amplituda[i] = Complex.Abs(x[i]);
            }

            return Amplituda;

        }
    }
}

Jest tutaj użyta dana biblioteka dla C#.
Wykorzystywana w tym miejscu:

MathNet.Numerics.IntegralTransforms.Fourier.Forward(x);

Chcę przerobić powyższy kod z C# na Delphi.
Na stronie znajduje się biblioteka, lecz nie jestem w stanie użyć jej tak jak w powyższym przykładzie.

http://www.alglib.net/fasttransforms/fft.php

Ktoś byłby mi w stanie pomóc przy zamianie kodu na Delphi?
Bardzo proszę o pomoc.

Jeśli posiadasz liczbę próbek równą potędze 2 to skorzytaj z tego modułu:
https://github.com/teragonaudio/Convolver/blob/master/Third-Party/FFTReal-2.11/delphi/fftreal.pas

Przykład użycia:

 
procedure FFT;
var
  n, i:integer;
  TabOut, TabIn: pflt_array;
  FFTReal: TFFTReal;
begin
    GetMem(TabOut, n * sizeof_flt);
    GetMem(TabIn, n * sizeof_flt);
 
    //uzupełnianie tablicy wejściowej
    for i := 0 to n - 1 do
    TabIn[i] := random(500);

    FFTReal := TFFTReal.Create(n);
    FFTReal.do_fft(@TabOut[0], @TabIn[0]);
    FreeAndNil(FFTReal);
	
    for i := 0 to (n div 2) - 1 do
    begin
    // obliczanie amplitudy
    end;
	
    FreeMem(TabOut);
    FreeMem(TabIn);
end;

FFTReal.do_fft daje w wyniku połowę próbek, bo wykres FFT jest symetryczny, dlatego w TabOut do indeksu równego n div 2 - 1 znajduje się część rzeczywista, a dalej są odpowiednie części urojone.

Jak nie znasz liczby próbek i nie chcesz uzupełniać ich zerami do potęgi dwójki to możesz wykorzystać tą DLLke:
http://www.fftw.org/install/windows.html

Przykład użycia i interfejs znajdziesz tutaj:
ftp://ftp.fftw.org/pub/fftw/fftw3-delphi.zip

A więc przekazuje do tego algorytmu obrazek w postaci liczb zespolonych i na wyjściu dostaje tablice tej samej długości po przekształceniu.

**Próbek **tzn jako parametr ta procedura miała by otrzymywać tablice 25-elementową:

 dane[24]

.
W takim razie można użyć pierwszego przykładu ?

Ten pierwszy moduł nadaje się tylko do oblicza FFT z liczb rzeczywistych. Dodatkowo wymaga, żeby liczba próbek/elementów wejściowych była jedną z potęg liczby 2. Jak chcesz liczyć FFT z dowolnej ilości próbek/elementów liczb zespolonych to skorzystaj z tej drugiej biblioteki. Tutaj znajdziesz więcej informacji: http://www.fftw.org/fftw3_doc/Complex-DFTs.html#Complex-DFTs

//Edit
Dopiero teraz zauważyłem, że część urojona na wejściu zawsze jest równa 0, czyli na wejściu masz dane rzeczywiste. Jak dodasz 7 zer na końcu danych wejściowych to będziesz mógł skorzystać z pierwszego modułu/biblioteki.

A jakiej użyc biblioteki odnośnie tego pierwszego przykładu, ponieważ przy deklaracji tablic typu pflt_array oraz TFFTreal zaznacza mi błąd:
undeclared identifier

Do sekcji uses musisz dodać nazwę modułu, czyli w tym przypadku "FFTReal". Z tym że w module musisz jeszcze zmienić deklarację typu

flt_t = single;

na flt_t = Double;

Robiłeś coś kiedykolwiek w delphi albo pascalu? Bo dołączanie modułów i błąd "undeclared identifier" to podstawy podstaw.

Przepraszam za moje niedopatrzenie.
Nie zwróciłem uwagi na link, który podałeś.
Przetestuję i mam nadzieję że zadziała.
Dziękuję za pomoc!

@Tajiri
A jeśli mógłbyś napisać,czy ta linijka

 FFTReal.do_fft(@TabOut[0], @TabIn[0])

Odpowiada za to samo co to(C#):

MathNet.Numerics.IntegralTransforms.Fourier.Forward(x); 

Bo w przykładzie z C# jest przekazywana cała tablica z wszystkimi danymi.

Czy należałoby zrobić coś takiego

for i:=0 to length(TabIn) do
begin
  FFTReal.do_fft(@TabOut[i], @TabIn[i])
end;
 

I jeszcze chciałbym zapytać,jak w tym przypadku miałaby wyglądać amplituda?
Czy tak?

 
TabOut[i] = Abs(TabIn[i]);

Patrząc na przykład z C#.

Metoda "do_fft" przyjmuje dwa parametry: wskaźnik na tablicę wyjściową i wskaźnik na tablicę wejściową, dlatego wpisane jest "@TabOut[0], @TabIn[0]". Znaczek "@" oznacza wskaźnik, więc cały zapis "@TabOut[0]" to wskaźnik na pierwszy element tablicy wyjściowej. Analogicznie zapis "@TabIn[0]" daje wskaźnik na pierwszy element tablicy wejściowej. Rozmiar obu tablic ustawiany jest w konstruktorze - linijka "TFFTReal.Create(n)". Za rozmiar i typ elementu tablicy odpowiada linijka "flt_t = Double;" w module "FFTReal.pas"

Wydaje mi się, że ta metoda "MathNet.Numerics.IntegralTransforms.Fourier.Forward" liczy normalną transformatę fouriera z tablicy danych o dowolnej długości i w wyniku zwraca pełną liczbę próbek. Pełną czyli dane są symetryczne, co wynika z samej transformaty.

Ta linijka "FFTReal.do_fft(@TabOut[0], @TabIn[0]);" nie jest do końca równoważna. Pierwsza różnica to liczba danych wejściowych. W metodzie "do_fft" musi być równa potędze dwójki, dlatego na końcu danych musisz dodać kilka zer. W danych wynikowych spowodują one zwiększenie rozdzielczości FFT. Druga różnica to sama konstrukcja danych wyjściowych. Tak jak pisałem wcześniej FFTReal.do_fft daje w wyniku połowę próbek, bo wykres FFT jest symetryczny, dlatego w TabOut do indeksu równego n div 2 - 1 znajduje się część rzeczywista, a dalej są odpowiednie części urojone.

Delphi automatycznie nie wspiera liczb zespolonych, więc amplitudy musisz policzyć ręcznie:

// Abs(a+ib) = sqrt(a^{2 + b}2)

function GetFFTAmplitude(realPart: Double; imaginaryPart: Double): Double;
begin
	  Result := Sqrt(Sqr(realPart) + Sqr(imaginaryPart));
end;

SetLength(AbsTab, n);
for i := 0 to (n div 2) - 1 do
begin
	AbsTab[i] := GetFFTAmplitude(TabOut[i], TabOut[(n div 2) + i]);
	AbsTab[n - i - 1] := AbsTab[i]; // załatwia sprawę z symetrią wykresu
end;

Po wyjściu z pętli powinieneś dostać symetryczny wykres FFT, ale ze zwiększoną rozdzielczością, więc nie będzie identyczny z tym w C#. Z punktu widzenia samej analizy FFT takie dane będą poprawne, ale jak masz przepisać ten kod jeden do jednego to musisz skorzystać z biblioteki FFTW.

Nie musi być jeden do jednego.
Tablica dane posiada tylko wartości -1; 1; -1; 1; -1...
A więc mam taki kod:

procedure FFT.przeksztalc(var dane : array of double);
var
  n, i:integer;
  TabOut, TabIn: pflt_array;
  FFTReal: TFFTReal;
begin
    n:=length(dane);
    GetMem(TabOut, n * sizeof_flt);
    GetMem(TabIn, n * sizeof_flt);
 
    //uzupełnianie tablicy wejściowej
    for i := 0 to n - 1 do
    begin
       TabIn[i] := dane[i]+0.0000000;
    end;
 
    FFTReal := TFFTReal.Create(n);
    for i:=0 to N-1 do 
    begin
       FFTReal.do_fft(@TabOut[i], @TanIn[i]);
    end;

    setLength(Amplituda,N);
    for i:=0 to (N div 2) -1 do
    begin
       Amplituda[i]:=GetFFTAmplitude(TabOut[i],TabOut[ (n div 2) + i]);
       Amplituda[n -i -1] := Amplituda[i];
    end;
end;
 

Otrzymuję wartość Amplitudy bardzo dużą np:13156,5789118707 ; 79495,2368...
Tak samo po wywołaniu

FFTReal.do_fft(@TabOut[i], @TanIn[i]);
 

Otrzymuję w TabOut:
-3; -36,0406..; -11,2938..; -32,6561..; -127,1033..;

W C# amplituda jest bliska 0.
0,360.. ;

Czy ta funkcja jest w porządku napisana?

Z tymi zerami to chodziło mi o uzupełnianie tabeli o dodatkowe zera. Pisałeś, że tabela "dane" ma 25 elementów, więc musisz dodać 7 zer. (-1; 1; -1; 1; -1; ... ; 0; 0; 0; 0; 0; 0; 0)

function FFT.przeksztalc(var dane : array of double):Double;
var
  n, i:integer;
  TabOut, TabIn: pflt_array;
  FFTReal: TFFTReal;
begin
	n:=length(dane) + 7;
    GetMem(TabOut, n * sizeof_flt);
    GetMem(TabIn, n * sizeof_flt);
 
	// zerowanie tablic
	for i := 0 to n - 1 do
	begin
		TabOut[i] := 0;
		TabIn[i] := 0;
	end;
	
    //uzupełnianie tablicy wejściowej
    for i := 0 to length(dane) - 1 do
    begin
       TabIn[i] := dane[i];
    end;
	
	FFTReal := TFFTReal.Create(n);
    FFTReal.do_fft(@TabOut[0], @TabIn[0]);
 
    setLength(Amplituda,N);
    for i:=0 to (N div 2) -1 do
    begin
       Amplituda[i]:=GetFFTAmplitude(TabOut[i],TabOut[ (n div 2) + i]);
       Amplituda[n -i -1] := Amplituda[i];
    end;

    // zwalnianie pamięci
    FreeMem(TabOut);
    FreeMem(TabIn);
    FreeAndNil(FFTReal);
end;

Co do samych wartości to zrób prosty test: w C# dla znanych danych wejściowych zobacz jakie wartości dostajesz po wyjściu z metody liczącej FFT. Dla tych samych danych wejściowych zobacz co dostajesz w Delphi. Wartości powinny być podobne tzn. różne wartości ale mniej więcej tego samego rzędu. Być może metoda z C# robi wewnętrznie jakąś normalizacje tego nie wiem, bo z niej nie korzystałem.