[c++] Fast Fourier Transform

[Verwijderd]

Ik probeer een programma te maken dat een opgenomen geluidsgolf weergeeft op het scherm, zowel als golf en als fft (Fast Fourier Transform).
Een fft laat per frequentie de intensiteit van een signaal zien.

Nu ben ik al zover dat ik een geluid kan opnemen en weergeven als golf, maar nu moet ik een fft maken.

Via google heb ik een hele goede site gevonden:
http://www.relisoft.com/Science/Physics/fft.html
Maar hierbij zit zo veel wiskunde dat ik het bos niet meer door de bomen kan zien.
Zo staat erbij bij de c++ code nergens waar ik het invoersignaal moet meegeven / specificeren.
Ook heb ik al veel uur verkloot aan de fft source-code die je kan downloaden, want dit lukt gewoon niet zoals ik het wil.

Het enige wat ik moet hebben is een soort van functie waarbij ik een array mee geef met samples, (en een evt. lengte van die sample-array), zodat de functie een fft-array kan produceren.
Dus zo'n soort header zou die functie kunnen hebben:
void fftTransform(short *input, int n_samples, short *output);
waarbij input mijn sample-array is (dus van ongeveer -32260 tot 32260: dus 16-bits unsigned), n_samples is het aantal samples (bij fft moet het een macht van 2 zijn dacht ik), en output gaat de fft-waarden bevatten.

Begrijpt iemand fft dusdanig zodat diegene mij kan helpen.

[Soultaker]

Fouriertransformaties zijn best ingewikkeld en eigenlijk bijna onmogelijk te begrijpen zonder redelijk diepgaande wiskundige achtergrondkennis (met Wiskunde op middelbare school nivo ben je er nog lang niet). Ik kan je dan ook niet 'even' uitleggen hoe een 'gewone' Fouriertransformatie werkt, laat staan hoe een FFT werkt.

Als je doel simpelweg is om een FFT algoritme te gebruiken, dan is het een kwestie van zoeken naar een geschikte library. In dat geval zou ik graag meer specifiek horen waar je niet uitkomt, want we gaan natuurlijkg een functie zoals jij die beschrijft voor je maken. FFT is gelukkig een vrij standaard algoritme dus met Google is genoeg code te vinden, lijkt me.

[Verwijderd]

Dat is het hele punt: ik hoef/wil het niet te begrijpen. Maar ik weet wel wat een fft voor resultaat geeft, en daar gaat het mij om.

Dus, als ik een array heb van samples dus:

short samples[1024]


Dit stop ik in een functie en ik wil een array krijgen die de fft-waarden bevat.
Het hoeft niet eens perse broncode te zijn, als ik een library oid kan gebruiken vind ik het ook goed

[alienfruit]

Volgens mij worden die Fast Fourier Transforms ook gebruikt bij Kalman Filtering. Daar snap ik zelf na drie maanden nog steeds de ballen van

[Soultaker]

Afterdawn: Nou, doe eens Google'en dan; dan vind je al snel iets als FFTW: "FFTW is a C subroutine library for computing the discrete Fourier transform (DFT) in one or more dimensions, of arbitrary input size, and of both real and complex data (as well as of even/odd data, i.e. the discrete cosine/sine transforms or DCT/DST). We believe that FFTW, which is free software, should become the FFT library of choice for most applications."

alienfruit: Ik geloof dat zo'n beetje alle signaalfilters wel op de een of andere manier wat met Fouriertransformaties te maken hebben.

[Verwijderd]

Bedankt! Ik had al wel gezocht alleen deze nog niet gevonden! Staan ook instructies bij voor hoe en wat met visual c++. Ik zal het snel uitproberen, nu ga ik slapen.

[Zoijar]

Numerical Recipes in C, heel hoofdstuk erover (12): http://www.library.cornell.edu/nr/bookcpdf.html

