Chciałbym zrobić implementację filtracji FIR do plików WAV, które mogą zawierać dźwięk, ale też mogą symulować sygnał szerokopasmowy (np. sygnał próbkowany z częstotliwością 22050000Hz może być zapisany w pliku z próbkowaniem 44100Hz, przy czym w odtwarzaczach audio taki plik będzie odtwarzany w 500-krotnie zwolnionym tempie).
Zauważyłem, że filtracja trwa długo i próbuję ją zoptymalizować czasowo. Algorytm jest oparty o bufor pierścieniowy, dla filtrów symetrycznych (na przykład https://en.wikipedia.org/wiki/Sinc_filter ) w rzeczywistości czynność wykonuje się na połowie okna, przy czym każdy współczynnik za wyjątkiem zerowego jest wykorzystywany dwa razy (zerowy jest używany jeden raz).
Czy coś jeszcze da się zoptymalizować czasowo?
Plik testowy maja ok. 500MB, jest to jednokanałowy WAV zawierający ok. 300 milionów próbek 16-bitowych. Jak się zakomentuje kod od komentarza START, do komentarza STOP, to kod trwa ok. 10 sekund. Jest to do zaakceptowania, więc chciałbym się skupić na optymalizacji samej filtracji. Dla tego przykładu trwa ponad 10 minut dla filtru zawierającego 500 współczynników, co daje rzeczywisty filtr mający okno o długości 1001 współczynników.
Można interpretować plik testowy jako sygnał o próbkowaniu 22 mln próbek na sekundę, jak się "rozciągnie" plik źródłowy 500 razy. Po takim rozciągnięciu trzeba zrobić filtr dolno przepustowy, bardzo wąskie filtry przepustowe wymagają dużej liczby współczynników w celu dokładnej filtracji.
Niezaprzeczalnym dowodem na to, że da się filtrować takie sygnały komputerowo w czasie rzeczywistym jest oprogramowanie do odbiorników radiowych SDR, które oprócz samej filtracji wykonują wiele innych czynności z sygnałem o częstotliwości próbkowania liczonej w megahercach i to wszystko w czasie rzeczywistym. W takim razie musi się dać to jeszcze bardziej zoptymalizować. Nie chodzi o optymalizację samego filtra (zmniejszenie liczby współczynników kosztem dokładności), ale o optymalizację filtrowania zadanym filtrem FIR. Interesuje mnie optymalizacja kodu pomiędzy komentarzami START i STOP.
Poniżej kod procedury dla filtrów symetrycznych, okno filtru jest w tablicy FilterWindow typu llong[], w której na pozycji 0 jest współczynnik środkowy, a na kolejnych pozycjach są współczynniki kolejno dalsze od środka. Innymi słowy, ta tablica zawiera połowę okna filtru wraz z elementem środkowym.
void AppCore::ProcessPassFilter(string FileI, string FileO)
{
// Rzeczywista długość filtru dla danej liczby wspolczynnikow
// W tym przypadku kazdy filtr ma okno symetryczne o nieparzystej liczbie probek
// FilterSize to liczba wspolczynnikow okna filtru nie liczac wspolczynnika zerowego
int FilterWindowWorkSize = (FilterSize * 2) + 1;
// Zerowanie bufora roboczego, ktory jest typu short[] i bedzie zawierac
// dokladnie tyle ostatnich probek z pliku zrodlowego, jak dlugi jest filtr
for (int I = 0; I < FilterWindowWorkSize; I++)
{
FilterWindowWork[I] = 0;
}
// Przygotowanie plikow, DataI to plik zrodlowy, DataO to plik docelowy
EdenClass::WaveFile DataI;
EdenClass::WaveFile DataO;
DataI.FileOpen(FileI, true);
DataO.FileNew(FileO, DataI.SampleRate, 1, 0);
// Dlugosc pliku zrodlowego
llong FileSize = DataI.Length;
// Bufor do odczytu i zapisu danych, ChunkSize jest typu llong
// i ma wartosc 1024, ustawiana przed wywolanie procedury
short *Buf = new short[ChunkSize];
// Miejsce odczytu i zapisu w pliku
llong FileOffset = 0;
// Wskaznik miejsca w buforze roboczym wskazujacy miejsce zapisu odczytanych probek
int WindowPointer = 0;
// Wskazniki miejsc w buforze roboczym wskazujace miejsca odczytu probek przy filtracji
int WindowPointer1 = 0;
int WindowPointer2 = 0;
// Petla glowna
while (ChunkSize > 0)
{
// Odczyt porcji danych z pliku zrodlowego
DataI.ReadM(Buf, FileOffset, ChunkSize);
// START
// Przetwarzanie odczytanej porcji danych, jedna iteracja to jedna probka
for (int II = 0; II < ChunkSize; II++)
{
Zapisanie probki do bufora roboczego
FilterWindowWork[WindowPointer] = Buf[II];
// Ustawianie wskaznikow odczytu probek do filtru na miejsce srodkowe miedzy pierwsza a ostatnia odczytana probka
WindowPointer1 = WindowPointer - FilterSize - 0;
if (WindowPointer1 < 0) { WindowPointer1 += FilterWindowWorkSize; }
WindowPointer2 = WindowPointer1;
// Obliczenie wartosci probki na podstawie wspolczynnika zerowego (znajdujacego sie w osi symetrii okna filtru)
llong Temp = FilterWindowWork[WindowPointer1];
llong XXX = (Temp * FilterWindow[0]);
// Przesuniecie wskaznikow odczytu o jedna probke, jeden do przodu, drugi do tylu,
// ale z zachowanie pierscieniowego charakteru bufora
WindowPointer1++;
if (WindowPointer1 == FilterWindowWorkSize) { WindowPointer1 = 0; }
if (WindowPointer2 == 0) { WindowPointer2 = FilterWindowWorkSize; }
WindowPointer2--;
// Uwzglednianie pozostalych wspolczynnikow filtru
for (int III = 1; III < FilterSize; III++)
{
// Obliczanie skladowych wartosci probki dla pozostalych wspolczynnikow
// i dodawanie obliczonych wartosci do wartosci na podstawie wspolczynnika zerowego
Temp = FilterWindowWork[WindowPointer1];
XXX = XXX + (Temp * FilterWindow[III]);
Temp = FilterWindowWork[WindowPointer2];
XXX = XXX + (Temp * FilterWindow[III]);
// Przesuniecie wskaznikow odczytu o jedna probke, jeden do przodu, drugi do tylu,
// ale z zachowanie pierscieniowego charakteru bufora
WindowPointer1++;
if (WindowPointer1 == FilterWindowWorkSize) { WindowPointer1 = 0; }
if (WindowPointer2 == 0) { WindowPointer2 = FilterWindowWorkSize; }
WindowPointer2--;
}
// WindowAmplitude to wartosc stala 65536
XXX = XXX / (WindowAmplitude);
// Ograniczanie wartosci obliczonej probki w celu unikniecia bledu z powodu zaokraglen
if (XXX > 32767) { XXX = 32767; }
if (XXX < -32768) { XXX = -32768; }
// Zapisanie obliczonej probki do bufora do zapisu danych
Buf[II] = XXX;
// Przesuniecie wskaznika zapisu do bufora roboczego o jedna probke do przodu
// z zachowaniem pierscieniowego charakteru bufora
WindowPointer++;
if (WindowPointer >= FilterWindowWorkSize) { WindowPointer = 0; }
}
// STOP
// Zapis porcji danych do pliku docelowego
DataO.WriteM(Buf, FileOffset, ChunkSize);
// Przesuwanie miejsca w pliku o dlugosc ostatniej porcji danych
FileOffset += ChunkSize;
// Ostatnia porcja danych moze byc krotsza w przypadku, gdy dlugosc pliku
// nie jest podzielna przez dlugosc porcji, a po przetworzeniu ostatniej porcji
// jej nowa wielkosc wynosi 0 i to powoduje zakonczenie petli
if ((FileOffset + ChunkSize) > FileSize)
{
ChunkSize = FileSize - FileOffset;
}
}
// Niszczenie bufora do odczytu i zapisu danych
delete[] Buf;
// Zamkniecie plikow
DataI.FileClose();
DataO.FileClose();
}
Poniżej kod procedury dla dowolnego filtru FIR, w komentarzach są opisane tylko różnice w stosunku do powyższej procedury.
W tym przypadku, tablica FilterWindow zawiera całe okno filtru.
void AppCore::ProcessUserFilter(string FileI, string FileO)
{
// W tym przypadku FilterSize to rzeczywista dlugosc okna filtru
for (int I = 0; I < FilterSize; I++)
{
FilterWindowWork[I] = 0;
}
EdenClass::WaveFile DataI;
EdenClass::WaveFile DataO;
DataI.FileOpen(FileI, true);
DataO.FileNew(FileO, DataI.SampleRate, 1, 0);
llong FileSize = DataI.Length;
short *Buf = new short[ChunkSize];
llong FileOffset = 0;
int WindowPointer = 0;
// Jest tylko jeden wskaznik odczytu z bufora roboczego
int WindowPointer1 = 0;
ProgressAll = FileSize / ChunkSize;
while (ChunkSize > 0)
{
DataI.ReadM(Buf, FileOffset, ChunkSize);
// START
for (int II = 0; II < ChunkSize; II++)
{
// Zapisanie probki do bufora roboczego
FilterWindowWork[WindowPointer] = Buf[II];
// Ustawianie wskaznika odczytu probek do filtru na miejsce za ostatniej odczytana probka,
// bedzie wskazywac wartosc pierwszej (najstarszej) odczytanej probki
WindowPointer1 = WindowPointer - FilterSize - 0;
if (WindowPointer1 < 0) { WindowPointer1 += FilterSize; }
llong Temp = 0;
llong XXX = 0;
if (WindowPointer1 == FilterSize) { WindowPointer1 = 0; }
for (int III = 0; III < FilterSize; III++)
{
// Obliczanie skladowych wartosci probki dla wszystkich wspolczynnikow
// i dodawanie obliczonych wartosci do wartosci na podstawie wspolczynnika zerowego
Temp = FilterWindowWork[WindowPointer1];
XXX = XXX + (Temp * FilterWindow[III]);
WindowPointer1++;
if (WindowPointer1 == FilterSize) { WindowPointer1 = 0; }
}
XXX = XXX / (WindowAmplitude);
if (XXX > 32767) { XXX = 32767; }
if (XXX < -32768) { XXX = -32768; }
Buf[II] = XXX;
WindowPointer++;
if (WindowPointer >= FilterSize) { WindowPointer = 0; }
}
// STOP
DataO.WriteM(Buf, FileOffset, ChunkSize);
FileOffset += ChunkSize;
if ((FileOffset + ChunkSize) > FileSize)
{
ChunkSize = FileSize - FileOffset;
}
ProgressCurrent += ProgressStep;
}
delete[] Buf;
DataI.FileClose();
DataO.FileClose();
}