Názory k článku Použití instrukcí SSE a AVX pro zrychlení bitových operací

Co se týká autovektorizace, tak na tomto kódu je o dost lepší clang než gcc. Clang zvládne docela dobře autovektorizovat i tu původní úplně neoptimalizovanou verzi a ještě udělá loop unrolling. Možno vyzkoušet zde: https://godbolt.org/z/4qnqsd9E7

Výkonově se naivní ale moc nezlepšila, imho to zabíjí ty vextractps místo aby to rovnou uložilo do cílové destinace.

Ale tu „ručně rozbalenou“ to už zoptimalizuje pěkně a zvládne to uložit přímo vektor. Pro vektory dlouhé 30000 mám na i5-8265U tohle: 106us naivní, 48us ručně rozbalená, 35us ruční SSE, 28us ruční AVX.

Super - pěkné skloubení teorie s praxí!

a proc nezkusit tu paralelizaci?
kazde vlakno/proces by si samo nezavisle zpracovalo nejaky usek dat, treba usek 250 ms a jen vysypaly nekam vysledek.

Protože teď mám jednoduché „vytáhnu prvek z fronty, zpracuju, vyklopím do ZeroMQ socketu“. Když bude workerů víc, doběhnout v náhodném pořadí a někdo to bude muset serializovat (aby se to do toho socketu zapisovalo ve správném pořadí). Samozřejmě to není nic neřešitelného, ale vypadalo to jako víc práce než vyřešit podivně pomalé čtení int16.

A taky jak píšu, tohle není jedná věc co běží, z těch výsledných dat se ještě musí vyrobit netriviální produkty (to nevím jak snadno rychle popsat a žádný intro článek jsem nenašel; ve zkratce získané odpovědi od cíle z jednoho azimutu se dají vedle sebe a pak se nad tím dělají různé operace do určí Dopplerův posuv, odstraní věci s nulovou Dopplerovskou frekvencí (ground clutter), a dělají nějaké šílené heuristiky pro odhadnutí věcí jako KDP kterým vůbec nerozumím ale naštěstí máme v týmu lidi co to umí :) a to generování produktů by mělo na tom počítači zvládnout běžet dvakrát nebo ideálně třikrát současně (produkční/vý­vojová).

Ocenujem clanok na takuto low-level tematiku, ale popravde, robi viac skody nez osohu.

Ano, je skvele, ze dnesne procesory maju plno sikovnych instrukcii. Ale co tak si spomenut, ze prvy krok sa vola optimalizacia algoritmu?

Dal som tomu 15 minut svojho casu:

for (int i = 0; i < SAMPLES*4; i += 4) {
uint64_t* raw = (uint64_t*)&inbuf[i + 0];
uint64_t mask = (*raw >> 1) & 0x1000100010001000;
*raw |= mask;

const int i2 = i >> 1;

outbuf1[i2] = (float)inbuf[i + 0];
outbuf1[i2 + 1] = (float)inbuf[i + 1];

outbuf2[i2] = (float)inbuf[i + 2];
outbuf2[i2 + 1] = (float)inbuf[i + 3];
}

Ako vidno, absolutne som sa nesnazil to optimalizovat na registre ci davat kompilatoru nejake hinty a vysledok je takmer 2x rychlejsi kod.

Nemozem uverit, ze autor radsej vyzaduje najnovsiu generaciu procesorov, ktore mu usporiadaju 4 16-bitove slova do jedneho miesto toho, aby si to napisal sam.

Nehovoriac o tom, ze ak by si autor ten bit #12 (to cislovanie 1az 13 v clanku je dost... nestandardne) dal na trochu vhodnejsie miesto (cize predpripravil data vo vhodnejsom formate), tak to otvara dvere dalsim optimalizaciam.

Bez jedinej SSE instrukcie by sa ta slucka dala urychlit 3x - 4x.

Ten tvůj program porušuje strict aliasing rule, nemůžeš do paměti ukládat přes uint64_t pointer a pak to číst jako int16_t pointer, to je undefined behavior, není to korektní kód.

To, že je to undefined behaviour v daném jazyce, neznamená, že je to undefined na cílové/ých architekturách.

Tohle s architekturou nesouvisí. C++ překladač nemusí invalidovat data v registrech (načítat je znovu), pokud se do paměti ukládá přes pointer jiného typu. Je to optimalizace. Uvedený program tohle porušuje, ukládá do paměti int64 hodnotu a potom ze stejného místa čte int16 hodnotu. Překladač ale z paměti int16 hodnotu nemusí načíst, protože strict aliasing rule říká, že se nesmí aliasovat pointery různého typu.