Een fourier transformatie is eigenlijk gewoon een overgang op een (al dan niet oneindige) lineaire basis, van sin en cos, of complexe e-machten wat hetzelfde is. (
(Je kan dan bv de extreem hoog frequente sinussen eruit filteren omdat die toch weinig aan je signaal toevoegen, en dan terug transformeren (het was een lineaire xform, dus heeft een inverse). Verliest iets aan signaal, maar compressed wel beter.)

Dan krijg je de DFT, die kan werken met eindige discrete data, zoals computer samples. Maar die heeft complexiteit O(N²).

Dan heb je de FFT die hetzelfde doet als DFT, maar van order O(n log n) is. Sneller dus voor grote data.

Als je er over wilt lezen moet je het in deze volgorde doen. Op zich hoef je het niet diep wiskundig door te spitten, maar het is wel handig als je een beetje een idee hebt van wat je eigenlijk doet.

[ Voor 82% gewijzigd door Zoijar op 09-04-2004 11:17 ]

[Verwijderd]

zo en zo mis ik ook nog een imagionair gedeelte in je FFT, en je zal het nooit per frequentie krijgen ten minste, meestal in frequentie banden, veelgaand een stuk of 16 etc. Maar als je er echt effectief mee om wil gaan en met signalen in het algemeen zou ik toch maar wat er over lezen. En misschien is matlab wel eens een aardig programma om te bekijken.

[Verwijderd]

Zoijar schreef op 09 april 2004 @ 10:48:
Numerical Recipes in C, heel hoofdstuk erover (12): http://www.library.cornell.edu/nr/bookcpdf.html

Een fourier transformatie is eigenlijk gewoon een overgang op een (al dan niet oneindige) lineaire basis, van sin en cos, of complexe e-machten wat hetzelfde is. (
(Je kan dan bv de extreem hoog frequente sinussen eruit filteren omdat die toch weinig aan je signaal toevoegen, en dan terug transformeren (het was een lineaire xform, dus heeft een inverse). Verliest iets aan signaal, maar compressed wel beter.)

Dan krijg je de DFT, die kan werken met eindige discrete data, zoals computer samples. Maar die heeft complexiteit O(N²).

Dan heb je de FFT die hetzelfde doet als DFT, maar van order O(n log n) is. Sneller dus voor grote data.

Als je er over wilt lezen moet je het in deze volgorde doen. Op zich hoef je het niet diep wiskundig door te spitten, maar het is wel handig als je een beetje een idee hebt van wat je eigenlijk doet.

Ok dit had ik nodig, bedankt en nu kom ik er wel uit!

[Verwijderd]

Afterdawn....je bent duidelijk een grote optimist.

Als je gebruik gaat maken van de routines uit Numerical recipes test dan zeker met bekende functies waarvan je weet wat de uitkomst moet zijn. Ik heb ze ooit geimplementeerd in C++ en het heeft me toen erg veel hoofdbrekens gekost om uit te vogelen hoe de data in die routine moet worden gestopt en hoe het er vervolgens weer uit komt.
(en bij een aantal van die routines kreeg ik althans de uitkomst niet kloppend)

Ik vond dat boek wat dat betreft buitengewoon onduidelijk.
En ik weet genoeg van FFTs om te snappen wat ik aan het doen ben.

[ Voor 8% gewijzigd door Verwijderd op 14-04-2004 21:46 ]

[Janoz]

Verwijderd schreef op 09 april 2004 @ 10:57:
zo en zo mis ik ook nog een imagionair gedeelte in je FFT, en je zal het nooit per frequentie krijgen ten minste, meestal in frequentie banden, veelgaand een stuk of 16 etc. Maar als je er echt effectief mee om wil gaan en met signalen in het algemeen zou ik toch maar wat er over lezen. En misschien is matlab wel eens een aardig programma om te bekijken.

Het imaginaire deel wordt alleen gebruikt voor de fase. Voor een leuke equalizer scherm heb je alleen de amplitude van het signaal nodig. En dat wordt in het reele deel opgeslagen.

Over die banden heb je inderdaad wel gelijk .

[Verwijderd]

Janoz, als je een perfecte sinus aanbied. Dus een sinus (die dus start op 0 in je data) is de FFT volledig imaginair (het reele deel is 0). Cosinus is precies omgekeerd.

Je wil inderdaad de amplitude laten zien in je equalizer scherm, maar dan moet je hier wel even over nadenken. (reele deel is namelijk niet hetzelfde als de amplitude) Helemaal als je dit wil doen met numerical recipies.

ps) en dan heb ik het nog niet eens over windowing effecten bij korte samples waar je toch zeker wel voor moet oppassen.

[ Voor 13% gewijzigd door Verwijderd op 14-04-2004 22:52 ]

[MSalters]

Verwijderd schreef op 14 april 2004 @ 21:45:
Afterdawn....je bent duidelijk een grote optimist.

Als je gebruik gaat maken van de routines uit Numerical recipes test dan zeker met bekende functies waarvan je weet wat de uitkomst moet zijn. Ik heb ze ooit geimplementeerd in C++ en het heeft me toen erg veel hoofdbrekens gekost om uit te vogelen hoe de data in die routine moet worden gestopt en hoe het er vervolgens weer uit komt.
(en bij een aantal van die routines kreeg ik althans de uitkomst niet kloppend)

Ik vond dat boek wat dat betreft buitengewoon onduidelijk.
En ik weet genoeg van FFTs om te snappen wat ik aan het doen ben.

Heb je een recente druk? De oude NR was berucht om het feit dat de C code geporte FORTRAN was, en technisch illegaal. De geheugenallocaties klopten bijvoorbeeld niet, wat kwam omdat er een fout was gemaakt bij de overgang van FORTRAN-stijl 1-based indices naar de C-stijl 0-based. De code is op een gegeven ogenblik gefixt.

Oh, en gebruik NR niet in commerciele programma's. Het ding heeft een licentiepolitiek die zeg maar het slechtste combineert van de GPL en de MS EULA.

[Verwijderd]

@moesasji1975nl:

Ik heb de eerste functie geprobeerd (de four1(...) functie uit hfd 12.2).

Deze functie gaat van 1 tot 2n (opletten: de 0 wordt genegeerd omdat de functie oorspronkelijk van een Fortran-achtige taal komt bij welke arrays beginnen op index 1!).
Hierbij is het eerste reele element op index 1, het eerste imaginaire element op index 2, het tweede reele element staat op index 3, imaginaire op index 4 etc.

Dus ik heb mijn oorspronkelijke short array (dus bevat waarden van ongeveer 32000 tot -32000) van bijvoorbeeld 1024 samples.

Nu kopieer ik dit in een nieuwe array geschikt voor de four1() functie. Het moet een float array zijn omdat de functie dit eist. Omdat ik alleen reele waarden heb, moet ik de imaginaire gedeeltes op 0 initialiseren.

Dus:

code:

float *nieuw_array = new float[1 + 2 *1024]; //reserveer ruimte (lengte = 2049)
for (int i = 0; i < 1024; i++)
{
    nieuw_array[1 + 2*i] = oud_array[i]; //kopieer sample in reele gedeelte
    nieuw_array[2 + 2*i] = 0; //zet imiginaire gedeelte op nul;
}

Nu gooi ik nieuw_array in de functie op deze manier:
four1(nieuw_array, 1024, 1);

En hiep hiep hoera, hij geeft geen errors, zowel in compile- als run-time.
Maar nu wordt mijn array dus getransformeerd in een fft-array. Tot zover lukt alles. Maar vervolgens weet ik niet hoe ik de waarden moet gebruiken. Een (co)sinusvorming signaal is te schrijven als een complex getal dus met 2 floats. De functie zorgt ervoor dat ik 1024 (in mijn geval) complexe getallen krijg.

Hoe interpreteer ik deze array?

[Tomatoman]

Een klein stukje theorie, zonder rekenwerk.

Als je een signaal f door een filter leidt, komt hier een ander signaal uit. Het uitgangssignaal het signaal is nogal ingewikkeld te berekenen, daarvoor moet je het signaal en het filter convolueren. Het convolutieteken is een sterretje, soms met een cirkeltje eromheen. Het uitgangssignaal ziet er als volgt uit:


Om je daar een beter beeld van te vormen kun je hier zien wat er gebeurt als je een blauw en een rood signaal convolueert, het resultaat is het groene signaal. Twee voorbeelden:



Aan de hand van de plaatjes kun je wel zien dat convolutie een ingewikkeld proces is dat erg veel rekenkracht kost. Daarom is er een prachtige oplossing gevonden, die gebruikt maakt van Fouriertransformaties. Als je de signalen f en g transformeert, kun je de ingewikkelde convolutie-bewerking vervangen door een eenvoudige vermenigvuldiging. Het resultaat transformeer je terug met de omgekeerde Fouriertransformatie en ziedaar: het uitgangssignaal.

Stel dat je
h =
wilt uitrekenen. De Fouriergetransformeerde van f is F, de Fouriergetransformeerde van g is G. H is de fouriergetransformeerde van H. Nu geldt:
H = F ×G
In de praktijk blijkt dat het minder (computer)rekenwerk kost om f en g te transformeren naar F en G, dan F en G te vermenigvuldigen (dit resulteert in H) en H terug te transformeren naar h.

Als je een tijdafhankelijk signaal hebt, is de Fouriergetransformeerde de representatie van datzelfde signaal in het frequentiespectrum. Daarom maken digitale signaalbewerkers (DSP's, digitale equalizers) ook vaak gebruik van Fouriertransformaties.


Tot zover de theorie. Bij AMD kun je gratis een bibliotheek met functies t.b.v. Fouriertransformaties downloaden. Voorbeelden in C zijn bijgevoegd in de download. Het is de AMD Core Math Library, kortweg ACML.

Ik heb de fourierfuncties laatst in Delphi getest en het resultaat is veelbelovend, de functies werken erg eenvoudig. Helaas loop ik tegen een vaag probleem aan dat de applicatie die de ACML aanroept altijd crasht bij afsluiten. Blijkbaar doe ik nog iets verkeerd, alleen weet ik niet wat.

Nog een aanvulling. De Fast Fourier Transform is een vorm van de Fouriertransformatie die veel efficiënter rekent. Voorwaarde in de meeste implementaties is dat het aantal samples een macht van twee is (16, 32, 64, 128...). Vaak ligt het maximale aantal samples bij 1024, 2048 of 4096. Jij wilt een kleine 65000 samples transformeren, dat is veel te veel. Ten eerste is dat waarschijnlijk onnodig voor het doel dat je wilt bereiken, ten tweede is het geen macht van 2 (hoewel je de array zou kunnen aanvullen met nullen), ten derde vinden veel algoritmes zoveel samples te veel en ten vierde is het onnodig (onwerkbaar?) rekenintensief.

[ Voor 17% gewijzigd door Tomatoman op 15-04-2004 02:30 ]

[.oisyn]

offtopic:
Waar lees jij dat ie 65k samples wilt? Ik lees 1024

[Verwijderd]

@afterdawn, om de uitkomst van de FFT te gebruiken moet je de amplitude berekenen.
Dit doe je door de re^2 + im^2 te nemen en dan de wortel. Als voorbeeld een stukje uit mijn code hieronder. Dit moet je nog wel aanpassen..ik maak gebruik van de C++ variant van NR, die dus wel start vanaf 0 en een vector+array klasse heeft om de library te gebruiken.

code:

void SpmData::FourierX(int &Row, int &ISign){
  NRVec<DP> Line(2*SizeX());                //storage for complex Fourier data

  for (int j=0; j<SizeX(); ++j) {
    Line[2*j] = Data2D[j][Row];
    Line[2*j+1] = 0;               //make imaginary part zero for real functions
  }
  NR::four1(Line, ISign);          //call the complex fourier, could be done by the realft!!
  //now store it in the original data
  for (int k=0; k<SizeX(); ++k) {
    Data2D[k][Row] = sqrt(SQR(Line[2*k])+SQR(Line[2*k+1]));
  }
}

Nu heb je een array met de amplitude. Bedenk dat er zowel iets aan de 0 kant zit als aan de N kant van de array door de manier waarop deze routine werkt.

@MsAlters: Uit pure frustratie over het gebruiken van de C code uit NR in C++ heb ik toen vrij snel de 1e druk van de C++ versie gekocht. Die is een stuk makkelijker in gebruik.
--> wat ik probeerde aan te geven is dat het compileren geen probleem is, maar het op de juiste wijze de data erin en eruit krijgen. De NR-libraries willen echt alles zo compact mogelijk doen. Niet altijd even logisch. Dus goed testen

ps) wat doe ik fout met de code-tag

[ Voor 15% gewijzigd door Verwijderd op 16-04-2004 18:24 ]

[Janoz]

Verwijderd schreef op 14 april 2004 @ 22:50:
Janoz, als je een perfecte sinus aanbied. Dus een sinus (die dus start op 0 in je data) is de FFT volledig imaginair (het reele deel is 0). Cosinus is precies omgekeerd.

Je hebt gelijk. Ik ben in de war. Ik had in gedachten het complexe getal al herschreven naar poolcoordinaten . Het is al weer een tijdje terug dat ik met het frequentie domein heb gewerkt, en die stap was ik even vergeten .

[Verwijderd]

@tomatoman:
Waar lees jij waar ik 65000 samples heb; ik heb er 1024, en mijn bereik is 16 bits unsigned dus van ongeveer -32000 tot 32000.


@moesasji1975nl:
Inderdaad, ik had al zo'n vermoeden dat ik de wortel van het reele gedeelte kwadraat plus het imaginaire gedeelte kwadraat moest nemen, maar wist het niet zeker. Ik ga het nu proberen. Dus als ik een array van 1024 punten heb (reeel + imaginair) dan heb ik dus 1024 frequentie banden met de bijbehorende amplitude.
Maar welk frequentie band hoort bij welke positie in de array. Op de NR site staat van oplopend van de kleinste frequentie tot de grootste negatieve frequntie ofzo.. hier wordt ik niet wijs van.

edit:
@markieman:
Ja klopt foutje van mij. Het bedoelde signed maar typte unsigned.

[ Voor 12% gewijzigd door Verwijderd op 15-04-2004 17:57 ]

[Markieman]

<sub>"... ik heb er 1024, en mijn bereik is 16 bits unsigned dus van ongeveer -32000 tot 32000. "

Mierenneuk: -32000 tot 32000 is volgens mij signed en niet unsigned</sub>

[Verwijderd]

De schaling is vrij simpel, je moet de 'echte' lengte van je oorspronkelijke array nemen.
Frequentie is dan de pixelwaarde gedeeld door die lengte. (kijk even naar de units, zit ik altijd mee te knoeien en test even met een bekende frequency voor factoren 2pi)

Je moet daarmee alleen de eerste helft van de array gebruiken.
Tweede helft bevat in principe zelfde info.

**) pas trouwens op. Slechts een beperkt gedeelte van de array is betrouwbaar (nyquist criterium). Verder moet je je realiseren dat je ook altijd de FFT van je window terug krijgt, in dit geval dus een FFT van een blok --> sync fynctie over je echte data.
Je kunt daarbij rare dingen kijken als je golven afgekapt worden door dat window.
--> als je hier meer over wil lezen, kijk dan eens naar bijvoorbeeld een Hamming window

