Accélérer les calculs FFT clairsemés

Question

J'espère que quelqu'un pourra revoir mon code ci-dessous et vous donner des conseils pour accélérer la section entre tic et toc. La fonction ci-dessous tente d'effectuer un IFFT plus rapidement que la fonction intégrée de Matlab puisque (1) presque tous les bots de coefficients fft sont nuls (c.-à-d. 10 à 1000 les bins sur 10M à 300M sont non nulles), et (2) seulement le tiers central des résultats de l'IFFT est conservé (le premier et le dernier tiers sont rejetés - donc pas besoin de les calculer en premier lieu).

Les variables d'entrée sont les suivants:

fftcoef = complex fft-coef 1D array (10 to 1000 pts long) 
bins = index of fft coefficients corresponding to fftcoef (10 to 1000 pts long) 
DATAn = # of pts in data before zero padding and fft (in range of 10M to 260M) 
FFTn = DATAn + # of pts used to zero pad before taking fft (in range of 16M to 268M) (e.g. FFTn = 2^nextpow2(DATAn)) 

Actuellement, ce code prend quelques ordres de grandeur plus que l'approche de la fonction de Matlab ifft qui calcule l'ensemble du spectre élimine alors 2/3 « s de celui-ci. Par exemple, si les données d'entrée pour fftcoef et les bacs sont 9x1 tableaux (par exemple seulement 9 coefficients fft complexes par bande latérale, 18 pts lorsque l'on considère les deux bandes latérales), et DATAn=32781534, FFTn=33554432 (c.-à-2^25), l'approche IFFT prend 1.6 secondes alors que la boucle ci-dessous prend plus de 700 secondes.

J'ai évité d'utiliser une matrice pour vectoriser la boucle nn car parfois la taille du tableau pour fftcoef et bacs pourrait être jusqu'à 1000 pts long, et une matrice 260Mx1K serait trop grand pour la mémoire à moins qu'il puisse être divisé en quelque sorte .

Un conseil est très apprécié! Merci d'avance.

function fn_fft_v1p0(fftcoef, bins, DATAn, FFTn) 

fftcoef = [fftcoef; (conj(flipud(fftcoef)))];  % fft coefficients 
bins = [bins; (FFTn - flipud(bins) +2)];   % corresponding fft indices for fftcoef array 

ttrend = zeros((round(2*DATAn/3) - round(DATAn/3) + 1), 1); % preallocate 

start = round(DATAn/3)-1; 

tic; 
for nn = start+1 : round(2*DATAn/3) % loop over desired time indices 
    % sum over all fft indices having non-zero coefficients 
    arg = 2*pi*(bins-1)*(nn-1)/FFTn; 
    ttrend(nn-start) = sum(fftcoef.*(cos(arg) + 1j*sin(arg)); 
end 
toc; 

end 

Answer 1

Vous devez garder à l'esprit que Matlab utilise une bibliothèque compilé fft (http://www.fftw.org/) pour ses fonctions fft qui, en plus d'exploitation beaucoup plus rapide puis un script Matlab, il est bien optimisé pour de nombreux cas d'utilisation. Donc, une première étape pourrait être d'écrire votre code dans c/C++ et de le compiler sous forme de fichier mex que vous pouvez utiliser dans Matlab. Cela va certainement accélérer votre code au moins un ordre de grandeur (probablement plus).

En outre, une optimisation simple que vous pouvez faire est en considérant 2 choses:

1. Vous assumez votre série temporelle est une valeur réelle, vous pouvez donc utiliser la symétrie des coeffs fft.
2. Votre série temporelle est généralement beaucoup plus longue que votre vecteur coeffs fft, il est donc préférable d'itérer sur les cases au lieu des points temporels (vectorisant ainsi le vecteur le plus long).

Ces deux points sont convertis à la boucle suivante:

nn=(start+1 : round(2*DATAn/3))'; 
ttrend2 = zeros((round(2*DATAn/3) - round(DATAn/3) + 1), 1); 
tic; 
for bn = 1:length(bins) 
    arg = 2*pi*(bins(bn)-1)*(nn-1)/FFTn; 
    ttrend2 = ttrend2 + 2*real(fftcoef(bn) * exp(i*arg)); 
end 
toc; 

Notez que vous devez utiliser cette boucle avant vous étendez bins et fftcoef, puisque la symétrie est déjà pris en compte. Cette boucle prend 8,3 secondes pour s'exécuter avec les paramètres de votre question, alors que cela prend mon ordinateur 141,3 secondes pour fonctionner avec votre code.

Answer 2

J'ai posté une question/réponse à Accelerating FFTW pruning to avoid massive zero padding qui résout le problème pour le cas C++ en utilisant FFTW. Vous pouvez utiliser cette solution en exploitant les fichiers mex.