Musí, protože dřív to (na dané architektuře a kompilátoru) fungovalo a staré kódy s tím počítají. Warning/error na stric aliasing si nastavuješ, abys odhalil dopředu problémy, kdybys chtěl portovat na jinou architekturu (např. na MIPS v grafické stanici SGI nefungoval a chvíli jsem se trápil, proč ten kód dělá nesmysl).

Je to tak. Preto ta poznamka, ze som tomu nevenoval viac usilia, bolo to len na demonstraciu, ze dokonca aj pri manipulacii s pamatou je ten algoritmus rychlejsi len s minimalnymi upravami.

Pointa toho prikladu je, ze keby sa autor trochu zamyslel miesto prehladavania intelackych tabuliek s instrukciami, tak moze mat zadanie napisane "vektorovo" aj v plain C-cku na obycajnej 64-bitovej architekture zpred 10 rokov.

23. 11. 2022, 13:18 editováno autorem komentáře

Ten tvůj program má dvě chyby, dělá něco jiného a je v něm undefined behavior.

Děláš tam OR do metadat bitu, ale ten bit předtím nevynuluješ, takže to nefunguje správně.

Pointa je, že nemá smysl srovnávat výkon korektní implementace a tvé implementace, která dělá něco jiného a ještě s undefined behavior.

> Děláš tam OR do metadat bitu, ale ten bit předtím nevynuluješ, takže to nefunguje správně.

Safra, to jsem pěkně rozbil testy. Nechtěl jsem do gitu commitovat velký soubor z /dev/urandom na kterém jsem to testoval, ani se mi nechtělo shánět deterministický RNG co bych tam dal (vím že v Pythonu můžu nastavit seed, ale nevěřím tomu, že to vygeneruje stejné věci napříč verzemi), tak mě napadl ten hack "uděláme nějaký 'seq' a to použijeme jako data". Jenže co čert nechtěl, zrovna v ASCII hodnotách výstupu není ten klíčový bit nastaven, takže testujeme jenom že se to správně přioruje, ale ne že se vymaže pokud nastaven je.

Nápady na deterministický RNG pro generování testovacích binárních dat? Změnil jsem to na samodomo věc co opakovaně počítá MD5 a vypisuje výsledek.

Dopsal jsem mu tam ten jeden AND s invertovanou maskou, takže už to počítá totéž a zpomalilo se to zanedbatelně.

Ten UB se dá celkem jednoduše obejít pomoci memcpy().

Třeba:

static inline uint64_t load64_unaligned(const void* src) {
  uint64_t x;
  memcpy(&x, src, sizeof(uint64_t));
  return x;
}

Díky, přidal jsem to do testovacího programu.

i7-9700TE, GCC 10 (z Debianu stable): numpy 163us, naivní 294us, ručně unroll 92us*, tvoje 73us, moje SSE 31us, moje AVX 24us.

* tohle je zajímavé, protože ve standalone testovacím programu (co je na githubu) ruční unroll pomohl hodně, ale v celém tom zpracování vůbec. IMHO se nějak blbě zpropagovaly restrict nebo další kvalifikátory a ve standalone se to podařilo optimalizovat. Ale výsledný kód se mi asi teď úplně zkoumat nechce.

Strict aliasing se dá řešit pomocí unions těch datových typů, ne?

A ještě k tomuhle:

> Nemozem uverit, ze autor radsej vyzaduje najnovsiu generaciu procesorov

Ve skutečnosti použitá SSE umí snad všechny procesory vydané po roce ~2008 a na ničem jiném fakt nepoběžíme, pokud teda nepřejdeme na nějaký M1 nebo tak něco :)

> Nehovoriac o tom, ze ak by si autor ten bit #12 (to cislovanie 1az 13 v clanku je dost... nestandardne) dal na trochu vhodnejsie miesto (cize predpripravil data vo vhodnejsom formate), tak to otvara dvere dalsim optimalizaciam.

To by vyžadovalo trochu větší hackování FPGA kódu v bladeRF, což je pro mě komplikované, protože FPGA nerozumím. A nevím moc jak by to pomohlo -- řekněme že bych to dal do nejvyššího bitu, furt musím udělat v podstatě totéž, vymaskovat a nahradit.

protože FPGA nerozumím
To by nebolo zle sa doucit, pretoze na FPGA urobite aj tu fourierovu transformaciu daleko efektivnejsie ako v strojaku...

