Inlägg

Pretty Fast FFT https://github.com/marton78/pffft/blob/master/fftpack.h
Denna har cffti, cfftf, cfftb och med dessa så kanske jag kan göra en FFT 1D på en 2D matris.
Men då är frågan hur jag skapar min komplexa array?
Hur ska arrayen c vara? Jag vet att den ska vara komplex, men om den ska vara n lång, så måste den väll vara 2*n lång om den ska vara komplex?

Jag har svårt att tolka hur arrayen c ska skapas. Det är av datatypen float. Så jag misstänker att antingen ska arrayen vara 2*n, eller så så blir det bara halva frekvensen?

Jag misstänker att FFTW3 använder sig av hårdvaran?

Är det någon som vet hur FFTPack är normaliserad? Enligt manualen så är den normaliserad. Men jag har ett behov utav att få den mer lik matlabs.

Jag är nyfiken på vilken mapp du pratar om då i FFTW3.

Men FFTW innehåller så mycket filer. Jag gillar när det är bara några .c filer som man kan bara klistra in. Då slipper man allt görm med att länka hit och dit.

Jag kollade upp FFT för MKL och den kräver C99 för FFT complex. Jag försöker alltid att vara halvt besvärlig igenom att strikt hålla mig till ANSI C. Detta har med att ANSI C är liksom oändlig och kod som är skrivet i ANSI C blir för evig.

Jag har löst detta problem nu. Jag bara använde cfft2f/b och multiplicerade varje tal med n*n. Då fick jag samma resultat som Matlab. Om det är värt det, är tveksamt. Det är ju bara en skalning.

Hur som helst så kanske det är värt att testa MLK för FFT, nu när FFTpack 5.1 är något långsamare än FFTPack 5.1.

Hur ser deras komplexa rutin ut? Kräver den C99?

Jag behöver tydligen inte normalisera när jag kör FFTPack 5.1, men kör jag FFTPack 5.0 (Det första optimerade) https://github.com/marton78/pffft/blob/master/fftpack.c så måste jag normalisera efter invers FFT.

Vi kan ta ett arbetsexempel.

#define row 3
#define column 3

int main() {
	clock_t start, end;
	float cpu_time_used;
	start = clock();

	/* Create X */
	float X[row * column] = { -0.6485,   0.2851,   1.3475,
					       -0.6743, -0.1499, -1.5549,
						1.4951,   0.4504, -0.4982 };

	/* Do FFT2 */
	fft2(X, row, column);

	/* Print */
	print(X, row, column);

	end = clock();
	cpu_time_used = ((float)(end - start)) / CLOCKS_PER_SEC;
	printf("\nTotal speed  was %f\n", cpu_time_used);

	return EXIT_SUCCESS;
}

Utskriften blir

0.0058111       0.0258111       -0.1242458
0.1611111       -0.2950723      0.3444953
0.3681955       -0.2190056      0.0864390

Kör jag

ifft2(X, row, column);

direkt efter

fft2(X, row, column);

så får jag samma matris X tillbaka.

Men om det inte spelar någon roll, så låter jag det vara.
Det viktigaste är att jag kan återskapa datat.

Jag gillar dock inte MLK's FFT implementering. Men den så måste man antingen välja complex.h, att inkludera, eller så måste man ha +2 längre data array för att andra elementet måste ha en noll-imaginär tal, vilket är rätt värdelöst. Då början och slutet på en FFT-array i frekvensplanet är alltid ett reellt tal.

Nu har jag testat FFT2 med FFTPack 5.1. Det fungerar.
Det blir dock inte samma MATLAB-resultat vid tidsplan -> frekvensplan.
Men det blir samma resultat vid tidsplan -> frekvensplan -> tidsplan.

Verkar som FFT inte skalar med kärnor, även med 100 000 samples används bara en kärna både på Ubuntu 22.04 och MacOS (det är då testat mer 3 olika ramverk, FFTW, Arm ComputeLibrary samt Apple Accelerate).

