Intercorrélation de deux signaux

Version imprimable

Voir 40 message(s) de cette discussion en une page

Bonsoir à tous,
Je dois coder un programme sous matlab qui me permettrait de comparer les fréquences de deux sons audio purs (donc de forme Acos(2*pi*f*t)), que j'ai préalablement générés. Cette partie est correcte d'après mon professeur.
On connait la fréquence d’échantillonnage, qui est de 16KHz. L'algorithme doit renvoyer un résultat binaire : 1 si les deux fréquences sont les mêmes, à une marge d’erreur près, 0 sinon.
Je dois pour cela utiliser l'intercorrélation dont la formule est ci-dessous. Voilà ce que j'ai écrit :
Code:

1 2 3 4 5 6 7 8 9 10 11 12 N = length(x1); Cxx = zeros (1,2*N); inter = zeros(1,N); for k = 1 : N for n = k : N-1 inter(n) = x1(n)*x2(n+k); end Cxx(k) = 1/N * sum(inter); end plot(Cxx) ;
Cependant, on m'affiche l'erreur : Index exceeds the number of array elements (501).
Error in Exercice (line 43)
inter(n) = x1(n)*x2(n+k);
Si je comprend bien, je pense que je n’ai pas bien indicé la boucle for, mais je n’arrive pas à trouver comment corriger ce problème.
Merci d'avance pour votre aide.

Bonjour,

Un essai :

Code:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
 clear
Te=1/12000;
A=1;f1=200;f2=200;
t=0:Te:0.01;
x1=A*cos(2*pi*f1*t);
x2=A*cos(2*pi*f2*t);
 
N = length(x1);
Cxx = zeros(1,2*N);
inter = zeros(1,N);
[c,lags]=xcorr(x1,x2);
figure(1)
plot(lags,c)
grid
for k = 1: N-1
    for n = 1 : N-k-1
        inter(n) = x1(n)*x2(n+k);
    end
    Cxx(k) = 1/N*sum(inter);
end
figure(2)
plot(Cxx) ;
grid
 
figure(3)
plot(lags,c/max(c))
hold on
plot(Cxx/max(Cxx),'r')
grid

06/01/2021, 16h49
tiop5

Bonjour,
Merci pour votre réponse !
Pouvez-vous m'expliquer pourquoi avoir utilisé la fonction xcorr, et comment puis-je renvoyer une réponse binaire par rapport à une marge d'erreur?
06/01/2021, 16h57
phryte

Bonjour,
J’ai utilisé xcorr pour comparer les résultats de matlab avec ton programme (voir figure 3)
Je regarde la suite.

Tu peux regarder la discussion ci-dessous qui traite le même sujet où on trouve un algorithme de mesure de période du signal.

Bonjour,

Un essai avec le test sur la période de la fonction d'intercorrélation :

Code:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
clear
Te=1/12000;
A=1;f1=200;f2=210;
t=0:Te:0.02;
x1=A*cos(2*pi*f1*t);
x2=A*cos(2*pi*f2*t);
N = length(x1);
for k = 1: N-1
    for n = 1 : N-k-1
        inter(n) = x1(n)*x2(n+k);
    end
    Cxx(k) = 1/N*sum(inter);
end
figure(1)
plot((1:N-1)*Te,Cxx) ;
grid
d=find(Cxx(1:end-1).*Cxx(2:end)<0);
dd=diff(d(2:end-1))
%dd=dd(1:end-1)
%diff(dd)
% Résultat binaire
if max(dd)-min(dd) > 2
    disp('signaux dfifférents')
else
    disp('signaux prochent')
end

07/01/2021, 12h08
tiop5

Bonjour,
Pouvez-vous m'expliquer ce que vous faites dans votre algorithme svp ?
07/01/2021, 12h46
phryte

Bonjour,
Rappel :
Propriété : si deux signaux x et y sont périodiques de même période T0 alors leur intercorrélation
Rxy(τ) est périodique de période T0.
Deux signaux x et y seront indépendants s’il n’y a pas de lien entre eux. S’il existe un lien alors celui-ci doit dépendre du décalage temporel entre ces deux signaux.

Je mesure les périodes du signal (écarts en zéros) s’il n’y a pas de variations importantes, les signaux ont des fréquences proches.
07/01/2021, 13h09
tiop5
Ah je comprend mieux comme ça merci !
Deux signaux x et y seront indépendants s’il n’y a pas de lien entre eux. S’il existe un lien alors celui-ci doit dépendre du décalage temporel entre ces deux signaux.

Citation:

Je mesure les périodes du signal (écarts en zéros) s’il n’y a pas de variations importantes, les signaux ont des fréquences proches.

Je ne comprend pas, comment vous avez fait ce que vous dites dans le code que vous m'avez envoyé ?
Code:

1 2 3 4 plot((1:N-1)*Te,Cxx) ; grid d=find(Cxx(1:end-1).*Cxx(2:end)<0); dd=diff(d(2:end-1))
En fait, je ne comprend pas la manière dont vous avez indicé, et ce que représentent les points détectés à la ligne 3
07/01/2021, 13h29
phryte

