Discussion :

[Gébix]

Bonjour à tous,

j'ai un problème de temps d'exécution avec matlab. En effet, j'ai réécrit une fonction FFT mais celle-ci est beaucoup plus lente que celle implémenté dans matlab. Y a-t-il un moyen d'accélérer tout ça ?

Au cas où il y aurai un problème dans mon code, je vous donne le code de ma FFT :

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
10
11
N = length(signal);
spectrum = zeros(N, 1);
 
for k=0:ceil((N+1)/2)
	for n=0:N-1
		spectrum(k+1) = spectrum(k+1) + signal(n+1)*exp(-2*pi*j*k*n/N);
	end
	if k>0
		spectrum(N-k+1) = conj(spectrum(k+1));
	end
end

[Jerome Briot]

La FFT sous MATLAB est implémentée grâce à la bibliothèque FFTW écrite en C. Je doute qu'une implémentation en MATLAB puisse rivaliser...

[Gébix]

L'utilisation d'un mex-file pourrait-il grandement arrangé mon problème ? (j'imagine que oui mais c'est juste pour confirmation vu que j'ai encore jamais utilisé un tel fichier)

[Jerome Briot]

Oui... en théorie.

[Gébix]

okai, vais voir ça. Merci pour cette prompte réponse

[shaiHulud]

Tu peux au moins vectoriser

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
for n=0:N-1
		spectrum(k+1) = spectrum(k+1) + signal(n+1)*exp(-2*pi*j*k*n/N);
	end

en

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
 
spectrum(k+1)= sum(signal.*exp(-2*pi*j*k*(0:N-1)/N);

et essayer de voir si tu peux vectoriser la seconde boucle (ca m'a l'air tout à fait faisable)

Ca sera surement mieux, mais dans tous les cas, tu n'atteindra pas la vitesse C

[Jerome Briot]

Il est possible de pousser un peu plus loin l'optimisation en remplaçant la fonction SUM par une simple multiplication matricielle :

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
N = length(signal);
spectrum = zeros(N, 1);
 
spectrum(1)= sum(signal);
for k=1:ceil((N+1)/2)
    spectrum(k+1)= signal*exp(-2*pi*1i*k*(0:N-1)/N)';
    spectrum(N-k+1) = conj(spectrum(k+1));
end

ou encore

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
N = length(signal);
spectrum = zeros(N, 1);
 
spectrum(1)= sum(signal);
for k=1:ceil((N+1)/2)
    spectrum(k+1)= exp(-2*pi*1i*k*(0:N-1)/N)*signal(:);
    spectrum(N-k+1) = conj(spectrum(k+1));
end

Attention, SPECTRUM est le nom d'une fonction de la Signal Processing Toolbox. Il serait préférable de modifier le nom de cette variable.