[C++] onduidelijkheden over dataformaat FFTW

[Arjan]

Niet gehinderd door wiskundige kennis was ik wat aan het experimenteren met DFT filters. Na wat gelezen te hebben besloot ik om de FFTW library te gebruiken voor het berekenen van de DFT en de inverse.

In principe ging dat voorspoedig, ik kan een geluid inlezen, naar freq. domein plotten en weer terug rekenen naar tijd domein.

Nu wilde ik echter wat gaan filteren en toen werd het wat lastiger.

FFTW slaat de data op in een array van complexe getallen, bestaande uit twee doubles.
de eerste is het reële gedeelte, de tweede het imaginaire als ik het goed heb.

aanvankelijk dacht ik dan gewoon een tweedimensionale array te krijgen waarbij [i][0] naar de reële- en [i][1] naar de imaginaire data wees.
Maar dit klopt maar deels, want ik heb inmiddels begrepen dat als je een array grootte van 2048 opgeeft, er slechts de helft van wordt gebruikt voor reële data

De data wordt gespiegeld op de helft, het nut hiervan ontgaat me, maar ik ben er dus wel in geslaagd om op die manier een high/low pass filter te maken door de eerste/laatste x kolommen te bewaren en de rest op 0 te zetten.

Ik snap alleen niet helemaal wat er nu gebeurd. Waar is al deze extra data voor nodig? en als ik wil filteren, moet ik dan ook het imaginaire gedeelte bewerken?

hopelijk kan iemand er wat licht op schijnen, want ik mis de benodigde wiskunde om dit op eigen houtje uit te vogelen

[FragFrog]

In principe krijg je een reele waarde uit een imaginair getal alleen correct als je de wortel neemt van het kwadraat van het imaginaire deel plus het kwadraat van het reele deel, als ik m'n colleges complexe analyse tenminste goed onthouden heb. In die zin kan het nuttig zijn dat je ook het imaginaire gedeelte meekrijgt en moet je dus niet alleen het reele gedeelte pakken voor de eindwaarde.

In hoeverre dit van toepassing is op deze filters kan ik je helaas niet zeggen though

[Verwijderd]

/me is ook niet gehinderd door veel wiskundige kennis dus alles onder voorbehoud.

je haalt 2 problemen door elkaar. dat complexe getal is gewoon wat eruit hoort te komen. je reëele waarde is gewoon het piefje in je frequentiespectrum, en de imaginaire waarde de fase.
wat betreft die spiegeling: dat zou wel eens doodordinaire aliasing kunnen zijn. dat zijn frequenties hoger dan de helft van je samplefrequentie die voor lijken te komen, maar in feite helemaal niet voorkomen (en die moet je idd wegfilteren )

[Arjan]

volgens mij is het spiegelen wel by design, de ene helft is de sinus en de tweede helft cosinus, om zo ook de fase van het signaal te kunnen berekenen.

wat ik daarbij echter niet snap is wat de imaginaire data dan nog toevoegt aan het geheel...

als ik de imaginaire data echter negeer bij het filteren dan lijkt het signaal amper aangepast bij het terugrekenen

[MSalters]

Spiegelen is vrij normaal in Fourier Transformaties. Uiteindelijk komt dat doordat cos(x)=cos(-x); die spiegelt ook.

Je array (van bijvoorbeeld 2048 punten) moet zowel het reele deel als het imaginaire gaan bevatten. Aangezien bij elke frequentie een complexe amplitude hoort is het logisch dat dan maar 1024 punten uit de array je reele delen zijn.

Je imaginaire delen bevatten fase-informatie voor die frequentie. Voor de frequenties die je bewaart moet je ze bewaren. De fase van een "golf" met amplitude 0 is irrelevant, dus die mag je bewaren of op 0 stellen, dat maakt niet meer uit: 0 * sin(x) = 0 * sin(x+c)

[Arjan]

Maar mag ik deze fase informatie dan gewoon behandelen zoals ik de amplitude behandel bij het filteren?

bijvoorbeeld dit low pas achtige ding dat ik even snel opgeschreven heb:

C++:

void doeiets() {

    if(filterToggle) return;

    float a = 0.1;

    for(int i=20; i<FFT_SIZE/2; i++) {
        
        out[i][0] *= a;
        out[i][1] *= a;
    }

    for(int i=FFT_SIZE/2; i<FFT_SIZE-20; i++) {

        out[i][0] *= a;
        out[i][1] *= a;
    }
}

C++:

out[i][1]

is dus het imaginaire deel, het 'filter' werkt wel, maar ik hoor tikjes in het geluid bij het afspelen, vandaar dat ik me afvraag of ik de bewerkingen correct doe

[Arjan]

hier de glitch visueel:

[Soultaker]

Arjan schreef op zondag 09 december 2007 @ 15:16:
wat ik daarbij echter niet snap is wat de imaginaire data dan nog toevoegt aan het geheel...

als ik de imaginaire data echter negeer bij het filteren dan lijkt het signaal amper aangepast bij het terugrekenen

Misschien wist je dat al, maar het klopt dat je voor geluid meestal de imaginaire as kunt negeren. Het menselijk oor is alleen gevoellig voor de sterkte van het signaal op verschillende frequenties (dat is de reële component) en niet voor de fase (dat is de imaginaire component). Als je de imaginaire component dus negeert zul je in de praktijk geen verschil horen, maar de imaginaire component is wél nodig om de exacte invoer te reconstrueren (voor zover rekenen met floating point getallen exact is, natuurlijk). Verder kun je door interferentie e.d. nog wel wat verschillen krijgen, maar daar heb je waarschijnlijk geen last van als je niet met een synthetische signaal begint.

De imaginaire component schalen zoals je in je voorbeeld doet lijkt me ook niet zo zinnig. Wat denk je daarmee te bereiken? Ik ken verder FFTW ook niet echt, maar heb je niet gewoon de dimensies van je array omgedraaid toevallig, waardoor je dus eerst 2048 reële waarden krijgt, en dan de 2048 imaginaire waarden?

[.oisyn]

Soultaker schreef op zondag 09 december 2007 @ 16:46:
[...]

Misschien wist je dat al, maar het klopt dat je voor geluid meestal de imaginaire as kunt negeren. Het menselijk oor is alleen gevoellig voor de sterkte van het signaal op verschillende frequenties (dat is de reële component) en niet voor de fase (dat is de imaginaire component).

Ho ho, dit klopt niet. De amplitude is |x| en de fase is atan2(x.imag, x.real). Als je de vectors converteert naar polaire vorm dan heb je idd de amplitude (de magnitude) en de fase (de hoek) los van elkaar, maar dat komt niet overeen met resp. het reële en het imaginaire deel van het complexe getal, want die zijn Euclidisch.

[ Voor 8% gewijzigd door .oisyn op 09-12-2007 18:07 ]

[Arjan]

Soultaker schreef op zondag 09 december 2007 @ 16:46:
[...]

De imaginaire component schalen zoals je in je voorbeeld doet lijkt me ook niet zo zinnig. Wat denk je daarmee te bereiken? Ik ken verder FFTW ook niet echt, maar heb je niet gewoon de dimensies van je array omgedraaid toevallig, waardoor je dus eerst 2048 reële waarden krijgt, en dan de 2048 imaginaire waarden?

als ik die imaginaire component niet schaal dan gaat m'n geluid haywire bij de inverse...
omdraaien van de array lijkt me ook niet, gezien de manier waarop je deze aanmaakt:

C:

         fftw_complex *in, *out;
         fftw_plan p;
         ...
         in = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * N);
         out = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * N);
         p = fftw_plan_dft_1d(N, in, out, FFTW_FORWARD, FFTW_ESTIMATE);
         ...
         fftw_execute(p); /* repeat as needed */
         ...
         fftw_destroy_plan(p);
         fftw_free(in); fftw_free(out);

fftw_complex == double[2]

ik gok dat mijn geluid-glitch voorkomt uit het gebrek aan windowing, daar ga ik eens even mee experimenteren.

[Soultaker]