[ Voor 6% gewijzigd door Verwijderd op 15-04-2004 13:31 ]

[Verwijderd]

Het werkt (bijna)!

Ik test mijn programma met een functiegenerator zodat ik een perfecte golf opneem. Hierdoor zie ik een mooie piek bij een (zelf in te stellen) frequentie.


Enige punt is nog dit:

Ik zie 2 pieken in mijn venster gespiegeld, nu kan ik dit wel oplossen door de halve array te nemen (zoals moesasji1975nl al zei).

(Ik sample trouwens op 44kHz.)
Nu zie ik nog 1 piek zoals het hoort. De piek begint mooi links en schuift gewoon naar rechts zodra ik de frequentie opschroef. Op een gegeven moment, als de piek op 11 kHz zit, en ik verhoog de frequentie nog meer, dan gaat de piek weer terug!! Tot 22kHz want dan is de piek weer aan de linkerkant!
(Bij een nog hogere frequentie zie ik niks meer / ruis. Dit is logisch, want ik weet dat op een samplerate van 44kHz geen frequenties boven de 22kHz kunnen worden gemeten.)


Waar mij het omgaat is dat na 11kHz de piek terug gaat in plaats van verder.
Ik heb 2 andere programma getest waarbij de ene programma (http://www.relisoft.com/freq.zip) inderdaad frequenties laat zien tot 11kHz (op 44kHz samplerate) en daarboven laat het programma NIETS zien ipv in mijn geval een gespiegelde golf.
Het andere programma ((http://techmind.org/audio/specan22.exe) laat daaraantegen wel tot 22kHz zien (waarbij de piek gewoon helemaal van links naar rechts gaat). Let op bij dit programma moet je wel via het menu de schaal even aanpassen.

Dus wat ik afvraag is of dit effect hierboven in mijn programma hoort.. of op een 1 of andere manier te omzeilen / oplossen valt. Ik wil in ieder geval proberen dat na 11 kHz mijn signaal niet meer terug gaat maar gewoon 'wegvalt'.

Wat is trouwens 'windowing effect' of dat 'sync'?

edit:

O ja 1024 samples werkt, net zoals 2048 en 4096 maar bij een hoger aantal samples crashed mijn programma (waarschijnlijk in de routine four1(...)).

[ Voor 6% gewijzigd door Verwijderd op 15-04-2004 18:31 ]

[Verwijderd]

Wat je hier ziet gebeuren is (volgens mij) een combinatie van effecten.

1) boven de 22kHz voldoe je niet langer voldoet aan het Nyquist criterium. Kort samengevat zegt dat criterium dat je voldoende punten op je oorspronkelijke golf moet hebben om een nauwkeurige FFT te kunnen doen. (volgens mij typisch meer dan drie punten per periode)

Dit kun je je vrij makkelijk voorstellen als je minder punten hebt kun je op een gegeven moment ook sinussen met een andere periode fitten aan de meetpunten. En dat is grofweg wat een FFT doet. Daardoor krijg je dus troep als resultaat voor te hoge frequenties.

Samengevat heb je dus hier boven een bepaalde gewoon niet genoeg samples. Zou je je aantal punten verhogen voor hetzelfde bereik, dan is een hogere frequentie haalbaar. Het verhogen van het aantal sample punten voor hetzelfde meetbereik is daarom ook altijd een cruciale test om te zien of een piek in het spectrum echt is.

2) Daarnaast heb je te maken met "windowing effecten". Dit heeft er mee te maken dat je dus niet alleen je oorspronkelijk golf transformeert, maar dus ook het window waarbinnen je meet. In dit window heb je data en daarbuiten niets. Daarom noemde ik dit al een blok-functie. Je transform is daardoor in het geval van het ideale geval een convolutie van je deltapiek met de FFT van een blok --> resultaat is dan een sync-functie.

ECHTER : je hebt hier niet het ideale geval. Je gebruikt waarschijnlijk sinussen met constante amplitude om te testen. Een FFT spiegelt wiskundig in feite je array van data-punten om een oneindige reeks te representeren. Als je dan dus sinussen hebt die geen geheel aantal keer in je window passen krijg je continu een fase-jump bij deze spiegeling. Effectief ziet je algoritme daardoor een heel andere periodiciteit (teken maar eens een array met steeds een halve periode gespiegeld) Het gevolg is dat je dan resultaten krijgt die dus nergens op slaan. Dit effect is vooral merkbaar als je heel weinig periodes in je sample overhoudt. Dus bij hoge frequenties. (maar dit treed op bij frequenties dieduidelijk nog beneden het nyquist criterium liggen)

Dat je de piek ziet teruggaan noem je terugvouwing en is dus een artefact van het ondersamplen/window effecten.

Oplossing is om je oorspronkelijke data-set te vermenigvuldigen met een andere window-functie zodanig dat de data op de randen van je array geleidelijk nul wordt. Dan heb je immers geen fase-jumps meer. Hierbij kun je kiezen voor veel verschillende filter-types wat een soort "zwarte kunst" is. Vaak volstaat vermenigvuldigen met een Gauss....daar heb je namelijk weinig last van omdat de FFT van een Gauss een Gauss oplevert.

ps) door gebruik te maken van een optimale window-functie heb je in principe wel voldoende aan het nyquist criterium van ten minste twee (of drie) punten per golf.
Maar in de praktijk wil je hier ver vandaan blijven. Hoe meer punten op je oorspronkelijke golf hoe beter

