[Alg / PHP] 2D FFT

[Cavorka]

Ik was laatst begonnen met wat basic image analysis (plaatjes vergelijken, plaatjes blurren etc.) met PHP (softwarematige, pixel-per-pixel stuff dus). Uiteraard kom je dan vrij snel 'vast' te zitten, en het leek mij wel aardig om de Fourier Transform van een plaatje te berekenen. Enerzijds uit programmeur-technisch oogpunt, anderzijds uit ervaringen met Fourier / Laplace transformaties in mijn studie.
Ook wanneer je plaatjes echt kwalitatief wilt gaan vergelijken is een fourier getransformeerde van je originele plaatje meestal onontbeerlijk.

Nou, allemaal rozegeur en maneschijn: ik op zoek, en kwam vrij snel op de FFT (Fast Fourier Transform) uit, na eerst de pseudo-hard-core DFT (Discrete Fourier Transform) bekeken te hebben.

Beetje Googlen later:
2D FFT, wiskundig en een beetje C source code
DFT uitleg door Wolfram Research (goed!!)

Vooral die eerste link is zeer nice.

Dit is dus de FT van f (genaamd 'F') en de inverse FT (van 'F' terug naar 'f') (met om overklaarbare redenen een j i.p.v. een i als sqrt(-1)):


Hier is een 1D FFT gemaakt in PHP (wat ik dus ook wil doen):
http://www.mathewbinkley.org/fft.php

Impressive, maar niet leuk genoeg, ik wil een 2D en eventueel 3D versie produceren. Op dit moment zit ik nog in lees / nadenk fase, dus erg veel code kan ik nog niet laten zien.
Waar ik ook nog niet helemaal uit ben is hoe het imaginaire deel nou precies zijn rol gaat spelen als ik het wil programmeren... aangezien de input (een plaatje in mijn geval, met als waarde per punt/pixel 0-256) reeel is.

PS: Ik weet dat er downloadbare FTT C libraries / functies zijn, maarja, where's the fun in that en ik wil het dus per se in PHP gaan doen.

[ Voor 15% gewijzigd door Cavorka op 06-02-2004 19:34 ]

[Limhes]

Cavorka schreef op 06 februari 2004 @ 19:25:
...
Impressive, maar niet leuk genoeg, ik wil een 2D en eventueel 3D versie produceren. Op dit moment zit ik nog in lees / nadenk fase, dus erg veel code kan ik nog niet laten zien.
Waar ik ook nog niet helemaal uit ben is hoe het imaginaire deel nou precies zijn rol gaat spelen als ik het wil programmeren... aangezien de input (een plaatje in mijn geval, met als waarde per punt/pixel 0-256) reeel is.

Ik snap niet helemaal wat je uberhaupt met een 2D of 3D fft wil. Het principe van een Fourier transformatie is dat je een frequentiespectrum berekent van een bepaald bereik binnen je signaal (in dit geval een plaatje). Nu kun je wel voor alle lijnen een aparte fft gaan doen, maar dan haal je het hele idee onderuit, aangezien de pixels onderling behoorlijk zwaar gecorreleerd zijn en dus niet in lijnen of whatsoever kunnen worden ondergebracht wilt het nog enigszins representatief zijn voor de betekenis van het signaal.

Wat de berekening betreft neem ik aan dat je weet dat je een complexe e-macht kan schrijven als een cosinus voor het reele deel en een sinus voor het imaginaire deel. Vervolgens neem je na transformatie de modulus van elke fft-sample. Zo hou je dus uiteindelijk geen complexe waardes meer over (zou je toch gehad moeten hebben op school, als je ook al Fourier en Laplace hebt gehad?).

Disclaimer: heb geen verstand van fft's over plaatjes, alleen over audio, vandaar dat ik er hier en daar wat naast kan zitten.

[ Voor 6% gewijzigd door Limhes op 06-02-2004 19:52 ]

[Cavorka]

Ehm... dude, ik ben 3.5e jaars student natuur/sterrenkunde, ben doodgegooid met de FT's / LT's.

Wat je met een 2D FFT kan? Ja wat niet, dat is meer de vraag. JPG is bijvoorbeeld gebaseerd op FFT (eigenlijk FCT). Alle beeldanalyse gaat zo'n beetje met de FFT van het origineel, noise reductie + compressie ook. Zonder FFT ben je nergens. Maar besides: het maakt niet zoveel uit wat ik ermee wil, ik wil het gewoon maken.

Het nemen van absolute waarde heb ik nog nergens terug gezien, maar dit is inderdaad wel misschien de manier.

[ Voor 14% gewijzigd door Cavorka op 06-02-2004 19:57 ]

[Limhes]

Ah kijk, ik zit blijkbaar _onder_ je niveau, en niet erboven Maar zoals de disclaimer al vermeldde doe ik alleen fft's op audio signalen, misschien dat ik daardoor niet inzag dat je er meer van wist dan ik dacht.

Wat betreft absolute waarde (stelling van Pythagoras):

absValue = sqrt(imagValue .^ 2 * realValue .^ 2);

Maar dit lijkt me dan toch wel weer heel basic...

[Cavorka]

No problem.

Het vage vind ik dus ook dat bijvoorbeeld in de source van deze gast helemaal geen cos noch sin voorkomt.
Hier dan weer wel.