.oisyn schreef op zondag 09 december 2007 @ 18:00:
Ho ho, dit klopt niet. De amplitude is |x| en de fase is atan2(x.imag, x.real). Als je de vectors converteert naar polaire vorm dan heb je idd de amplitude (de magnitude) en de fase (de hoek) los van elkaar, maar dat komt niet overeen met resp. het reële en het imaginaire deel van het complexe getal, want die zijn Euclidisch.

Oh ja, zo zat het. Ik haalde wat dingen door elkaar. Dan klopt het schalen ook gewoon.

Arjan schreef op zondag 09 december 2007 @ 18:44:
als ik die imaginaire component niet schaal dan gaat m'n geluid haywire bij de inverse...
omdraaien van de array lijkt me ook niet, gezien de manier waarop je deze aanmaakt

Ah, het was maar een suggestie. Is achteraf gezien ook logischer om die componenten bij elkaar te houden (efficienter uitrekenen).

[Arjan]

Na wat testen blijkt dat het probleem ontstaat door mijn 'filter'.
Hiermee kom ik weer opnieuw terug bij de vraag: hoe moet je de data bewerken om correct te filteren?

Het lijkt erop dat ik niet zomaar afzonderlijke frequenties kan versterken/verzwakken. (opzich ook niet zo raar natuurlijk aangezien je dan de fase van deze frequenties afzonderlijk van de rest aanpast)

[.oisyn]

Als je zowel het imaginaire deel als het reële deel met dezelfde waarde vermenigvuldigd dan blijft de fase juist gelijk. Hoe vul je eigenlijk de buffer, en hoe ga je weer terug naar het tijddomein? Denk je er wel aan om bij de input alle imaginaire delen op 0 te zetten? En de output te vermenigvuldigen met 1/N, en voor elke sample |x| te pakken ipv x.real?

[ Voor 6% gewijzigd door .oisyn op 09-12-2007 23:57 ]

[Arjan]

.oisyn schreef op zondag 09 december 2007 @ 23:57:
Als je zowel het imaginaire deel als het reële deel met dezelfde waarde vermenigvuldigd dan blijft de fase juist gelijk. Hoe vul je eigenlijk de buffer, en hoe ga je weer terug naar het tijddomein?

C++:

for(i=0; i<FFT_SIZE; i++) {
    in[i][0] = in[i][1] = input[i];
}
...
toFreq = fftw_plan_dft_1d(FFT_SIZE, in, intermediate, FFTW_FORWARD, FFTW_MEASURE);
toTime = fftw_plan_dft_1d(FFT_SIZE, intermediate, out, FFTW_BACKWARD, FFTW_MEASURE);

ik schrijf dus de float audio-data uit input naar in, roep

code:

fftw_execute(toFreq);

aan, hetgeen me het freq. domein oplevert in intermediate.
als ik vervolgens

code:

fftw_execute(toTime);

aanroep, dan heb ik in out weer m'n audio data.

Denk je er wel aan om bij de input alle imaginaire delen op 0 te zetten?

doe ik atm niet, FFTW heeft het er ook niet over in de docs. Zal testen of het wat uitmaakt.

En de output te vermenigvuldigen met 1/N, en voor elke sample |x| te pakken ipv x.real?

ik normalizeer de output idd, maar heb geen idee wat je met |x| bedoelt, (absolute waarde van een complex getal ) maar atm. pak ik de reële waarde van de sample.

[farlane]

Lengte van x neem ik aan ...

[Arjan]

waar hebben we het dan over, over het terugschrijven naar bruikbare audio data, of bij het plotten van het freq. domein naar een frequentiespectrum?
Bij dat laatste gebruik ik idd. de lengte, bij de eerste de reële data.

[.oisyn]

C++:

for(i=0; i<FFT_SIZE; i++) {
    in[i][0] = in[i][1] = input[i];
}

Da's dus niet goed. Je input is in het reële domein, dwz een complex getal zonder imaginaire. in[i][1] moet dus altijd 0 zijn (nou ja, technisch gezien maakt het niet zoveel uit, als de amplitude maar klopt, oftewel √(in[i][0]² + in[i][1]²) == input[i], maar het is wel practisch om het imaginaire gedeelte gewoon op 0 te zetten en de input te kopiëren naar het reële gedeelte).

