Szybka Konwolucja 1D bez FFT

Question

Potrzebuję 1K Konwolucji przeciwko 2 dużym tablicom. Używam tego kodu w języku C#, ale uruchomienie trwa krócej.

Wiem, wiem! Zwój FFT jest bardzo szybki. Ale w tym projekcie NIE MOŻNA go używać. Jest to ograniczenie projektu, aby nie używać FFT (proszę nie pytać dlaczego: /).

To jest mój kod w języku C# (przeniesiony z Matlab, nawiasem mówiąc):

var result = new double[input.Length + filter.Length - 1]; 
for (var i = 0; i < input.Length; i++) 
{ 
    for (var j = 0; j < filter.Length; j++) 
    { 
     result[i + j] += input[i] * filter[j]; 
    } 
} 

Więc ktoś zna żadnego szybki algorytm splot widthout FFT?

Answer 1

Splot jest liczbowo samo jako wielomian mnożenia z dodatkowym w pętli wokół etapu. Dlatego wszystkie algorytmy mnożenia wielomianu i dużej liczby całkowitej mogą być użyte do przeprowadzenia splotu.

FFT jest jedynym sposobem uzyskania szybkiego czasu pracy w trybie O (n log (n)). Ale nadal można uzyskać pod kwadratową pracę przy użyciu metod dziel i rządź, takich jak Karatsuba's algorithm.

Algorytm Karatsuby jest dość łatwy do wdrożenia, gdy zrozumiesz, jak to działa. Działa w O (n^1,585) i prawdopodobnie będzie szybszy niż próba superoptymalizowania klasycznego podejścia O (n^2).

Answer 2

Oto dwie możliwości, które mogą powodować drobne przyspieszenia, ale trzeba by przetestować, aby się upewnić.

1. Rozwiąż wewnętrzną pętlę, aby usunąć niektóre testy. Będzie to łatwiejsze, jeśli wiesz, że długość filtra zawsze będzie wielokrotnością, jeśli jakiś numer będzie odwracał kolejność pętli. Czy filter.length przechodzi przez całą tablicę. To powoduje mniej dereferencji w wewnętrznej pętli, ale może mieć gorsze zachowanie podczas buforowania.

Answer 3

Można zmniejszyć liczbę zaindeksowanych dostęp do result, jak również Length nieruchomości:

int inputLength = filter.Length; 
int filterLength = filter.Length; 
var result = new double[inputLength + filterLength - 1]; 
for (int i = resultLength; i >= 0; i--) 
{ 
    double sum = 0; 
    // max(i - input.Length + 1,0) 
    int n1 = i < inputLength ? 0 : i - inputLength + 1; 
    // min(i, filter.Length - 1) 
    int n2 = i < filterLength ? i : filterLength - 1; 
    for (int j = n1; j <= n2; j++) 
    { 
     sum += input[i - j] * filter[j]; 
    } 
    result[i] = sum; 
} 

---

Jeśli dodatkowo podzielone zewnętrznej pętli, można pozbyć się powtarzających niektórych warunkowych . (Zakłada 0 < filterLength ≤ inputLength ≤ resultLength)

int inputLength = filter.Length; 
int filterLength = filter.Length; 
int resultLength = inputLength + filterLength - 1; 

var result = new double[resultLength]; 

for (int i = 0; i < filterLength; i++) 
{ 
    double sum = 0; 
    for (int j = i; j >= 0; j--) 
    { 
     sum += input[i - j] * filter[j]; 
    } 
    result[i] = sum; 
} 
for (int i = filterLength; i < inputLength; i++) 
{ 
    double sum = 0; 
    for (int j = filterLength - 1; j >= 0; j--) 
    { 
     sum += input[i - j] * filter[j]; 
    } 
    result[i] = sum; 
} 
for (int i = inputLength; i < resultLength; i++) 
{ 
    double sum = 0; 
    for (int j = i - inputLength + 1; j < filterLength; j++) 
    { 
     sum += input[i - j] * filter[j]; 
    } 
    result[i] = sum; 
} 

Answer 4

Możesz użyć specjalnego filtra IIR. Następnie przetwarzać że jak:

y(n)= a1*y(n-1)+b1*y(n-2)...+a2*x(n-1)+b2*x(n-2)...... 

myślę, że to szybciej.