Я хотел бы сделать линейную свертку для сигнала длиной 4000 * 270 с ядром длиной 16000. Сигнал не фиксируется, пока ядро фиксировано. Это нужно повторять много раз для моей цели, поэтому я хочу улучшить скорость как можно скорее. Я могу реализовать эту свертку либо в R или C.
Сначала я попытался выполнить свертку в R, но скорость не может удовлетворить мою потребность. Я пытался делать это итерацией, и это было слишком медленно. Я также попытался сделать это с помощью FFT, но поскольку сигнал и ядро длинные, FFT не сильно улучшил скорость.
Тогда я решил итеративно выполнить свертку в Си. Но С, похоже, не в состоянии обрабатывать такое количество вычислений и сообщать об ошибке очень часто. Даже когда это работает, это все еще очень медленно. Я также пытался сделать свертку FFT в C, но программа всегда закрывалась.
Я нашел этот код у моего друга и не уверен насчет оригинального источника. Я удалю его, если возникнет проблема с авторскими правами. Это код C, который я использовал для выполнения fft в C, но программа не может обработать длинный вектор длиной 2097152 (наименьшая степень 2 больше или равна длине вектора сигнала ).
#define q 3 /* for 2^3 points */
#define N 2097152 /* N-point FFT, iFFT */
typedef float real;
typedef struct{real Re; real Im;} complex;
#ifndef PI
# define PI 3.14159265358979323846264338327950288
#endif
void fft( complex *v, int n, complex *tmp )
{
if(n>1) { /* otherwise, do nothing and return */
int k,m;
complex z, w, *vo, *ve;
ve = tmp;
vo = tmp+n/2;
for(k=0; k<n/2; k++) {
ve[k] = v[2*k];
vo[k] = v[2*k+1];
}
fft( ve, n/2, v ); /* FFT on even-indexed elements of v[] */
fft( vo, n/2, v ); /* FFT on odd-indexed elements of v[] */
for(m=0; m<n/2; m++) {
w.Re = cos(2*PI*m/(double)n);
w.Im = -sin(2*PI*m/(double)n);
z.Re = w.Re*vo[m].Re - w.Im*vo[m].Im; /* Re(w*vo[m]) */
z.Im = w.Re*vo[m].Im + w.Im*vo[m].Re; /* Im(w*vo[m]) */
v[ m ].Re = ve[m].Re + z.Re;
v[ m ].Im = ve[m].Im + z.Im;
v[m+n/2].Re = ve[m].Re - z.Re;
v[m+n/2].Im = ve[m].Im - z.Im;
}
}
return;
}
void ifft( complex *v, int n, complex *tmp )
{
if(n>1) { /* otherwise, do nothing and return */
int k,m;
complex z, w, *vo, *ve;
ve = tmp;
vo = tmp+n/2;
for(k=0; k<n/2; k++) {
ve[k] = v[2*k];
vo[k] = v[2*k+1];
}
ifft( ve, n/2, v ); /* FFT on even-indexed elements of v[] */
ifft( vo, n/2, v ); /* FFT on odd-indexed elements of v[] */
for(m=0; m<n/2; m++) {
w.Re = cos(2*PI*m/(double)n);
w.Im = sin(2*PI*m/(double)n);
z.Re = w.Re*vo[m].Re - w.Im*vo[m].Im; /* Re(w*vo[m]) */
z.Im = w.Re*vo[m].Im + w.Im*vo[m].Re; /* Im(w*vo[m]) */
v[ m ].Re = ve[m].Re + z.Re;
v[ m ].Im = ve[m].Im + z.Im;
v[m+n/2].Re = ve[m].Re - z.Re;
v[m+n/2].Im = ve[m].Im - z.Im;
}
}
return;
}
Я нашел эту страницу, говорящую о свертке длинных сигналов https://ccrma.stanford.edu/~jos/sasp/Convolving_Long_Signals.html
Но я не уверен, как использовать идею в этом. Любые мысли будут по достоинству оценены, и я готов предоставить больше информации по моему вопросу.