Ono se to FFT zvládá, výhodu ve zpracování přímo ve FPGA bych viděl spíš v tom, že by se obešel bottleneck propojení s počítačem -- aktuálně je to 20 MS/s * 4 bajty na vzorek * MIMO2x2 = 160 MB/s pro každý směr (RX a TX) což už se blíží tomu co dá USB 3.0 a přitom ta data obsahují podstatně méně informace, většina z toho jsou „nuly“, resp. oversampling. Navíc by to mohlo přímo posílat balíčky otagované azimutem, takže bych nemusel řešit to lepení metadat do nejvyšších bitů o kterém je článek.

Velikou nevýhodou pak ale je, že tohle bylo nabastlené za chvilku, zatímco to FPGA by klidně mohlo být několik měsíců práce - a já samozřejmě nemůžu řešit jenom tohle, protože tak nějak stavím celý radar (teď už samozřejmě ve víc lidech, ale stejně). Právě to, že jsme na rozdíl od konkurenčních firem použili jako jediní v podstatě nemodifikované hotové SDR nám nejspíš umožnilo vyrobit radar s omezenými prostředky v rekordním čase. Dále, hrábnout do mého C/Python kódu může kdokoli z týmu nebo i zákazníků (když si to třeba koupí research tým na univerzitě a chce si něco vyzkoušet), zatímco FPGA rozumí málokdo a dělat tam nějaké úpravy je peklo. Poslední myšlenka pak je, že my nevíme, jestli třeba nenajdeme nějaký lepší způsob filtrace/zpracování signálu - a opět, v kódu na počítači se to snadno vyzkouší a změní.

Ještě doplním jeden problém: když se zjistí že bladeRF / to FPGA co se v něm používá nejde koupit (to se teď stává doslova se vším a není výjimkou číst zprávy o tom že to zastavilo třeba celou ohromnou automobilku; od ostatních výrobců radarů slyšíme přesně totéž) a musím použít jiné rádio (což bude znamenat jinou řadu a třeba i jiného výrobce FPGA (Intel / Xilinx)), tak mě čeká velmi komplikovaná portace. Software oproti tomu běží všude.

Totéž pokud se zjistí že bladeRF nevyhovuje a je potřeba něco jiného.

Jj, to dava rozum, akurat mi prislo logicke na to pouzit to fpga. A ako som letmo pogooglil, tak zdaleka nie som sam, takze mesiace znovuvyliezania kolesa by nehrozily...

Parádní praktická ukázka. Moc díky.

skvely clanek. Diky za nej. Alespon je videt, jak konkretne dokazou tyto instrukce zoptimalizovat program a kolik casu se tim usetri.
K tomu prevodu int16 -> float, ja bych na to vyuzil look up table. Mas tady pouze 12 bitu plus znamenkovy bit. To je 8k float hodnot. Kdyz k tomu pridame i ten bit, co je nevyznamny (a v tabulce ho budeme ignorovat), budeme mit 16k float hodnot. To neni tak moc. Jestli mas chut zkusit, jak by to vypadalo casove, pripadne jestli by se dal i tento program zoptimalizovat pomoci SSE, napis potom, jak ten pokus dopadl.
Nezkousel jsem to prelozit, doufam, ze tam nemam nejaky preklep

float lut[0x4000];
void prepare_lut()
{
  for (int i=0; i<0x4000; i++)
  {
    int16_t nr=(int16_t)(((i & 0x2000)?0xf000:0)|(i & 0x0fff));
    lut[i]=(float)nr;
  }
}

void convert()
{
  float *op1=outbuf1;
  float *op2=outbuf2;
  int16_t *ip=inbuf;
  while (ip<inbuf+4*SAMPLES)
  {
    *op1++=lut[*ip++ & 0x3fff];
    *op1++=lut[*ip++ & 0x3fff];
    *op2++=lut[*ip++ & 0x3fff];
    *op2++=lut[*ip++ & 0x3fff];
  }
}

Lookup table má dva problémy, způsobuje cache thrashing a nejde moc dobře vektorizovat.

Cache thrashing bývá u takových algoritmů s požadavky na vysokou propustnost docela problém, protože je nutné, aby bylo pokud možno všechno v L1 cache, jinak hodně rostou latence, když se čeká na paměť.

Pro konverze int na float existují SIMD instrukce s relativně nízkou latencí, např. _mm_cvtepi32_ps, což je obvykle lepší volba, než lookup table.