En FYI: FFTW_EXHAUSTIVE/FFTW_MEASURE verkar hänga sig på ARM64 (i alla fall combon Ubuntu server 22.04 + Orange Pi 4, som har 2 st Cortex A72 (samma CPU som RPi4) och 4 st Cortex A53 (samma CPU som RPi3)).

Jag kan göra så här. Jag testar hur snabbt FFTpack 5.1 är, som är skrivet i Fortran 77. Sedan återkopplar jag i denna tråd. Då har vi något att jämföra med

Använde FFTW biblioteket från Homebrew, med det tar det ~0,90ms att köra FFT på en vektor med 100 000 element. Samma sak tar ~0,61ms med PyTorch (på CPU).

Jag vet att FFT3W är något bättre. Men för mig så är 0.001 riktigt snabbt.
Jag har lyckats konvertera om FFTPack 5.1 nu till C-kod. Men det verkar inte finnas någon seriös manual hur man använder FFTPack när det är konverterat från Fortran 77 till ANSI C.

Exempel:

 SUBROUTINE RFFT1F (N, INC, R, LENR, WSAVE, LENSAV, WORK, LENWRK, IER)

 INTEGER    N, INC, LENR, LENSAV, LENWRK, IER
 REAL       R(LENR), WSAVE(LENSAV), WORK(LENWRK)

int rfftf1_(integer *n, integer *in, real *c__, real *ch, real *wa, real *fac);

Den enda manualen jag hittar är: https://github.com/jlokimlin/fftpack5.1/blob/master/doc/FFTPA...

Låter riktigt bra.
Jag gjorde en 100 tusen array med normalfördelande siffror.
0.005000 sekunder tog det att beräkna FFT där.

Med FFTPack. Gör så här. Ladda ned dessa två filer.
https://cs.android.com/android/platform/superproject/+/master...
https://cs.android.com/android/platform/superproject/+/master...

Sedan hittar du på någon array x som har längden n
Kör denna kod

void fft(float x[], size_t n) {
	/* Init */
	float* wsave = (float*)malloc((2 * n + 15) * sizeof(float));
	int* ifac = (int*)malloc((n + 15) * sizeof(float));
	rffti(n, wsave, ifac);

	/* Forward transform */
	rfftf(n, x, wsave, ifac);

	/* Free */
	free(wsave);
	free(ifac);
}

Ja, det gör jag. MLK är känd på internet för att inte vara den snabbaste. FFTW3 är säkert den snabbaste, men FFTPack, trots sin ålder, så är den en av dom snabbaste, och mest lättanvändarvänligaste också.

Enda skillnaden mellan MLK och FFTPack var att FFTPack kunde ge talen 3.99999999999 medan MLK gav 4.0 exakt.

Nu sitter jag med ett problem att jag måste göra en beräkning med 2D med FFTpack. Bara FFTpack 5.1 har stöd för 2D FFT, men denna är skrivet i Fortran 90, och inte ANSI C som jag skriver min kod i.

Då har jag en plan!

1. Om jag har en matris X som är M x N och sedan tar jag transponatet på X för att göra om den till column-major.
2. Jag använder FFT på varje "rad", som är egentligen varje kolumn med tanke på transponantet.

#define row 8
#define column 3

float X[row * column] = { -0.678816,   1.363392, - 0.020819,
							   0.014361,   0.115923,   0.309136,
							  - 0.195106, - 1.376395, - 0.336756,
							   1.678252, - 0.234396, - 0.023950,
							   1.776940,   1.062231, - 0.580292,
							  - 0.723726, - 1.209755,   0.394779,
							  - 0.906230, - 1.120195,   1.207566,
							   2.242360, - 0.311222,   0.433339 };

Efter ha tagit FFT på matrisen X så får jag

3.2080352       -1.7104170      1.3830030
-1.5349649      1.1842327       0.8222664
-0.8341455      -0.7355201      1.9282329
2.1994600       4.9222131       -1.4719210
4.6299772       0.5482139       -0.2945260
-3.3765469      -0.5819106      0.2966796
0.5881023       -1.2479200      -1.1604111
-3.2144592      1.5684826       -0.8436050

