1D Convolución rápida sin FFT

Question

Necesito una Convolución 1D contra 2 grandes matrices. Estoy usando este código en C#, pero lleva mucho tiempo ejecutarlo.

Lo sé, lo sé! Las convoluciones de FFT son muy rápidas. Pero en este proyecto NO PUEDO usarlo. Es una limitación del proyecto no usar FFT (no pregunte por qué: /).

Este es mi código en C# (portado de Matlab, por cierto):

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]; 
    } 
} 

tanto, cualquier persona sabe cualquier algoritmo de convolución rápida widthout FFT?

Answer 1

convolución es numéricamente igual a una multiplicación de polinomios con un extra envolvente paso. Por lo tanto, todos los algoritmos de multiplicación de números enteros y polinomiales se pueden usar para realizar convolución.

FFT es la única manera de conseguir la rápida O (n log (n)) en tiempo de ejecución. Pero aún puede obtener tiempo de ejecución sub-cuadrático usando los enfoques de dividir y vencer como Karatsuba's algorithm.

algoritmo de Karatsuba es bastante fácil de poner en práctica una vez que entienda cómo funciona. Se ejecuta en O (n^1.585), y probablemente será más rápido que tratar de super-optimizar el enfoque O (n^2) clásico.

Answer 2

Aquí hay dos posibilidades que pueden dar aceleraciones menores, pero se necesitaría para poner a prueba para estar seguro.

1. Desenrolle el lazo interno para eliminar algunas pruebas. Esto será más fácil si conoce la longitud del filtro será siempre un múltiplo si alguna N.
2. invertir el orden de los bucles. Do filter.length pasa sobre toda la matriz. Esto hace menos desreferenciación en el ciclo interno, pero puede tener un peor comportamiento de almacenamiento en caché.

Answer 3

Se podría reducir el número de indexado accesos a result, así como los Length propiedades:

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; 
} 

---

Si divide aún más el lazo externo, usted puede deshacerse de algunos condicionales que se repiten . (Esto supone 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

Se puede utilizar un filtro IIR especial. entonces el proceso que a medida como:

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

Creo que es más rápido.