Ik zoek nog even verder, of anders ga ik gewoon proggen van scratch.

Hehe, vooral het ontbreken van enig comment in beide sources is nogal frustrerend, want al dat gebit-shift maakt het totaal onleesbaar (comments als /* Do the bit reversal */ helpen niet echt ).

[ Voor 38% gewijzigd door Cavorka op 06-02-2004 20:14 ]

[Limhes]

Cavorka schreef op 06 februari 2004 @ 20:10:
Hehe, vooral het ontbreken van enig comment in beide sources is nogal frustrerend, want al dat gebit-shift maakt het totaal onleesbaar (comments als /* Do the bit reversal */ helpen niet echt ).

Nee inderdaad, en de duidelijke naamgeving van de variabelen maakt het er ook nog eens een stuk fijner op.
Ik het trouwens wel het vermoeden dat de sinus en cosinus termen mbv. reeksontwikkeling worden berekend, maar zeker weten doe ik t niet, heb niet zo'n analyse-knobbel.

[.oisyn]

Cavorka schreef op 06 februari 2004 @ 19:25:
(met om overklaarbare redenen een j i.p.v. een i als sqrt(-1))

Een j wordt heel vaak gebruik hoor, bij mij op school gaven ze het ook met een j. Is geloof ik vooral in de natuurkunde en electronica zo, omdat i (of I) wordt gebruikt om de current aan te geven

.edit: oh, misschien begreep ik je opmerking verkeerd, maar dan snap ik niet waarom je zegt dat het onverklaarbaar is

[ Voor 13% gewijzigd door .oisyn op 06-02-2004 22:47 ]

[RickN]

Bit-reversed adressing.

Bij een conventionele implementatie van een FFT wil je de elementen het liefst in een specifieke volgorde adresseren. B.v. bij een 16-punts fft:

0,8,4,12,2,10,6,14,1,9,5,13,3,11,7,15

Als je dit binair opschrijft zul je zien dat het de normale 0,1,...,15 volgorde is, maar dan met de bits gereversed. Omdat dit zoveel gebruikt wordt beschikken DSP's vaak over de mogelijkheid om in deze volgorde een array te adresseren.
Er zijn ook zogenaamde selfsorting FFT's waar deze vorm van adresseren niet nodig is, maar daarvan is de interne structuur niet zo mooi. Een plaatje van een FFT, met zogenaamde butterflies en twiddle factors maakt een hoop duidelijk, dus ik zou daar eens naar opzoek gaan...

[ Voor 76% gewijzigd door RickN op 07-02-2004 01:25 ]

[Cavorka]

Limhes schreef op 06 februari 2004 @ 21:30:
[...]
Ik het trouwens wel het vermoeden dat de sinus en cosinus termen mbv. reeksontwikkeling worden berekend, maar zeker weten doe ik t niet, heb niet zo'n analyse-knobbel.

Zou kunnen, maar dat zuo wel heel vaag zijn, aangezien de cos / sin functies in programma's zelf al gerepresenteerd worden door reeksen (Chebychev polynomen (argh, de spelling)). Lijkt me makkelijker gewoon de ingebouwde functie te gebruiken...

.oisyn schreef op 06 februari 2004 @ 22:46:
[...]Een j wordt heel vaak gebruik hoor, bij mij op school gaven ze het ook met een j. Is geloof ik vooral in de natuurkunde en electronica zo, omdat i (of I) wordt gebruikt om de current aan te geven
.edit: oh, misschien begreep ik je opmerking verkeerd, maar dan snap ik niet waarom je zegt dat het onverklaarbaar is

'Heel vaak', . Maybe, maar zoals je zegt niet in de wis- en natuurkunde waar het vandaan komt. En de J wordt juist weer gebruikt voor stroomdichtheid in de natuurkunde...

RickN schreef op 07 februari 2004 @ 00:21:
Bit-reversed adressing.

Bij een conventionele implementatie van een FFT wil je de elementen het liefst in een specifieke volgorde adresseren. B.v. bij een 16-punts fft:

0,8,4,12,2,10,6,14,1,9,5,13,3,11,7,15

(...)

Thanks voor de uitleg. Ik had inderdaad ergens zoiets gelezen (van de vorm: 000, 001, 010, 011, 100, 101, 110, 111 naar 000, 100, 010, 110, 001, 101, 011, 111 wat jij dus zegt).

[ Voor 15% gewijzigd door Cavorka op 07-02-2004 09:59 ]

[Verwijderd]

Limhes schreef op 06 februari 2004 @ 20:04:
Ah kijk, ik zit blijkbaar _onder_ je niveau, en niet erboven Maar zoals de disclaimer al vermeldde doe ik alleen fft's op audio signalen, misschien dat ik daardoor niet inzag dat je er meer van wist dan ik dacht.

Wat betreft absolute waarde (stelling van Pythagoras):

absValue = sqrt(imagValue .^ 2 * realValue .^ 2);

Maar dit lijkt me dan toch wel weer heel basic...

Het moet trouwens sqrt(imagValue^2 + realValue^2) zijn!

Hier heb je toch geen nachten van wakker gelegen toch

[ Voor 13% gewijzigd door Verwijderd op 16-04-2004 21:06 ]