docela to vymlati L1 cache, takze by to fakt chtelo otestovat v realnem benchmarku. Toto se delavalo na osmibitech, dneska vetsinou takove LUT dopadaji hur (divna doba :)

Jeste bych to upravil, maskoval bych to 2ffff, tim padem se bude uplne ignorovat ten extra bit a nebude se pouzivat polovina tabulky. Tim by se mel snizit narok na cache.

Stačí maskovat 0x0FFF, celé hodnota je určena nejnižšími 12 bity (ADC vrací 4096 možností a kódování do „int12“ to nerozbije). Takže tabulka má jen 16KiB a běží to skoro tak rychle jako SSE implementace. Dal jsem to také do toho gitového repozitáře.

Hustý, implementoval jsem to (s tou optimalizací že tabulka stačí 4096 velká, protože metadata při lookupu vyandujeme) a je to skoro stejně rychlé jako SSE! (35 vs. 31us)

Tak se mi líbí jak jednoduše vyjadřuješ ty adresy pro destinaci, mě tohle nikdy nenapadne :)

Tabulka má jenom 16 KiB (4096 * sizeof(float)) takže s cache taky v pohodě.

Není tam chyba?


int idx = i/2 + imod - 3;
//vs:
outbuf2[i2 + 1] = (float)raw3;


IMHO v tom prvním má být imod-2 . Nezkoušel jsem, jestli to nějak zlepší výsledek autovektorizace.

Ohledně poznámky výše a SSE vs aarch64 - na Apple M1 netřeba, stačí zkusit Graviton3 v AWS - podporuje Neon i SVE. S těmi 2+2 konverzemi si nejsem jistý, jestli to zvládne efektivně. Ještě by mělo smysl zkusit AVX-512, které je taky už nějakou dobu standard, aspoň u Intel.

PS: Pěkný článek a hezké ne zcela běžné využití SIMD.

Jen doplním, že M1 má "plnotučná" jádra, která si poradí i s neoptimálním kódem (podobně jako x86 jádra od Intel a AMD). Zatímco vše ostatní jsou jádra Cortex, která tak chytrá (ale taky velká) nejsou. Příkladem jsou optimalizace, které nemají na M1 vliv (neoptimální kód tam běží stejně rychle jako optimální). Samozřejmě pro Amazon a cloud obecně je lepší více menších hloupějších jader a optimalizovat software na míru (variable vs fixed costs).

Je pravda, že u násobení matic mi přišly třeba Graviton na Neon pomalejší než M1, se SVE ale byly rychlejší.

Nicméně, zkusil jsem pro zajímavost M1 a Graviton3:

#include <stdint.h>
#include <stdlib.h>

#ifdef __aarch64__
#include <arm_neon.h>
#endif


void radar_ref(float *outbuf1, float *outbuf2, int16_t *inbuf, int SAMPLES)
{
    for (size_t i = 0; i<4*SAMPLES; i++) {

        // fix sign-extend
        int16_t raw = inbuf[i];
        raw = (raw & 0xEFFF) | ((raw & 0xE000)>>1);

        // make float from int16
        float raw_f = (float)raw;

        // decide where to put it - horizontal or vertical buffer
        int imod = i % 4;
        if(imod < 2) {
            int idx = i/2 + imod;
            outbuf1[idx] = raw_f;
        } else {
            int idx = i/2 + imod - 3;
            outbuf2[idx] = raw_f;
        }
    }
}

void radar_Neon(float *outbuf1, float *outbuf2, int16_t *inbuf, int SAMPLES)
{
    int16x4_t four16_EFFF = vdup_n_s16(0xEFFF);
    int16x4_t four16_E000 = vdup_n_s16(0xE000);

    for (size_t b = 0; b < SAMPLES; b += 1) {

        // fix sign-extend
        int16x4_t raw4 = ((int16x4_t *)inbuf)[b];
        int32x4_t adjusted = vmovl_s16(vorr_s16(vand_s16(raw4, four16_EFFF), vshr_n_s16(vand_s16(raw4, four16_E000), 1)));

        // make float from int16
        float32x4_t flt = vcvtq_f32_s32(adjusted);

        vst1_f32(outbuf1+b*2, vget_low_f32(flt));
        vst1_f32(outbuf2+b*2, vget_high_f32(flt));
    }
}

// main.c :
// (too smart compiler may completely remove code doing nothing)
#include <stdint.h>
#include <stdlib.h>