Als je het nu niet meer snapt gewoon doorvragen of in de boeken duiken.
Het heeft bij mij ook erg lang geduurd voordat ik het begreep.

ps) sync-functie is gewoon een wiskundige functie. Is moeilijk uit te leggen in woorden helaas. Waarschijnlijk geeft google je wel een plaatje.

ps2) het feit dat je de helft van de array niet nodig heeft is uit te leggen met het feit dat je positieve frequencies (eerste helft) en negatieve frequencies (tweede helft) hebt. Die negatieve frequencies horen dus eigenlijk voor de andere geplaatst te worden op de frequentie schaal. Dit is alleen relevant bij terug-transformaties en dat doe je hier niet. Dit is dus de reden dat je alleen naar 1 van de helften hoeft te kijken simpelweg omdat de fase van de FFT je niet interesseert.

[ Voor 42% gewijzigd door Verwijderd op 15-04-2004 20:29 ]

[Ralluph]

klok, klepel...

Verwijderd schreef op 15 april 2004 @ 19:51:
1) boven de 22kHz voldoe je niet langer voldoet aan het Nyquist criterium. Kort samengevat zegt dat criterium dat je voldoende punten op je oorspronkelijke golf moet hebben om een nauwkeurige FFT te kunnen doen. (volgens mij typisch meer dan drie punten per periode)

