[C++] Een variatie op de bit-reversal methode voor de FFT

[Verwijderd]

Ik ben al een tijdje aan het stoeien met het volgende probleem. Zoals jullie misschien wel weten wordt in het FFT algoritme gebruik gemaakt van 'bit reverse order'. Dit houdt in dat de elementen van een array herschikt worden volgens de bit-reversal methode.

Bijvoorbeeld.
We hebben een array van lengte 2⁴=16. Het element op index 3 (=0011) komt op index 12 (=1100), het element op index 7 (=0111) komt op index 14 (=1110), het element op index 1 (=0001) komt op index 8 (=1000) enzovoort..

Nu wil ik graag een loopje hebben die dit voor mij doet maar met een kleine variatie. In plaats van dat we kijken naar *alle* bits in de indexen wil ik alleen de bits omdraaien die op plek p t/m q staan.

Bijvoorbeeld voor p=1 en q=2 worden de volgende indexen verwisselt.
(we tellen de bits van rechts naar links, beginnende by nul)

C++:

3210    3210
============
0000 -> 0000 (hier gebeurt dus niks, maar staat hier voor de volledigheid)
0010 -> 0100
0100 -> 0010
0110 -> 0110 (gebeurt ook niks)

De rest van de indexen blijven ongemoeid.


Of voor p=0 en q=2:

C++:

3210    3210
============
0000 -> 0000 (hier gebeurt dus niks, maar staat hier voor de volledigheid)
0001 -> 0100
0010 -> 0010 (gebeurt niks)
0011 -> 0110
0100 -> 0001
0101 -> 0101 (gebeurt niks)
0110 -> 0011
0111 -> 0111 (gebeurt niks)

Bij voorbaat hartelijk dank!

[ Voor 4% gewijzigd door Verwijderd op 15-06-2006 23:04 ]

[Invisible_man]

Om te beginnen zou je de bitposities die je niet wilt verplaatsen er kunnen filteren doormiddel van masken. Die plaats je dan alvast in de nieuwe array.

Vervolgens filter je met wederom een masker (maar dan geinverteerd aan de oude) de bits uit de oude array.

Stel ik wil p=0 tot q=2, dan wordt hou je na het masken van bijv. 1110, 110 over (6 decimaal). De meest signifikante bit heeft een waarde van 4, dus je begint met te delen door 4. Levert dit (integer) 1 op, weet je dat die bit 1 is, levert dit 0 op, is die bit ook 0. Msb wordt lsb, dus +1 naar de eind-array schrijven in geval dat bitje 1 is en de waarde 4(decimaal) van 101(binair) aftrekken, wordt 01(binair). Dit doe je ook voor msb-1 maar dan met de waarde twee. Na het aftrekken van 4 hield je in het voorbeeld 2 over. 2 / 2 is 1, dus ook bit msb-1 is een 1(binair). Msb-1 wordt lsb+1, dus +2 bij het eindresultaat en 2 weer aftrekken, hou je 0 over. 0/1 = 0, dus msb-2 is een 0, dus 0*lsb+2. Enzovoorts, enzovoorts.

[ Voor 4% gewijzigd door Invisible_man op 15-06-2006 23:19 ]

[Peter_B]

Verwijderd schreef op donderdag 15 juni 2006 @ 23:01:
Ik ben al een tijdje aan het stoeien met het volgende probleem. Zoals jullie misschien wel weten wordt in het FFT algoritme gebruik gemaakt van 'bit reverse order'. Dit houdt in dat de elementen van een array herschikt worden volgens de bit-reversal methode.

Bijvoorbeeld.
We hebben een array van lengte 2⁴=16. Het element op index 3 (=0011) komt op index 12 (=1100), het element op index 7 (=0111) komt op index 14 (=1110), het element op index 1 (=0001) komt op index 8 (=1000) enzovoort..

Nu wil ik graag een loopje hebben die dit voor mij doet maar met een kleine variatie. In plaats van dat we kijken naar *alle* bits in de indexen wil ik alleen de bits omdraaien die op plek p t/m q staan.

Bijvoorbeeld voor p=1 en q=2 worden de volgende indexen verwisselt.
(we tellen de bits van rechts naar links, beginnende by nul)

C++:

3210    3210
============
0000 -> 0000 (hier gebeurt dus niks, maar staat hier voor de volledigheid)
0010 -> 0100
0100 -> 0010
0110 -> 0110 (gebeurt ook niks)

De rest van de indexen blijven ongemoeid.


Of voor p=0 en q=2:

C++:

3210    3210
============
0000 -> 0000 (hier gebeurt dus niks, maar staat hier voor de volledigheid)
0001 -> 0100
0010 -> 0010 (gebeurt niks)
0011 -> 0110
0100 -> 0001
0101 -> 0101 (gebeurt niks)
0110 -> 0011
0111 -> 0111 (gebeurt niks)

Bij voorbaat hartelijk dank!

Ik zou het doen dmv. het masken van bitposities, zoiets als dit:

C++:

for(i=0; i<ArraLen; i++)
{

    newindex = i;
    if(i&p)
    {
        //  bit positie p set
        newindex |= q;  //  or-en met q
    }
    else
    {
        //  bit positie p NOT set
        newindex &= ~q; //  and-en met de inverse van q
    }

    if(i&q)
    {
        //  bit positie q set
        newindex |= p;  //  or-en met p
    }
    else
    {
        //  bit positie p NOT set
        newindex &= ~p; //  and-en met de inverse van p
    }    
    arrdest[newindex] = arrsrc[i];
}

Je kijkt dus of de gewenste positie geset is of niet en set/unset dan op de andere plaats.

[Verwijderd]

Thanks! Ik wist nog niet van die bitmasks .

[Soultaker]

@Peter_B: zo verwissel je alleen de bits op posities p en q, en niet de bits daartussen! Merk op dat bladiebloe alle bits op posities p tot en met q wil omdraaien.

Ik zou zelf sowieso een look-up table gebruiken, anders moet je gewoon te veel berekenen en ontkom je niet aan het gebruik van een for-lus (wat traag is). Als je array klein is (bijvoorbeeld altijd <256 elementen) dan kun je alle resultaten in de tabel stoppen. Anders zou ik de omgekeerde bytes opslaan en daarmee integers met omgekeerde ranges in elkaar knutselen; zo bijvoorbeeld:

C++:

unsigned char rbytes[256]; // rbytes[x] == byte x in reverse

// reverses a 32-bit integer
inline unsigned reverse(unsigned x)
{
    return rbytes[(x>> 0)&255]<<24 |
           rbytes[(x>> 8)&255]<<16 |
           rbytes[(x>>16)&255]<< 8 |
           rbytes[(x>>24)&255]<< 0;
}

// reverses the bits with positions p to q (inclusive) in x
inline unsigned reverse(unsigned x, unsigned p, unsigned q)
{
    unsigned mask = ((1<<(q+1))-1) ^ ((1<<p)-1);
    return x&~mask | reverse((x&mask)>>p) >> 31-q;
}

#include <iostream>
int main()
{
    // Initialiseer rbytes
    for(unsigned n = 0; n < 256; ++n)
    {
        rbytes[n] = 0;
        for(int b = 0; b < 8; ++b)
            rbytes[n] |= ((n>>(7-b))&1)<<b;
    }

    // Voorbeeld
    for(int n = 0; n < 100; ++n)
    {
        int m = reverse(n, 1,7);

        // Print resultaat
        for(int b = 31; b >= 0; --b)
            std::cout << ((n>>b)&1);
        std::cout << " -> ";
        for(int b = 31; b >= 0; --b)
            std::cout << ((m>>b)&1);
        std::cout << std::endl;
    }
}

[Soultaker]

Overigens kan het nog makkelijker en waarschijnlijk ook efficiënter (dat zou je moeten testen) door gebruik te maken van het feit dat de 'nieuwe' index steeds uit een vaste permutatie van de bits van de 'oude' index bestaat. Bijvoorbeeld: als je bits 1 en 2 wisselt, dan corresponderen de bits 0, 1, 2 en 3 uit de invoer met respectievelijk 0, 2, 1 en 3 uit de uitvoer. Als er bits veranderen in de eerste index, dan hoef je alleen de corresponderende bits in de tweede aan te passen.

