Discussion :

[vampirella]

Bonjour à tous,

Le titre a l'air barbare, mais le vrai problème est de pouvoir enlever la boucle "for" de l'algorithme qui est présenté ci-dessous :

Code :

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

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
function estimated = kernel_smoothing( time, diracVect, Sf, KW)
 
Ts = 1/Sf;    %Sampling Time
estimated = zeros(length(time), 1);
N = length(time);
n = 1:N;
%check that time starts from 0 second
rescaled_time = time - time(1);
 
for k = 1:N
 
    inside_b = abs( k*Ts - rescaled_time(n) )/KW;
    b_res = Gaussian_kernel_func(inside_b);
 
    numerator = sum(diracVect .* b_res);
    denominator = sum(b_res);
 
    estimated(k) = numerator/denominator;
end
 
end

En effet, plus le signal 1D (fourni ici dans le paramètre "diracVect") à traiter est grand, plus l'algorithme prends du temps à s'exécuter, ce qui est normal ... Mais cela peut atteindre plusieurs dizaines de minutes pour moins de 10 minutes de signal !!

J'ai essayé plusieurs tentatives depuis hier, mais aucune ne m'a apporté satisfaction (au niveau temps ou résultat final).

Ca serait pas mal si quelqu'un pouvait ne serait-ce que m'aiguiller vers une bonne solution


PS : pour l'histoire, la fonction gaussienne est la suivante :

Code :

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

1
2
3
4
5
6
7
8
9
10
11
12
13
function r = Gaussian_kernel_func(u)
%
%        [ exp((-u^2)/2))    if -5 <= u <= 5
% b(u) = {
%        [ 0    otherwise
%
 
in_bound = find( (-5 <= u) & (u <= 5));
 
r = zeros(size(u));
 
r(in_bound) = exp((-u(in_bound).^2)/2);
end

[duf42]

Bonjour,

A première vue j'ai comme l'impression que ta boucle ne sert à rien ou plutôt que tu effectue N fois le même calcul (ton N n'est pas utilisé dans la boucle)

Il y a sans doute un truc que je n'ai pas vu mais as-tu essayé en enlevant purement et simplement la boucle.

Duf

[vampirella]

Bonjour duf42,

En effet, N n'est pas réutilisé tel quel, et ne sert qu'à marquer l'indice du dernier point.

Pour résumer l'action de la fonction, chaque nouveau point k du nouveau signal estimé ("estimated") est calculé par un moyennage gaussien du signal fourni en entrée ("diracVect") dans les alentours de k.

Ainsi, ce nouveau point en k dépendra des valeurs à environ une seconde avant et après ce même point k du signal d'entrée.

J'ai donc besoin de la boucle (dans l'état actuel des choses) pour calculer, pas à pas, les nouveaux points du signal estimé. Mais comme la boucle me prends trop de temps ...


* Je ne sais pas si je suis suffisamment clair, demandez-moi plus de précisions si besoin *

[duf42]

Ah oui autant pour moi, comme d'habitude j'ai lu trop vite...

Aurais-tu un échantillon de données pour tester?

[vampirella]

Bien sûr

Compter plus de 20 secondes de calcul pour 22k5 points échantillonés ... il faut être patient


EDIT : La pause café m'a donné une idée, qui serait d'utiliser "meshgrid" de la façon suivante :

Code :

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

1
2
3
4
5
6
7
8
9
10
11
12
%check that time starts from 0 second
rescaled_time = time - time(1);
 
[meshG_time, meshG_kTs] = meshgrid(n.*Ts, rescaled_time(n));
 
inside_b = abs( meshG_kTs - meshG_time )/KW;
b_res = Gaussian_kernel_func(inside_b);
 
numerator = sum(diracVect * b_res, 2);
denominator = sum(b_res, 2);
 
estimated = numerator/denominator;

Je pense que cette méthode devrait résoudre mon problème de temps. Cependant, MatLab me fournit une erreur de mémoire (alors que je suis en version 7.5 et que je ne dépasse pas 2^31 en taille , comme préconisé ici !!)

Code :

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

1
2
3
Maximum variable size allowed by the program is exceeded.
 
Error in ==> meshgrid at 44

[Jean Dumoncel]

Salut,

Cependant, MatLab me fournit une erreur de mémoire (alors que je suis en version 7.5 et que je ne dépasse pas 2^31 en taille , comme préconisé ici !!)

Code :

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

1
2
3
Maximum variable size allowed by the program is exceeded.
 
Error in ==> meshgrid at 44

C'est spécifié que tu peux atteindre 2^31 éléments sur un système 64 bit.

Avec meshgrid tu cherche à créer des données de tailles :
25500*25500*8 = environ 4Gb ce qui est au-dessus de la capacité des systèmes 32 bit.

[vampirella]

En effet, cela fonctionne encore un peu au ralenti, mais cela donne quelque chose.

Pour la fin mot de l'histoire, j'ai trouvé une autre implémentation de kernel smoothing qui, je doit bien l'avouer, est plus élégante que tout ce que je pourrais essayer de bricoler ...

J'ai tout de même repris ton idée, Magelan, et avec des calculs de kernel sur 2*Fe points (Fe = Fréquence d'échantillonnage), je tombe à moins de 15 secondes de calculs pour 25 minutes d'enregistrements à 125 Hz.

C'est résolu pour ma part.