Het Nyquist-criterium zegt dat je een signaal moet sampelen met een frequentie die minimaal 2x zo hoog is als de maximale frequentie in het signaal om het nog eenduidig te kunnen reproduceren. Kijk eens wat google je hierover vertelt (eerste hit bijvoorbeeld).

(...)
2) Daarnaast heb je te maken met "windowing effecten". Dit heeft er mee te maken dat je dus niet alleen je oorspronkelijk golf transformeert, maar dus ook het window waarbinnen je meet. In dit window heb je data en daarbuiten niets. Daarom noemde ik dit al een blok-functie. Je transform is daardoor in het geval van het ideale geval een convolutie van je deltapiek met de FFT van een blok --> resultaat is dan een sync-functie.

Maak de TS niet gek. Hij heeft aangegeven dat hij er voldoende aan heeft als hij netjes een FFT kan implementeren. Alle standaard FFT-algoritmen gebruiken overigens een rechthoekig window.

ECHTER : je hebt hier niet het ideale geval. Je gebruikt waarschijnlijk sinussen met constante amplitude om te testen. Een FFT spiegelt wiskundig in feite je array van data-punten om een oneindige reeks te representeren. Als je dan dus sinussen hebt die geen geheel aantal keer in je window passen krijg je continu een fase-jump bij deze spiegeling. Effectief ziet je algoritme daardoor een heel andere periodiciteit (teken maar eens een array met steeds een halve periode gespiegeld) Het gevolg is dat je dan resultaten krijgt die dus nergens op slaan. Dit effect is vooral merkbaar als je heel weinig periodes in je sample overhoudt. Dus bij hoge frequenties. (maar dit treed op bij frequenties dieduidelijk nog beneden het nyquist criterium liggen)

Dat je de piek ziet teruggaan noem je terugvouwing en is dus een artefact van het ondersamplen/window effecten.

Je legt eerst uit dat het genoemde effect niet wordt veroorzaakt door terugvouwing (aliasing), en daarna ineens weer wel?

Oplossing is om je oorspronkelijke data-set te vermenigvuldigen met een andere window-functie zodanig dat de data op de randen van je array geleidelijk nul wordt. Dan heb je immers geen fase-jumps meer. Hierbij kun je kiezen voor veel verschillende filter-types wat een soort "zwarte kunst" is. Vaak volstaat vermenigvuldigen met een Gauss....daar heb je namelijk weinig last van omdat de FFT van een Gauss een Gauss oplevert.

Ik neem aan dat je met Gauss de klokvormige kromme bedoelt die de meesten hier wel zullen kennen als de kansdichtheidsfunctie van de normale verdeling. Als dat zo is: deze heeft hier helemaal niets mee te maken. Veel toegepaste windows zijn: Hanning, Hamming en flat-top.

Als je het nu niet meer snapt gewoon doorvragen of in de boeken duiken.
Het heeft bij mij ook erg lang geduurd voordat ik het begreep.

Ik zou nog eens doorlezen zou ik zeggen.

ps) sync-functie is gewoon een wiskundige functie. Is moeilijk uit te leggen in woorden helaas. Waarschijnlijk geeft google je wel een plaatje.