ik normalizeer de output idd, maar heb geen idee wat je met |x| bedoelt, (absolute waarde van een complex getal ) maar atm. pak ik de reële waarde van de sample.

Wat ik al eerder zei, de amplitude van een sample is gelijk aan |x|, en daarin ben je geïnteresseerd. Je kunt een complex getal voorstellen als een 2d vector. De amplitude is de lengte van de vector, en de fase is de hoek die de vector maakt met (1, 0) (oftewel 1 + 0i, oftewel 1). |x| is idd het absolute (of de magnitude of de modulus) van een complex getal x, en is te berekenen door √(x.real² + x.imag²). Ik weet trouwens even zo snel niet of je filter ervoor zorgt dat het imaginaire gedeelte ooit niet 0 zou zijn, zo niet dan gaat het dus gewoon goed, maar je kunt beter het zekere voor het onzekere nemen

Je kunt je output idd normalizeren, maar om exact eruit te krijgen wat je erin stopte (als je een forward FFT direct gevolgd door een backward FFT zou doen) moet je vermenigvuldigen met 1/N.

[ Voor 12% gewijzigd door .oisyn op 10-12-2007 02:08 ]

[Arjan]

hmm, ok, dan snap ik niet hoe ik een normale sinus terugvorm...
√(x.real2 + x.imag2) levert toch alleen maar positieve getallen op?

[.oisyn]

Hmm, good point, helemaal vergeten. Kun je niet beter gewoon een complex to real conversie gebruiken in FFTW?

[Arjan]

een complex to real conversie?

wat ik niet helemaal snap is de reden dat er tikjes in m'n output komen. Zoals al aangegeven was het gelukt om van/naar freq/tijd domein te converteren.

Als ik alle frequenties met dezelfde waarde schaal dan heb ik effectief een enorm ingewikkelde volume setting, maar zodra ik één specifieke frequentie schaal gaat de output tikken omdat de afzonderlijk geconverteerde gedeelten niet goed op elkaar aansluiten...

[voodooless]

Even tussendoor glippen: Al deze dingen heb je helemaal niet nodig eigenlijk Zoek eens op hoe je FIR filters maakt. Effectief doen die hetzelfde als je FFT -> filter -> invFFT, alleen veel sneller, en veel simpeler Nouja, simpeler in de zin van berekeningen uitvoeren. De filter coefficienten berekenen is weer een complexer verhaal

Neemt natuurlijk niet weg dat het zeer zeker interessant is om het op een andere manier te doen

[ Voor 34% gewijzigd door voodooless op 11-12-2007 12:45 ]

[.oisyn]

Arjan schreef op dinsdag 11 december 2007 @ 12:01:
een complex to real conversie?

Check de documentatie van FFTW. Je input en output zijn reals (reële getallen), maar de FFT die je nu gebruikt werkt met complexes. FFTW ondersteunt ook gewoon real to complex FFTs en weer terug.

wat ik niet helemaal snap is de reden dat er tikjes in m'n output komen. Zoals al aangegeven was het gelukt om van/naar freq/tijd domein te converteren.

Als je je reals ook naar de imaginairen kopieert dan klopt het gewoon niet en is het logisch dat je eindresultaat ook niet klopt. Overigens moet je ook eens bijlezen over circular convolution, want met een goede conversie ben je er nog niet. Je deelt je stream op in blokjes, die voor het doel van de FFT worden behandeld alsof ze periodiek zijn. Het is dus niet raar dat je na filtering clicks krijgt. Je zou bijvoorbeeld een grotere window kunnen gebruiken die je steeds maar voor de helft opschuift, en daarvan alleen de middelste helft gebruikt (evt. iets groter dan de helft en crossfading)

Overigens zijn FIRs en IIRs idd een stuk simpeler

[Arjan]

.oisyn schreef op dinsdag 11 december 2007 @ 12:59:
[...]

Check de documentatie van FFTW. Je input en output zijn reals (reële getallen), maar de FFT die je nu gebruikt werkt met complexes. FFTW ondersteunt ook gewoon real to complex FFTs en weer terug.