extern void radar_ref(float *outbuf1, float *outbuf2, int16_t *inbuf, int SAMPLES);
extern void radar_Neon(float *outbuf1, float *outbuf2, int16_t *inbuf, int SAMPLES);

#define SAMPLES 30000
#define ITERATIONS 1000000


int main(void)
{
    float outbuf1[SAMPLES*2], outbuf2[SAMPLES*2];
    int16_t inbuf[SAMPLES*4];

    for (int i = 0; i < ITERATIONS; ++i) {
        radar_ref(outbuf1, outbuf2, inbuf, SAMPLES);
    }
    return 0;
}

Nejspíš by šlo vylepšit zpracováním osmi prvků najednou, podobně jako v původním AVX. Případně šestnácti prvků a využít load and store pair.
Neznám tak dobře shuffle instrukce na Neon, abych to dal z hlavy.

Apple M1 Pro (clang 14.0.0) :
Letmým pohled, clang byl schopen autovektorizovat aspoň operace po dvou prvcích.
ref: 44 us
Neon: 20 us

Graviton3 c7g.medium (gcc 12.1.0) :
Přímý kód, žádná autovektorizace.
ref: 183 us
Neon: 41 us

Pěkný článek. Pro daný účel bych zkusil výpočet na GPU ;-). Rád programuji v C++, x64 asm a AVX2.

Aktuálně se pokouším naprogramovat prolomení WiFi hesla o délce 8 znaků (malá/velká písmena a číslice) hrubou silou, což je 62^8 tj. cca 200 bilionů kombinací. Zkušebně v C++ CUDA programu (Visual Studio 2022) všechna hesla vygeneruji za cca 40 ms (notebook má NVIDIA GeForce GTX 960M). Umím už v C++ spočítat SHA-1 a chystám se doplnit HMAC a PBKDF2, abych našel správné heslo v rozumném čase.

#include "cuda_runtime.h"
#include "device_launch_parameters.h"

#include <stdio.h>
#include <iostream>
#include <chrono>

using namespace std;
using namespace std::chrono;

# define blocks 4
# define threads 992
# define characters 8

cudaError_t cudaStatus;

// "abcdefghijklmnopqrstuvwxyzABCDEFGHIJKLMNOPQRSTUVWXYZ0123456789"
__constant__ char1 charset[] = { 0x61, 0x62, 0x63, 0x64, 0x65, 0x66, 0x67, 0x68, 0x69, 0x6a, 0x6b, 0x6c, 0x6d, 0x6e, 0x6f, 0x70, 0x71, 0x72, 0x73, 0x74, 0x75, 0x76, 0x77, 0x78, 0x79, 0x7a, 0x41, 0x42, 0x43, 0x44, 0x45, 0x46, 0x47, 0x48, 0x49, 0x4a, 0x4b, 0x4c, 0x4d, 0x4e, 0x4f, 0x50, 0x51, 0x52, 0x53, 0x54, 0x55, 0x56, 0x57, 0x58, 0x59, 0x5a, 0x30, 0x31, 0x32, 0x33, 0x34, 0x35, 0x36, 0x37, 0x38, 0x39 };

__global__ void blackcat(void) {

        char1 password[characters];

        uint8_t counters[characters];
        uint64_t n = (pow(62, characters) / threads);   // Number of search cycles per thread

        // Nastavení počátečních hodnot hesla pro každé vlákno, odkud se mají začít generovat
        for (int i = characters - 1; i >= 0; i--) {
                counters[i] = (n * threadIdx.x / (uint64_t)pow(62, characters - 1 - i) % 62);
        }

        while (n > 0) {

                bool flag = false;
                for (int i = characters - 1; i >= 0; i--) {
                        password[i] = charset[counters[i]];
                        if (i == characters - 1) {
                                counters[i]++;
                                if (counters[i] > 61) {
                                        counters[i] = (uint8_t)0;
                                        flag = true;
                                }
                        }
                        else {
                                if (flag) {
                                        counters[i]++;
                                        if (counters[i] > 61) {
                                                counters[i] = (uint8_t)0;
                                        }
                                        else {
                                                flag = false;
                                        }
                                }
                        }
                }
                // Po odkomentování vypíše poslední 3 generované hesla
                //if (threadIdx.x == threads - 1 && blockIdx.x == blocks - 1 && n < 4) {
                //      printf("Thread[%d]",threadIdx.x);
                //      for (int i = 0; i < characters; i++) {
                //              printf(" %c", password[i]);
                //      }
                //      printf("\n");
                //}

                /* Test zda jsme našli password,
                pokud ano vypíšeme password, ukončíme všechna vlákna a předčasně se vrátíme z funkce,
                možná bude dobré občas vypsat čas běhu, abychom věděli, že program stále běží */
                n--;
        }
}