De sync-functie is niets meer en niets minder dan sin(x)/x.

ps2) het feit dat je de helft van de array niet nodig heeft is uit te leggen met het feit dat je positieve frequencies (eerste helft) en negatieve frequencies (tweede helft) hebt. Die negatieve frequencies horen dus eigenlijk voor de andere geplaatst te worden op de frequentie schaal. Dit is alleen relevant bij terug-transformaties en dat doe je hier niet. Dit is dus de reden dat je alleen naar 1 van de helften hoeft te kijken simpelweg omdat de fase van de FFT je niet interesseert.

De TS is slechts geïnteresseerd in amplituden en niet in de fasen. Hij wil vanuit een gegeven tijdrepresentatie een frequentierepresentatie van zijn signaal zien. Met terugtransformeren heeft hij niets te maken. De meeste FFT-algoritmen leveren overigens slechts het positieve spectrum.

@TS: ik heb hier nog wel wat blaadjes liggen met pseudocode voor diverse FFT-algoritmen. Ik heb echter geen zin en tijd om deze nu over te typen. Zoeken op google naar bijvoorbeeld [google]fft pseudocode OR "pseudo code"[/google] levert vast wat moois en begrijpelijks op.

[Verwijderd]

@Ralluph.....volgens mij moet jij nog eens HEEL goed lezen wat ik heb opgeschreven.
Je laat mij namelijk niet zien waarom je zegt klok...klepel. Erg veel van wat je opschrijft is namelijk hetzelfde, maar dan in andere woorden.

Het is duidelijk dat jij het kennelijk "perfect" snapt. (gezien een studie electro-techniek zou het moeten) Als dat zo is zou je ook begrijpen dat ik heb geprobeerd om het principe achter de FFT duidelijk te maken. Om aan de hand daarvan duidelijk te maken wat hij nu te zien krijgt uit de FFT-routine die hij gekozen heeft (dus een complexe FFT uit NR) De routine heeft TS overduidelijk al werkend namelijk, hij verbaast zich echter over de uitkomst. Die zijn echter inherent aan de logica van een FFT.

Dat je niet begrijpt dat een Gaussisch window (en dus inderdaad die klok-vormige kromme) erg goed bruikbaar is en heel makkelijk te gebruiken is als window zegt w.m.b. genoeg. (Inderdaad heb je volledig gelijk dan Hamming, flat-top e.d. veel gebruikte windows zijn, wat ik al aangaf in mijn eerdere posts in dit topic)

En ik heb inderdaad proberen uit te leggen dat hij te maken heeft met 2 effecten.
Als je even met de FFTs zou spelen die je hebt liggen zou je gauw genoeg zien wat ik bedoel. Probeer maar eens net voldoende sample punten te nemen volgens het Nyquist criterium terwijl je de perfecte sinus op een ongelukkige fase afkapt en dit zonder window-functie. (vergelijk het verschil tussen precies passende periode en ongelukkige fase, lijkt me duidelijk?)

ps) ik doe dit soort dingen dagelijks voor mijn werk veel doorlezen gaat niet meer.
ps2) TS vraagt wat windowing is.....mag ik het dan niet uit proberen te leggen
Zoniet leg jij dan maar even uit aan TS wat hij ziet tussen de 11kHz en de 22kHz.
Want DAT doe je dus ook niet (en hij sampled op 44kHz).

[ Voor 52% gewijzigd door Verwijderd op 15-04-2004 22:24 ]

[Ralluph]

Verwijderd schreef op 15 april 2004 @ 21:25:
@Ralluph.....volgens mij moet jij nog eens HEEL goed lezen wat ik heb opgeschreven.
Je laat mij namelijk niet zien waarom je zegt klok...klepel.

Het is duidelijk dat jij het kennelijk "perfect" snapt. (gezien een studie electro-techniek zou het moeten)

Allerminst, ik doe eigenlijk niets meer met deze materie (zie bijvoorbeeld mijn signature). Maar dat doet natuurlijk niets af aan het feit dat ik het nog zou moeten beheersen en vanuit die invalshoek geef ik graag mijn mening. Ik claim niet dat ik de onbetwiste expert ben op dit gebied, dat doe ik wel weer eens in een ander topic

Als dat zo is zou je ook begrijpen dat ik heb geprobeerd om het principe achter de FFT duidelijk te maken. Om aan de hand daarvan duidelijk te maken wat hij nu te zien krijgt uit de FFT-routine die hij gekozen heeft (dus een complexe FFT uit NR) De routine heeft TS overduidelijk al werkend namelijk, hij verbaast zich echter over de uitkomst. Die zijn echter inherent aan de logica van een FFT.

Jouw post was in zoverre duidelijk dat ik het op het eerste gezicht niet begreep. De TS geeft aan niet veel van de achterliggende mathematiek te (willen) begrijpen, dus ik kan mij voorstallen dat het hem helemaal duizelt. Mijn reactie was puur bedoeld als richtlijn voor de TS om weer wat duidelijkheid in de discussie te scheppen.

Dat je niet begrijpt dat een Gaussisch window (en dus inderdaad die klok-vormige kromme) erg goed bruikbaar is en heel makkelijk te gebruiken is als window zegt w.m.b. genoeg. (Inderdaad heb je volledig gelijk dan Hamming, flat-top e.d. veel gebruikte windows zijn.)

Het wat de boer niet kent dat eet hij niet-syndroom. Ik heb nog nooit een gaussische venster toegepast, vandaar dat ik niet wist dat het gebruikelijk is dit te begrijpen. Een typisch leermomentje, zie het PDF'je aan het einde van deze post. Wat het voordeel is om met een gaussische kromme te convolueren, boven andere vormen is mij overigens nog niet duidelijk.

edit:
Er is zelfs een MATLAB-functie voor: gausswin