[...]

Even tussendoor glippen: Al deze dingen heb je helemaal niet nodig eigenlijk Zoek eens op hoe je FIR filters maakt. Effectief doen die hetzelfde als je FFT -> filter -> invFFT, alleen veel sneller, en veel simpeler Nouja, simpeler in de zin van berekeningen uitvoeren. De filter coefficienten berekenen is weer een complexer verhaal

Neemt natuurlijk niet weg dat het zeer zeker interessant is om het op een andere manier te doen

De insteek van deze oefening was enerzijds om een werkend filter te krijgen, maar belangrijker, om wat informatie over DFT/FFT bewerkingen op te doen
Maar ik zal zeker nog eens kijken naar FIR filters, zeker als het efficiënter is, dan zou ik het kunnen gebruiken op m'n gp2x.

.oisyn schreef op dinsdag 11 december 2007 @ 12:59:
[...]

Check de documentatie van FFTW. Je input en output zijn reals (reële getallen), maar de FFT die je nu gebruikt werkt met complexes. FFTW ondersteunt ook gewoon real to complex FFTs en weer terug.


[...]

Als je je reals ook naar de imaginairen kopieert dan klopt het gewoon niet en is het logisch dat je eindresultaat ook niet klopt. Overigens moet je ook eens bijlezen over circular convolution, want met een goede conversie ben je er nog niet. Je deelt je stream op in blokjes, die voor het doel van de FFT worden behandeld alsof ze periodiek zijn. Het is dus niet raar dat je na filtering clicks krijgt. Je zou bijvoorbeeld een grotere window kunnen gebruiken die je steeds maar voor de helft opschuift, en daarvan alleen de middelste helft gebruikt (evt. iets groter dan de helft en crossfading)

Overigens zijn FIRs en IIRs idd een stuk simpeler

Ah, dan zat ik toch goed met mijn vermoeden dat het probleem bij windowing lag. Ga ik vanavond eens proberen met alleen het middenstuk van het bewerkte gedeelte te gebruiken

[Arjan]

ok, nu ik slechts het middenstuk gebruik werkt het perfect!
hier de code voor een ieder die ook met dit speelgoed aan de slag wil:
het FFTW gedeelte

C++:

    in = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * FFT_SIZE);
    out = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * FFT_SIZE);
    intermediate = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * FFT_SIZE);
    toFreq = fftw_plan_dft_1d(FFT_SIZE, in, intermediate, FFTW_FORWARD, FFTW_MEASURE);
    toTime = fftw_plan_dft_1d(FFT_SIZE, intermediate, out, FFTW_BACKWARD, FFTW_MEASURE);

de SDL audio callback

C++:

void audioCallback(void *data, Uint8 *buffer, int size) {

    int samples = size / 2;
    int p = 0;
    int offset = 0;
    int i;
    int start = FFT_SIZE / 4;
    int end = FFT_SIZE - start;

    double maxVolume = (double)0x7FFF;

    while( samples ) {

        offset = 0;
        for(i=0; i<FFT_SIZE; i++) {

            if(pos+i-offset > totalSamples) offset = pos+i;

            in[i][0] = input[pos+i-offset];
            in[i][1] = 0.0;
        }

        fftw_execute(toFreq);

        // doe hier het filteren in frequency domein

        fftw_execute(toTime);

        for(i=0; i<FFT_SIZE; i++) {
            out[i][0] /= FFT_SIZE;
            out[i][1] /= FFT_SIZE;
            if(out[i][0] > 1.0) {
                out[i][0] = 1.0;
            } else if(out[i][0] < -1.0) {
                out[i][0] = -1.0;
            }
            if(out[i][1] > 1.0) {
                out[i][1] = 1.0;
            } else if(out[i][1] < -1.0) {
                out[i][1] = -1.0;
            }
        }

        for(i=start; i<end && samples; i++) {

            ((int16_t*)buffer)[p] = (int16_t)( out[i][0] * maxVolume );

            p++;
            if(pos < totalSamples) {
                pos++;
            } else {
                pos = 0;
            }
            samples--;
        }
    }

    return;
}