Discussion :

[tiop5]

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 :

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

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.

[phryte]

Bonjour,

Un essai :

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
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

[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?

[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.

[phryte]

Bonjour,

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

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
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

[tiop5]

Bonjour,
Pouvez-vous m'expliquer ce que vous faites dans votre algorithme svp ?

[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.

[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.

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 :

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

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

[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.

[tiop5]

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

Code :

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

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

et pas :

Code :

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

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 ?

[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.

[tiop5]

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 :

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
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

[phryte]

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

Code :

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

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

[phryte]

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
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

[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 ?

[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

[phryte]

Bonjour,

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

Pour la solution en testant les amplitudes, un essai :

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
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

[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

[phryte]

Bonjour,

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

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
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 :

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
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

[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?