La ligne 3 détecte le passage à zéro du signal par le changement de signe de l’instant n-1 à l’instant n.
07/01/2021, 13h36
tiop5

C'est vraiment une bonne idée merci beaucoup !
Pour la ligne 4, pourquoi avoir écrit :

Code:

dd=diff(d(2:end-1))

et pas :

Code:

dd=diff(d)

J'ai vu qu'on obtenait deux valeurs de plus dans le vecteur dd, mais est ce que cela change quelque chose dans les résultats ?
07/01/2021, 14h21
phryte

diff(d(2:end-1)) a pour but d’éliminer les mesures extrêmes. Selon le pas choisi il peut y avoir des erreurs.
D’autre part, si les signaux sont bruités il faudra les filtrer.

En fait, voilà mon code final pour générer deux signaux sinusoidaux purs, puis vérifier s'ils ont la même fréquence. Pour faire le test, j'ai écrit que f1=f2, mais avec ce code, ça renvoie que les fréquences sont différentes... Je pense que l'erreur peut venir du fait que le signal x2 a une phase non nulle (c'est obligatoire dans l'exercice, x1 a une phase nulle, x2 une phase non nulle)
Code:

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 f1 = 440 ; f2 = 440 ; Fe = 16000; Te = 1/Fe; T1 = 1/f1; T2 = 1/f2; phi2 = 2*pi*rand(1); A1 = randi([1,10]); A2 = randi([1,10]); Nb_Ech = 500; alpha = (Nb_Ech*f1)/Fe ; D = alpha /f1 ; t = 0:Te:D; x1 = A1*cos(2*pi*f1*t); x2 = A2*cos(2*pi*f2*t + phi2); subplot(3,1,1) plot(t,x1, 'b') subplot(3,1,2) plot(t, x2, 'r') N = length(x1); for k = 1: N-1 for n = 1 : N-k-1 inter(n) = x1(n)*x2(n+k); end Cxx(k) = 1/N*sum(inter); end subplot(3,1,3) plot((1:N-1)*Te,Cxx) grid d=find(Cxx(1:end-1).*Cxx(2:end)<0); disp(d); dd = diff(d(2:end-1)) ; %dd=dd(1:end-1) %diff(dd) % Résultat binaire if max(dd)-min(dd) > 2 disp('signaux différents') else disp('signaux prochent') end

Il faut éliminer plus de mesures de la fin du signal qui tend vers zéro:

Code:

1
2
3
4
5
6
7
8
9
10
d=find(Cxx(1:end-1).*Cxx(2:end)<0);
dd = diff(d) ;
dd=dd(2:end-5);
disp(dd);
% Résultat binaire
if max(dd)-min(dd) > 3
    disp('signaux différents')
else
    disp('signaux prochent')
end

Code:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
clear
f1 = 440 ;
f2 = 460 ;
Fe = 16000;
Te = 1/Fe;
T1 = 1/f1;
T2 = 1/f2;
phi2 = 2*pi*rand(1);
A1 = randi([1,10]);
A2 = randi([1,10]);
Nb_Ech = 500;
alpha = (Nb_Ech*f1)/Fe ;
D = alpha /f1 ;
t = 0:Te:D;
x1 = A1*cos(2*pi*f1*t);
x2 = A2*cos(2*pi*f2*t + phi2);
subplot(3,1,1) 
plot(t,x1, 'b')
grid
subplot(3,1,2)
plot(t, x2, 'r')
grid
N = length(x1);
for k = 1: N-1
    for n = 1 : N-k-1
        inter(n) = x1(n)*x2(n+k);
    end
    Cxx(k) = 1/N*sum(inter);
end
subplot(3,1,3)
plot((1:N-1)*Te,Cxx) 
grid
d=find(Cxx(1:end-1).*Cxx(2:end)<0);
dd = diff(d) ;
dd=dd(1:end-5);
disp(dd);
% Résultat binaire
if max(dd)-min(dd) > 3
    cprintf('red','signaux différents')
else
    cprintf('blue','signaux prochent')
end

11/01/2021, 10h14
tiop5

Bonjour,
Je n'avais pas vu votre réponse, désolée de ne répondre que maintenant...
J'ai essayé de faire tourner le programme sur matlab, mais il y a une erreur à cause de "cprintf". Je n'ai jamais utilisé cette commande, donc je ne sais pas d'ou cette erreur peut provenir.
Aussi, mon professeur a dit qu'il serait surement mieux d'utiliser le maximum de l'intercorrélation. En effet, on sait que le maximum de la fonction d’intercorrélation indique le décalage
temporel entre deux signaux.
J'avais aussi une question sur la facon dont vous avez indicé la boucle qui correspond à la formule de l'intercorrélation, pourquoi avoir mis n=0:N-k-1, alors que dans la formule, n=0:N-1 ?
11/01/2021, 10h23
tiop5

Aussi, lorsque vous avez comparé avec xcorr, on voit bien que la fonction n'est pas représentée graphiquement en entière, mais que à partir de 0