Je kunt je zo veel werk besparen omdat het uitzoeken welke bits waar staan één keer gedaan hoeft te worden en je die informatie gedurende de lus kunt gebruiken. Er wordt wel weer een extra lus geïntroduceert, dus of het ook sneller is durf ik niet te zeggen, maar de code is in ieder geval simpeler:

C++:

#include <algorithm>
void transform(int *a, int *b, int size, int p, int q)
{
    // Uitrekenen welke bits in n corresponderen met bits in m
    int bits[33];
    for(int n = 0; n < 32; ++n)
        bits[n] = 1<<n;
    std::reverse(bits + p, bits + q + 1);
    
    for(int n = 0, m = 0; n < size; ++n)
    {
        a[m] = b[n];
        
        // m ophogen
        int b = 0;
        do { m ^= bits[b]; } while((m&bits[b++])==0);
    }
}

[ Voor 11% gewijzigd door Soultaker op 16-06-2006 01:12 ]

[Verwijderd]

Soultaker schreef op vrijdag 16 juni 2006 @ 00:42:
@Peter_B: zo verwissel je alleen de bits op posities p en q, en niet de bits daartussen! Merk op dat bladiebloe alle bits op posities p tot en met q wil omdraaien.

<snip>

Dank je Soultaker! De array die wil gaan hershuffelen is vrij groot, zeg maar van lengte 10⁶. Ik heb wel een paar vragen over je code. Ik zou dus mijn hele array moeten aflopen, steeds m=reverse(n) moeten aanroepen en dan een swap moeten uitvoeren op de elementen array[m] en array[n] ?

En de array rbytes[] bevat dus de reverse van een byte x? En waarmee vul je precies rbytes[] in regel 22?

Edit: ik check nu ff je 2e bericht.

Edit2: Het werkt perfect met je code hierboven! Ik wil het liefst geen extra array gebruiken om de nieuwe vector in op te slaan, dus ik ben nu even aan het kijken over welke elementen de loop moet springen zodat ik niet twee keer elementen swap.

[ Voor 48% gewijzigd door Verwijderd op 16-06-2006 01:34 ]

[Soultaker]

Verwijderd schreef op vrijdag 16 juni 2006 @ 01:13:
Dank je Soultaker! De array die wil gaan hershuffelen is vrij groot, zeg maar van lengte 10⁶. Ik heb wel een paar vragen over je code. Ik zou dus mijn hele array moeten aflopen, steeds m=reverse(n) moeten aanroepen en dan een swap moeten uitvoeren op de elementen array[m] en array[n] ?

Dat was het idee.

En de array rbytes[] bevat dus de reverse van een byte x?

Juist. Dus rbytes[0] = 0, rbytes[1] = 128 (want 00000001₂ = 1 en 10000000₂ = 128), rbytes[2] = 64, rbytes[3] = 192, enzovoorts.

En waarmee vul je precies rbytes[] in regel 22?

Daarmee dus. Ik gebruik een simpel lusje om de bits op de goede plek te zetten, maar omdat dat maar één keer in het hele programma gebeurt kan dat wel. (Je zou ook de lijst kunnen hardcoden; het zijn tenslotte maar 256 waarden.)

Edit2: Het werkt perfect met je code hierboven! Ik wil het liefst geen extra array gebruiken om de nieuwe vector in op te slaan, dus ik ben nu even aan het kijken over welke elementen de loop moet springen zodat ik niet twee keer elementen swap.

Een simpele oplossing is elementen n en m alleen verwisselen als n<m (of natuurlijk m>n); dat zorgt ervoor dat als je een paar n,m twee keer genereert (wat altijd zo is) je 'm één van de twee keer overslaat:

C++:

for(int n = 0, m = 0; n < size; ++n) 
{ 
    if(n < m)
        std::swap(a[n], a[m]);

    // m ophogen 
    int b = 0; 
    do { m ^= bits[b]; } while((m&bits[b++])==0); 
}

Volgens mij kun je nog wel wat beters verzinnen, omdat je nu steeds blokken uit de array verwisselt ter grootte van (1<<p); dat moet nog efficiënter kunnen. Misschien dat ik daar morgen nog op terug kom (als je dat wil) - ik vind het daar nu wat te laat voor.

[ Voor 7% gewijzigd door Soultaker op 16-06-2006 02:02 ]