En ik heb inderdaad proberen uit te leggen dat hij te maken heeft met 2 effecten.
Als je even met de FFTs zou spelen die je hebt liggen zou je gauw genoeg zien wat ik bedoel. Probeer maar eens voldoende sample punten te nemen terwijl je de perfecte sinus op een ongelukkige fase afkapt en dit zonder window-functie. (vergelijk het verschil tussen precies passende periode en ongelukkige fase, lijkt me duidelijk?)

ps) ik doe dit soort dingen dagelijks voor mijn werk veel doorlezen gaat niet meer.

Altijd blijven lezen hoor...

ps2) TS vraagt wat windowing is.....mag ik het dan niet uit proberen te leggen
Zoniet leg jij dan maar even uit aan TS wat hij ziet tussen de 11kHz en de 22kHz.
Want DAT doe je dus ook niet.

@TS: ik kan me voorstellen dat het zo langzamerhand begint te duizelen. Gelukkig levert google mij een documentje van onze vrienden bij Technische Natuurkunde waarin het een en ander over windowing staat uitgelegd. Wellicht wat te mathemtisch, maar er staan begrijpelijke plaatjes in. Te zien is dat de sinc-functie de Fourier getransformeerde is van de eenheidsblok-functie en vica versa. Er staan overigens ook voorbeelden in van een Gaussisch window.

[ Voor 4% gewijzigd door Ralluph op 15-04-2004 22:37 ]

[Verwijderd]

@Ralluph, complimented voor het gevonden documentje.....
hopelijk verhelderd het iets voor TS.
Plaatjes zeggen nu eenmaal meer dan duizend woorden.

specifiek voor Ralluph :
het voordeel van het gebruik van het Gauss window is dat je zowel in het tijdsdomein als frequency domein een gauss hebt. Vooral in de optica (waar ik het veel gebruik) is dit erg praktisch omdat het dan makkelijk mee te nemen is in de formules die toch al vol zitten met Gaussen.
(verder is het voor een fysicus snel te implementeren anders moet je ergens opzoeken hoe je het moet doen voor een Hamming window; wat beter en netter is overigens)

ps) lezen zal ik zeker blijven doen, maar na vier jaar kan ik geen FFT meer zien.

[ Voor 27% gewijzigd door Verwijderd op 15-04-2004 22:50 ]

[Verwijderd]

Sorry maar snap er niet echt wat van...

Maakt ook niet uit want in principe werkt het programma. Dit is gewoon een 'effect' waar de gebruiker maar rekening mee moet houden .

En ik vraag het ook nog even aan mijn leraar. In ieder geval bedankt!

edit:
ff wachten zal nog ff dat pdfje doorspitten

[ Voor 13% gewijzigd door Verwijderd op 15-04-2004 23:09 ]

[Verwijderd]

oke ik heb het gelezen. ik denk dat mijn probleem niet echt ligt in dat 'windowing effect'. Ik heb weliswaar ruis (wat veroorzaakt wordt door het windowing effect; want ik sample een perfecte golf met een constante frequency dmv de functie-generator).
Volgens mij zie ik gewoon net iets over het hoofd, of maak ik een logisch foutje ofzo...

K zal nog even laten zien hoe ik mijn fft met bijvoorbeeld 2048 punten (dus 2048 reeel + 2048 imaginair) precies laat zien:

code:

//dc is de device-context van mijn tekenvlakje,
//voor het gemak is mijn tekenvlak even 1024 pixels breed en 300 pixels hoog.

//complexSignaal bevat mijn signaal

int max = complexSignaal.getNumberOfComplexPoints / 2; //1024 dus (past precies op mijn vlak dus)
int hoogte_aanpassing = complexSignaal.getGrootsteAmplitude() / 300;

dc.MoveTo(0, 300) //begin linksonder
for (int x = 0; x < max; x++)
{
    //functie getAmplitude berekent de amplitude dmv wortel(reeel^2 + imaginair^2);
    int y = 300 - complexSignaal.getAmplitude(x) / hoogte_aanpassing;
    dc.LineTo(x, y);
}

edit:

Oja dit deze code laat dus 1 piek zien.

[ Voor 15% gewijzigd door Verwijderd op 15-04-2004 23:54 ]

[Verwijderd]

Aan deze code zie ik niets vreemds. Dit is puur display.

Bij het testen heb je een perfecte golf met een bepaalde frequentie.
Kies die frequentie even zodanig dat je een geheel aantal periodes in de lengte van je window past. Kijk dan eens of je nog steeds vreemde dingen ziet.

ps) even voor alle zekerheid : je hebt wel steeds alternerend reeel, imaginair, reeel, ...in je arrays gebruikt

[ Voor 18% gewijzigd door Verwijderd op 16-04-2004 00:01 ]

[Verwijderd]

Ik heb nu geen functie-generator tot mijn beschikken, alleen op school. Maar omdat het bij mij gaat om frequenties boven de 11kHz (dus waarbij de piek weer terugloopt), en daarbij de geluidsgolf op het superwazig is, kan ik nooit gehele periodes in mijn sample-array hebben.

Nog even voor de duidelijkheid het gaat me dus niet om de ruis; het principe van de fft is gelukt. Het gaat me om het feit dat boven de 11kHz de piek terugloopt.

[Verwijderd]

Maar waarom test je dit dan niet met een array die je zelf aanmaakt in je programma?

Dus gewoon een array vullen met sin(x) en die dan in je FFT stoppen.
Dan weet je precies wat je erin stopt en kun je er ook makkelijk mee spelen.

[Verwijderd]

Verwijderd schreef op 15 april 2004 @ 23:53:

ps) even voor alle zekerheid : je hebt wel steeds alternerend reeel, imaginair, reeel, ...in je arrays gebruikt

Ja hoor als je dit bedoelt: getAmplitude(x) berekent in feite dit:

code:

float getAmplitude(int point)
{
    return sqrt( pow(data[1 + 2*point],  2) + pow(data[2 + 2*point], 2) );
}

[Verwijderd]

Perfect, dit is precies zoals het moet zijn.
(Dit is echter precies zo'n stommiteit om over het hoofd te zien dus vandaar dat ik het nog even vroeg voor alle zekerheid)

[ Voor 32% gewijzigd door Verwijderd op 16-04-2004 00:36 ]

[Verwijderd]

Ok nou dan ga ik morgen (denk ik) expirimenteren met dat handmatig een sinus in de array stoppen dmv sin(). Ik houd jullie nog op de hoogte

[Zoijar]

Verwijderd schreef op 16 april 2004 @ 00:13:
[...]


Ja hoor als je dit bedoelt: getAmplitude(x) berekent in feite dit:

code:

float getAmplitude(int point)
{
    return sqrt( pow(data[1 + 2*point],  2) + pow(data[2 + 2*point], 2) );
}

Amp is gewoon de lengte van de vector in de 2d re/im lineaire ruimte. Maar wat ik eigenlijk wilde zeggen...pow(x, 2)... Vermenigvuldig gewoon met hetzelfde, waarschijnlijk een stuk sneller...

[Verwijderd]

Zoijar schreef op 16 april 2004 @ 03:28:
[...]

Amp is gewoon de lengte van de vector in de 2d re/im lineaire ruimte. Maar wat ik eigenlijk wilde zeggen...pow(x, 2)... Vermenigvuldig gewoon met hetzelfde, waarschijnlijk een stuk sneller...

De lengte van de vector wordt toch berekend (dmv pythagoras)?

Oja en over dat pow(x, 2)... is er geen kwadraat functie?

[MSalters]

Ja, X * X

[Verwijderd]

Dat klopt! Maar ik bedoelde eigenlijk gewoon een standaard functie, maar goed dan doe ik dit wel


@moesasji1975nl:
Regel 10 van je eigen code die je liet zien, daarvan neem ik aan dat k niet goed is. Dit moet toch j zijn?

edit:
Oh nee ik zie het verkeerd, de for-lus wordt al afgesloten. Dan klopt trouwens het aantal brackets niet.

edit 2:
Nu zie ik het echt duidelijk ; het is een stukje van je code dus er mist nog wat.. bijv. een for(...){ of iets dergelijks.

[ Voor 38% gewijzigd door Verwijderd op 16-04-2004 17:58 ]

[Verwijderd]

Inderdaad had ik in mijn oorspronkelijke code iets te enthousiast niet terzake doende dingen gedelete. Daardoor was het begin van een for-loop weggevallen.
Nu is het wel goed (even aangepast)

Sorry, niet goed opgelet.

[Verwijderd]

Oke, tijdje geleden weer maar hier zijn resultaten:

Dit is code om een array te vullen met een sinus, hiermee heb ik mijn FFT-algoritme getest:

code:

const int numberSamples = 4096;
short data[numberSamples];
const double amplitude = 10000;
const double PI = 3.14159265358979323846;
const double w = 0.4; //aan te passen voor gewenste frequente

for (int i = 0; i < numberSamples; i++)
{
    data[i] = (int) (amplitude * sin(w * PI * i);)
}

Oja deze array wordt uiteraard nogwel geconverteerd naar een geschikte array voor de FFT.

Doordat ik w * PI doe krijg ik altijd hele periodes (dit klopt ook; want dit zie ik op mijn oscilloscoop ) en door w te variëren kan ik de frequentie aanpassen.

Nu laat ik w variëren van 0.001 tot 0.5 waarbij de piek gewoon van links tot rechts verschuift. Boven de 0.5 loopt de piek gewoon weer terug.

Mijn conclusie hieruit is dat de functiegenerator en de geluidskaart geen invloed hebben op mijn "probleem" (dat een aantal posts hierboven wordt beschreven).

[Verwijderd]

Mijn eerste reaktie als ik deze code zie....

waarom gebruik je het data-type short voor je data
Een data-type double (precision) lijkt me beter op zijn plaats voor iets wat je in een FFT wil stoppen.

ps) typecasting van variabelen, probeer je dit niet aan te leren in de vorm (int).
Ik weet dat het makkelijk is, maar het is niet voor niets afgeschaft
(gebruik het in de vorm van static_cast<int> e.d. dan zie je in een oogopslag dat er een cast zit)

ps2) verder nu nog even geen nuttige suggesties, kom er nog wel op terug.
Nu geen tijd en energie om hierover na te denken.

[ Voor 10% gewijzigd door Verwijderd op 21-04-2004 11:28 ]

[Verwijderd]

Over dat ik de data in een short heb: Er staat in mijn post dat ik de data in een geschikte array kopieer (die idd gewoon float is). En dat converteren zal ik even aanpassen.

[Verwijderd]

Oh nee k zie al wat je bedoelt: Om met een sinus de fft te testen, kan ik beter de sinus diréct al in een float-array stoppen. Daar had ik even niet aan gedacht.

[Verwijderd]

Woei woei!

Hij doet het! Eindelijk werkt het nu naar behoren. Ik heb tijdens het kopieren van de array in een geschikte array voor de fft-functie de sample-waarde gekopieerd met een Hamming-functie.

Nu zie ik gewoon als ik de frequentie opschroef mijn piek van links naar rechts gaan (in het bereik van 0 tot 22kHz!!! )

Oké iedereen bedankt voor zijn aandacht en hulp.

edit:
Kleine correctie: Ik heb dat met dat hamming-window eens weg gehaald..en toen deed hij het alsnog goed!!! Wat blijkt.. ik had voor dat hamming-gebeuren de code voor het kopieren van de array is een geschikte fft-array herschreven. Daar zat blijkbaar een fout in maar nu niet meer.

[ Voor 28% gewijzigd door Verwijderd op 28-04-2004 18:08 ]