int main() {

        auto start = high_resolution_clock::now();

        cudaSetDevice(0);
        cudaStatus = cudaGetLastError();
                if (cudaStatus != cudaSuccess) {
                fprintf(stderr, "cudaSetDevice failed!  Do you have a CUDA-capable GPU installed?");
         }
        blackcat << <blocks, threads >> > ();
        cudaStatus = cudaGetLastError();
                if (cudaStatus != cudaSuccess) {
                fprintf(stderr, "Kernel launch failed: %s\n", cudaGetErrorString(cudaStatus));
        }
        cudaDeviceSynchronize();
        cudaStatus = cudaGetLastError();
                if (cudaStatus != cudaSuccess) {
                fprintf(stderr, "cudaDeviceSynchronize returned error code %d after launching addKernel!\n", cudaStatus);
         }

        auto stop = high_resolution_clock::now();
        auto duration = duration_cast<microseconds>(stop - start);
        printf("\nTime  = %llx (HEX)\n", duration.count());

        return 0;
}

Pre zaujímavosť som skúsil Visual Studio 2022 na Windows. Nejaké sse/avx (void*) som musel fixnúť na korektné casty ako napr. _mm_stream_lo­ad_si128((con­st __m128i*)(inbuf + i)).
Gcc tam nemal warningy?
Vysledky na mobilnom AMD Ryzen 7 PRO 2700U 2200Mhz:

naive:188
unroll:102
lut:65
unroll2: 55
sse:42
avx:37

Unroll2 je moje zjednodušenie unroll ():

int oi = 0;
for (int i = 0; i < 4 * SAMPLES; i += 4) {
int16_t raw0 = inbuf[i++];
raw0 = (raw0 & 0xEFFF) | ((raw0 & 0xE000) >> 1);

int16_t raw1 = inbuf[i++];
raw1 = (raw1 & 0xEFFF) | ((raw1 & 0xE000) >> 1);

int16_t raw2 = inbuf[i++];
raw2 = (raw2 & 0xEFFF) | ((raw2 & 0xE000) >> 1);

int16_t raw3 = inbuf[i++];
raw3 = (raw3 & 0xEFFF) | ((raw3 & 0xE000) >> 1);

outbuf1[oi] = (float)raw0;
outbuf2[oi] = (float)raw2;
outbuf1[++oi] = (float)raw1;
outbuf2[oi] = (float)raw3;
}

Kompiler tam pouzil nejake sse, je tam 4x instrukcia cvtdq2ps xmm0,xmm0

Základem pro rychlé zpracování dat je použít GPU.

Do vyhledávače google zadat "GPU CUDA Radar MIMO 2×2" a výsledky hledání omezit na 1 rok.

Viz. například
MIMO Radar Parallel Simulation System Based on CPU/GPU Architecture
Design of high-speed software defined radar with GPU accelerator

Pro lepší vektorizaci ručně rozvinutého kódu by mělo stačit zbavení se závislosti i2 na i, třeba nezávislým výpočtem ve for:

#pragma GCC ivdep
for(int i = 0, i2 = 0; i<4*SAMPLES; i+=4, i2+=2) {
  int16_t raw0 = inbuf[i+0];
  int16_t raw1 = inbuf[i+1];
  int16_t raw2 = inbuf[i+2];
  int16_t raw3 = inbuf[i+3];

  raw0 = (raw0 & 0xEFFF) | ((raw0 & 0xE000)>>1);
  raw1 = (raw1 & 0xEFFF) | ((raw1 & 0xE000)>>1);
  raw2 = (raw2 & 0xEFFF) | ((raw2 & 0xE000)>>1);
  raw3 = (raw3 & 0xEFFF) | ((raw3 & 0xE000)>>1);

  outbuf1[i2    ] = (float)raw0;
  outbuf1[i2 + 1] = (float)raw1;

  outbuf2[i2    ] = (float)raw2;
  outbuf2[i2 + 1] = (float)raw3;
}

Případně zpracovává-li se sekvenčně velké množství dat, která se nevejdou do L2 cache, tak by měl rychlejšímu běhu pomoci gcc přepínač -fprefetch-loop-arrays.