Bonjour,

cprintf ne sert qu’à afficher un résultat en couleur.

Pour la solution en testant les amplitudes, un essai :

Code:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
clear
Te=1/12000;
A=1;f1=440;f2=450;
phi2=2*pi*rand(1)
t=0:Te:0.02;
A1=randi([1,10]);
A2=randi([1,10]);
x1=A1*cos(2*pi*f1*t);
x2=A2*cos(2*pi*f2*t+phi2);
% Normalisation
x1=A1*cos(2*pi*f1*t)/max(x1);
x2=A2*cos(2*pi*f2*t+phi2)/max(x2);
N = length(x1);
% Test du maximum
for k = 1: N-1
    for n = 1 : N-k-1
        inter(n) = x1(n)*x1(n+k);
    end
    Cxx(k) = 1/N*sum(inter);
end
% Maximum de référence
MaxCxx=max(Cxx)
for k = 1: N-1
    for n = 1 : N-k-1
        inter(n) = x1(n)*x2(n+k);
    end
    Cxx(k) = 1/N*sum(inter);
end
% Test des maxima - Seuil à ajuster
if MaxCxx - max(Cxx) > 0.1
 
    disp('signaux dfifférents')
else
    disp('signaux prochent')
end
 
figure(1)
plot(t,x1,t,x2)
grid
figure(2)
plot((1:N-1)*Te,Cxx) ;
grid

12/01/2021, 12h46
tiop5

Bonjour,
J'ai testé votre programme mais je crains qu'il ne fonctionne pas correctement... En effet, en choisissant f1 = f2 = 500 Hz, le programme retourne parfois que les signaux sont proches, et d'autres fois que les signaux sont différents, or il devrait toujours retourner que les fréquences sont proches

Bonjour,

On peut durcir l'algorithme en limitant les signaux à un nombre d'alternances entier :

Code:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
clear
Te=1/12000;
A=1;f1=500;f2=500;
phi2=2*pi*rand(1);
t=0:Te:0.02;
A1=randi([1,10]);
A2=randi([1,10]);
x1=A1*cos(2*pi*f1*t);
x2=A2*cos(2*pi*f2*t+phi2);
% Normalisation
x1=A1*cos(2*pi*f1*t)/max(x1);
x2=A2*cos(2*pi*f2*t+phi2)/max(x2);
N = length(x1);
% Test du maximum
for k = 1: N-1
    for n = 1 : N-k-1
        inter(n) = x1(n)*x1(n+k);
    end
    Cxx(k) = 1/N*sum(inter);
end
% Recherche des maxima
 MM=find(x1==1);
% Maximum de référence
MaxCxx=max(Cxx);
for k = 1: N-1
    for n = 1 : N-k-1
        inter(n) = x1(n)*x2(n+k);
    end
    S=inter(1:MM(end-1));% Troncature pour alternances entières
    SL=length(S);
    Cxx(k) = 1/(SL)*sum(S);
end
% Test des maxima - Seuil à ajuster
if MaxCxx - max(Cxx) > 0.1
    disp('signaux dfifférents')
else
    disp('signaux prochent')
end
 
figure(1)
plot(t,x1,t,x2)
grid
figure(2)
plot((1:N-1)*Te,Cxx) ;
grid

L'algorithme en testant les périodes fonctionne mieux ?

J'ai refais un essai avec cette version et je n'ai pas ton erreur :

Code:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
clear
Te=1/12000;
A=1;f1=500;f2=500;
phi2=2*pi*rand(1);
t=0:Te:0.02;
A1=randi([1,10]);
A2=randi([1,10]);
x1=A1*cos(2*pi*f1*t);
x2=A2*cos(2*pi*f2*t+phi2);
% Normalisation
x1=A1*cos(2*pi*f1*t)/max(x1);
x2=A2*cos(2*pi*f2*t+phi2)/max(x2);
N = length(x1);
% Test du maximum
for k = 1: N-1
    for n = 1 : N-k-1
        inter(n) = x1(n)*x1(n+k);
    end
    Cxx(k) = 1/N*sum(inter);
end
% Maximum de référence
MaxCxx=max(Cxx);
for k = 1: N-1
    for n = 1 : N-k-1
        inter(n) = x1(n)*x2(n+k);
    end
    Cxx(k) = 1/N*sum(inter(1:end-3));
end
% Test des maxima - Seuil à ajuster
if MaxCxx - max(Cxx) > 0.1
    disp('signaux dfifférents')
else
    disp('signaux prochent')
end
 
figure(1)
plot(t,x1,t,x2)
grid
figure(2)
plot((1:N-1)*Te,Cxx) ;
grid

12/01/2021, 20h22
tiop5

Ah d'accord, en fait c'est quand f1 = f2 = 440 Hz que ça ne fonctionne pas, je ne comprend pas pourquoi d'ailleurs... Mais sinon ca marche merci ! Il me reste juste le problème de centrer l'intercorrélation en 0 comme il y a dans mon cours, avez vous une idée?

Voir 40 message(s) de cette discussion en une page