I matlab så skulle detta se ut som

>> fft(X)
ans =

   3.2080 +      0i  -1.7104 +      0i   1.3830 +      0i
  -1.5350 - 0.8341i   1.1842 - 0.7355i   0.8223 + 1.9282i
   2.1995 + 4.6300i   4.9222 + 0.5482i  -1.4719 - 0.2945i
  -3.3765 + 0.5881i  -0.5819 - 1.2479i   0.2967 - 1.1604i
  -3.2145 +      0i   1.5685 +      0i  -0.8436 +      0i
  -3.3765 - 0.5881i  -0.5819 + 1.2479i   0.2967 + 1.1604i
   2.1995 - 4.6300i   4.9222 - 0.5482i  -1.4719 + 0.2945i
  -1.5350 + 0.8341i   1.1842 + 0.7355i   0.8223 - 1.9282i

>>

Här får man vara vaksam! Notera att talet 1.5350 + 0.8341i uppträder två gånger i kolumnen och att datat är speglat. FFTPack ger oss endast halva spegeln.

   3.2080 +      0i  -1.7104 +      0i   1.3830 +      0i
  -1.5350 - 0.8341i   1.1842 - 0.7355i   0.8223 + 1.9282i
   2.1995 + 4.6300i   4.9222 + 0.5482i  -1.4719 - 0.2945i
  -3.3765 + 0.5881i  -0.5819 - 1.2479i   0.2967 - 1.1604i
  -3.2145 +      0i   1.5685 +      0i  -0.8436 +      0i

Och då är frågan: Hur ska jag använda FFT här på denna data för att göra 2D FFT?
Jo! Det jag måste göra, är att jag måste ta FFT för varje rad, inte kolumn.
Men det ställer till lite problem för mig.

Om jag ska göra FFT på denna rad

-1.5350 - 0.8341i   1.1842 - 0.7355i   0.8223 + 1.9282i

Så måste jag göra FFT på denna rad

-1.5349649      1.1842327       0.8222664
-0.8341455      -0.7355201      1.9282329

Dessa tal skulle fått byta plats. Andra problemet är att dessa är komplexa tal.
Om jag skulle köra FFT på raden ovan så skulle jag få detta resultat

>> fft([-1.5350 - 0.8341i   1.1842 - 0.7355i   0.8223 + 1.9282i])
ans =

   0.4715 + 0.3586i  -4.8451 - 1.7439i  -0.2314 - 1.1170i

Om jag kör med FFT Pack

float Y[6] = { -1.5350, -0.8341,  1.1842, 0.7355,   0.8223 , 1.9282 };
fft(Y, 6);
print(Y, 1, 6);

>>Utskrift: 2.3011003       -2.7266998      2.0788074       -2.3498001      2.7056365       -1.3581001

Ni ser att det fungerar inte att köra komplexa tal med FFTPack.

Körde MLK jämfört med FFTPack....stor skillnad. Det var nog 50% snabbare, men det handlade om dryga 0.01 sekunder.

FFTPack är ju riktig K&R-kod och uppfanns 1982.

Här är ett exempel där x är en 1D vektor av float och n-2 är storleken av 1D vektorn, Typ om vektorn är 12, så ska n vara 10.

	/* Create handle */
	DFTI_DESCRIPTOR_HANDLE handle = NULL;

	/* Prepare handle */
	DftiCreateDescriptor(&handle, DFTI_SINGLE, DFTI_REAL, 1, n);

	/* Commit */
	DftiCommitDescriptor(handle);

	/* Compute forward */
	DftiComputeForward(handle, x);

	/* Delete model*/
	DftiFreeDescriptor(&handle);

Notera att det är bara matrismultiplikationer som MLK är effektivast på. Övriga så som FFT så är en vanlig enkel FFT algoritm med en for-sats betydligt snabbare.

Vad är skillnaden på ström och effekt? Effektsnåla låter inte riktigt lika bra.