Filtrage par fft-ifft

**zugeur** · 31/03/2009, 20h39

bonjour à tous,
Je suis entrain de faire un programme sous matlab qui utilise un enregistrement bruité d'une note de musique (LA3 = 440Hz). Je fais tout d'abord sa FFT puis voudrait uniquement garder une bande de fréquence de 430 à 450 Hz par exemple, dans le but de faire la ifft et d'écouter uniquement le son de ma note, sans le bruit qui la couvrait dans le signal de départ.

Voici mon code:

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
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
%IMPORTATION DE L'ENREGISTREMENT A ETUDIER
[x,fs,Nbits]=wavread('son wav');
%Fréquence d'échantillonnage
fs;
%Nombre de bits
Nbits;
%Longueur du signal "L" et nombre de colonne "C" de "x"
[L,C]=size(x);
%Définition de la Période
T=1/fs;
%Définition du vecteur temps:
t=(0:L-1)*T;
%dessiner le signal (fichier .wav)
figure(1);plot(t(1:L),x(1:L))
%limiter l'étude de 0 à 1 seconde
xlim([0 1])
title('Signal bruité enregistré')
xlabel('Temps (secondes)')
%TRANSFORMEE DE FOURIER (Fast Fourier Transform)
Y=fft(x,L)/L;
%Définition du vecteur fréquence
f=linspace(0,fs/2,floor(L/2)); %->Shannon (sinon repliement)
%dessiner la fft
figure(2);plot(f(1:floor(L/2)),abs(Y(1:floor(L/2))))
%limiter l'étude
xlim([20 4000])
title('Fast Fourier Transform de x(t)')
xlabel('Frequency (Hz)')
ylabel('|Y(f)|');grid
 
%FILTRE fft-ifft
%dessin de la bande du spectre à conserver
ymax=max(ylim);
hold on
line([0 430 430 450 450 L],[ymax ymax 0 0 ymax ymax],'color','r')
%filtrage
for i=1:L;
g(i)=Y(i)*heavi(abs(Y(i))-seuila(i,L,ymax));
end
figure(3)
i=1:L;
plot(i,g(i));grid
 
 
function s=seuila(n,N,ymax)
plage=[ymax 0;0.5e-3 430;ymax 450];
 
if n>=plage(1,2)
s=plage(1,1);
end
if n>=plage(2,2)
s=plage(2,1);
end
if n>=plage(3,2)
s=plage(3,1);
end
 
function v=heavi(x)
if x<=0
v=0;
else
v=1;
end
??? Undefined function or method 'seuila' for input arguments of type 'double'.

Mon problème, comme vous pouvez le voir, est que j'ai un message d'erreur lors de la définition de la fonction seuila apparement, qui définit la portion de fréquence à garder (430-450Hz)...

J'ai lu qu'il fallait peut être definir cette fonction et l'enregistrer dans un fichier .mat à part, mais où faut-il placer le fichier ? je l'ai mis dans le dossier toolbox de matlab mais sans succès, l'erreur est toujours présente...

Pour vous donner une idée du signal, regardez la pièce jointe représentant la fft du signal.

Merci de votre aide, je dois finir cela dans les plus brefs délais

**Jerome Briot** · 31/03/2009, 23h23

Le problème vient du fait que ton fichier m est un fichier script.
Et que tu déclares à l'intérieur de ce script une sous-fonction.

C'est impossible. Seul les fichiers fonctions peuvent contenir d'autres fonctions.

Il faut donc transformer ton script en fonction en ajoutant tout simplement cette ligne tout au déut du fichier :

Code :

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

function fun

En remplaçant fun par le nom du fichier sans l'extension .m

**zugeur** · 31/03/2009, 23h47

Salut, merci de ta réponse, j'avais déjà lu ca quelque part mais j'ai toujours une 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
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
>> function TIPE
%IMPORTATION DE L'ENREGISTREMENT A ETUDIER
[x,fs,Nbits]=wavread('son wav');
%Fréquence d'échantillonnage
fs;
%Nombre de bits
Nbits;
%Longueur du signal "L" et nombre de colonne "C" de "x"
[L,C]=size(x);
%Définition de la Période
T=1/fs;
%Définition du vecteur temps:
t=(0:L-1)*T;
%dessiner le signal (fichier .wav)
figure(1);plot(t(1:L),x(1:L))
%limiter l'étude de 0 à 1 seconde
xlim([0 1])
title('Signal bruité enregistré')
xlabel('Temps (secondes)')
%TRANSFORMEE DE FOURIER (Fast Fourier Transform)
Y=fft(x,L)/L;
%Définition du vecteur fréquence
f=linspace(0,fs/2,floor(L/2)); %Shannon (sinon repliement)
%dessiner la fft
figure(2);plot(f(1:floor(L/2)),abs(Y(1:floor(L/2))))
%limiter l'étude
xlim([20 4000])
title('Fast Fourier Transform de x(t)')
xlabel('Frequency (Hz)')
ylabel('|Y(f)|');grid
 
%FILTRE fft-ifft
%dessin de la bande du spectre à conserver
ymax=max(ylim);
hold on
line([0 430 430 450 450 L],[ymax ymax 0 0 ymax ymax],'color','r')
%filtrage
for i=1:L;
g(i)=Y(i)*heavi(abs(Y(i))-seuila(i,L,ymax));
end
figure(3)
i=1:L;
plot(i,g(i));grid
 
 
function s=seuila(n,N,ymax)
plage=[ymax 0;0.5e-3 430;ymax 450];
 
if n>=plage(1,2)
s=plage(1,1);
end
if n>=plage(2,2)
s=plage(2,1);
end
if n>=plage(3,2)
s=plage(3,1);
end
 
function v=heavi(x)
if x<=0
v=0;
else
v=1;
end
??? function TIPE
    |
Error: Function definitions are not permitted at the prompt or in scripts.

En remplaçant fun par le nom du fichier sans l'extension .m

Tu parles d'un nom de fichier aléatoire (celui que je veux lui donner) ou bien celui d'un fichier qu'il faudrait que j'enregistre au préalable ?
que faut-il faire ?
merci beaucoup

**Jerome Briot** · 31/03/2009, 23h50

Pourquoi fais-tu un copier-coller des lignes de code dans le Command Window ?

Enregistre simplement ton code dans le fichier m et tape ensuite le nom du fichier m sans l'extension dans le Command window.

Donc si ta fonction tipe est enregistrée dans tipe.m, tu tape simplement tipe après l'invite de commande >>

**zugeur** · 01/04/2009, 00h03

Dsl pour ces questions stupides mais je ne comprend pas, j'ai l'impression de ne pas savoir enregistrer le fichier correctement

.
juste pour confirmation, il faut bien que je tappe le code suivant dans Commande Window :

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
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
>> function TIPE
%IMPORTATION DE L'ENREGISTREMENT A ETUDIER
[x,fs,Nbits]=wavread('son wav');
%Fréquence d'échantillonnage
fs;
%Nombre de bits
Nbits;
%Longueur du signal "L" et nombre de colonne "C" de "x"
[L,C]=size(x);
%Définition de la Période
T=1/fs;
%Définition du vecteur temps:
t=(0:L-1)*T;
%dessiner le signal (fichier .wav)
figure(1);plot(t(1:L),x(1:L))
%limiter l'étude de 0 à 1 seconde
xlim([0 1])
title('Signal bruité enregistré')
xlabel('Temps (secondes)')
%TRANSFORMEE DE FOURIER (Fast Fourier Transform)
Y=fft(x,L)/L;
%Définition du vecteur fréquence
f=linspace(0,fs/2,floor(L/2)); %Shannon (sinon repliement)
%dessiner la fft
figure(2);plot(f(1:floor(L/2)),abs(Y(1:floor(L/2))))
%limiter l'étude
xlim([20 4000])
title('Fast Fourier Transform de x(t)')
xlabel('Frequency (Hz)')
ylabel('|Y(f)|');grid
 
%FILTRE fft-ifft
%dessin de la bande du spectre à conserver
ymax=max(ylim);
hold on
line([0 430 430 450 450 L],[ymax ymax 0 0 ymax ymax],'color','r')
%filtrage
for i=1:L;
g(i)=Y(i)*heavi(abs(Y(i))-seuila(i,L,ymax));
end
figure(3)
i=1:L;
plot(i,g(i));grid
 
 
function s=seuila(n,N,ymax)
plage=[ymax 0;0.5e-3 430;ymax 450];
 
if n>=plage(1,2)
s=plage(1,1);
end
if n>=plage(2,2)
s=plage(2,1);
end
if n>=plage(3,2)
s=plage(3,1);
end
 
function v=heavi(x)
if x<=0
v=0;
else
v=1;
end

Que j'enregistre ce fichier en faisant File>Save Workspace As...>tipe.mat (dans quel dossier dois-je l'enregistrer svp ?)

Puis que j'ouvre une nouvelle Command window dans laquelle je tappe directement après >> tipe ...

Est-ce bien cette démarche ou n'ai-je pas compris le principe

?

Merci encore

**zugeur** · 01/04/2009, 00h16

Ah non excusez moi, il faut que j'enregistre mon code comme étant un M-File bien sur !
Voici le code de mon fichier TIPE.m

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
46
function TIPE
%**********************************************************************************************************
%        DETECTION SONORE
%**********************************************************************************************************
%IMPORTATION DE L'ENREGISTREMENT A ETUDIER
[x,fs,Nbits]=wavread('C:\Users\Hugo\Etudes\ECAM\TIPE detection sonore\Sons Wav\diapason 440 + dehors CUT');
%Fréquence d'échantillonnage
fs;
%Nombre de bits
Nbits;
%Longueur du signal "L" et nombre de colonne "C" de "x"
[L,C]=size(x);
%Définition de la Période
T=1/fs;
%Définition du vecteur temps:
t=(0:L-1)*T;
%dessiner le signal (fichier Wav)
figure(1);plot(t(1:L),x(1:L))
%limiter l'étude de 0 à 1 seconde
xlim([0 1])
title('Signal bruité enregistré')
xlabel('Temps (secondes)')
%TRANSFORMEE DE FOURIER (Fast Fourier Transform)
Y=fft(x,L)/L;
%Définition du vecteur fréquence
f=linspace(0,fs/2,floor(L/2)); %Shannon !!! sinon repliement
%dessiner la fft
figure(2);plot(f(1:floor(L/2)),abs(Y(1:floor(L/2))))
%limiter l'étude
xlim([20 4000])
title('Fast Fourier Transform de x(t)')
xlabel('Frequency (Hz)')
ylabel('|Y(f)|');grid
 
%FILTRE fft-ifft
%dessin de la bande du spectre à conserver
ymax=max(ylim);
hold on
line([0 430 430 450 450 L],[ymax ymax 0 0 ymax ymax],'color','r')
%filtrage
for i=1:L;
g(i)=Y(i)*heavi(abs(Y(i))-seuila(i,L,ymax));
end
figure(3)
i=1:L;
plot(i,g(i));grid

Un nouveau message apparait lorsque je tappe TIPE dans la command window

Code :

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

1
2
3
4
>> TIPE
Warning: Imaginary parts of complex X and/or Y arguments ignored 
> In TIPE at 48
>>

Ligne 48 :

Code :

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

plot(i,g(i));grid

Comment y remédier ?
Merci

**phryte** · 01/04/2009, 14h10

Bonjour.
g est complexe :

Code :

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

plot(i,abs(g(i)));grid

**zugeur** · 01/04/2009, 19h15

Bonjour, merci pour vos réponses.

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
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
function TIPE %transformer le script en fonction
 
%**********************************************************************************************************
%        DETECTION SONORE
%**********************************************************************************************************
 
%IMPORTATION DE L'ENREGISTREMENT A ETUDIER
[x,fs,Nbits]=wavread('C:\Users\Hugo\Etudes\ECAM\TIPE detection sonore\Sons Wav\diapason 440 + dehors CUT');
%Fréquence d'échantillonnage
fs;
%Nombre de bits
Nbits;
%Longueur du signal "L" et nombre de colonne "C" de "x"
[L,C]=size(x);
%Définition de la Période
T=1/fs;
%Définition du vecteur temps:
t=(0:L-1)*T;
%dessiner le signal (fichier .wav)
figure(1)
plot(t(1:L),x(1:L))
%limiter l'étude de 0 à 1 seconde
xlim([0 1])
title('Signal bruité enregistré')
xlabel('Temps (secondes)')
 
%TRANSFORMEE DE FOURIER (Fast Fourier Transform)
Y=fft(x,L)/L;
%Définition du vecteur fréquence
f=linspace(0,fs/2,floor(L/2)); %Théorème de Shannon (sinon repliement)
%dessiner la fft
figure(2)
plot(f(1:floor(L/2)),abs(Y(1:floor(L/2))))
%limiter l'étude
xlim([20 4000])
title('Transformée de Fourier rapide de x(t)')
xlabel('Fréquence (Hz)')
ylabel('|Y(f)|');grid
 
 
%FILTRAGE (fft-ifft)
%dessin du gabarit (=du contour du spectre à conserver)
ymax=max(ylim);
hold on
line([0 430 430 450 450 L],[ymax ymax 0 0 ymax ymax],'color','r')
%filtrage
for i=1:L;
g(i)=Y(i)*heavi(abs(Y(i))-seuila(i,ymax));
end
figure(3)
i=1:L;
plot(i,abs(g(i)))
xlim([20 4000])
title('FFT de x(t) filtrée')
xlabel('Fréquence (Hz)')
ylabel('|Y(f)|');grid
 
%Transformée de Fourier inverse
figure(4)
iY=ifft(abs(g),L);
plot(t(1:L),real(iY(1:L)));
title('Son produit par le diapason')
xlabel('Temps (secondes)')
 
 
 
%Définition des fonctions
function s=seuila(n,ymax)
plage=[ymax 0;0 430;ymax 450];
 
if n>=plage(1,2)
s=plage(1,1);
end
if n>=plage(2,2)
s=plage(2,1);
end
if n>=plage(3,2)
s=plage(3,1);
end
 
function v=heavi(x)
if x<=0
v=0;
else
v=1;
end

-Je n'ai cette fois plus aucune erreur mais le signal (de la fft) plotté est symétrique par rapport à la droite verticale passant par 0,5sec (donc la moitier du signal), pourquoi ?
-De plus, je voulais savoir pourquoi devait-on mettre "real" absolument lors du plot de la ifft ?
-Enfin, je voulais savoir comment faire pour soit écouter le signal final (donc filtré, note de 440 Hz environ) soit/et extraire ce son dans un fichier .wav comme celui d'orignine ?

Merci beaucoup pour vos aides

**Jerome Briot** · 01/04/2009, 20h09

Envoyé par zugeur

-Enfin, je voulais savoir comment faire pour soit écouter le signal final (donc filtré, note de 440 Hz environ) soit/et extraire ce son dans un fichier .wav comme celui d'orignine ?

Code :

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

1
2
doc wavplay
doc wavwrite

**zugeur** · 01/04/2009, 20h46

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
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
function TIPE %transformer le script en fonction
 
%**********************************************************************************************************
%        DETECTION SONORE
%**********************************************************************************************************
 
%IMPORTATION DE L'ENREGISTREMENT A ETUDIER
[x,fs,Nbits]=wavread('C:\Users\Hugo\Etudes\ECAM\TIPE detection sonore\Sons Wav\diapason 440 + dehors CUT');
%Fréquence d'échantillonnage
fs;
%Nombre de bits
Nbits;
%Longueur du signal "L" et nombre de colonne "C" de "x"
[L,C]=size(x);
%Définition de la Période
T=1/fs;
%Définition du vecteur temps:
t=(0:L-1)*T;
%dessiner le signal (fichier .wav)
figure(1)
plot(t(1:L),x(1:L))
%limiter l'étude de 0 à 1 seconde
xlim([0 1])
title('Signal bruité enregistré')
xlabel('Temps (secondes)')
 
%TRANSFORMEE DE FOURIER (Fast Fourier Transform)
Y=fft(x,L)/L;
%Définition du vecteur fréquence
f=linspace(0,fs/2,floor(L/2)); %Théorème de Shannon (sinon repliement)
%dessiner la fft
figure(2)
plot(f(1:floor(L/2)),abs(Y(1:floor(L/2))))
%limiter l'étude
xlim([20 4000])
title('Transformée de Fourier rapide de x(t)')
xlabel('Fréquence (Hz)')
ylabel('|Y(f)|');grid
 
 
%FILTRAGE (fft-ifft)
%dessin du gabarit (=du contour du spectre à conserver)
ymax=max(ylim);
hold on
line([0 430 430 450 450 L],[ymax ymax 0 0 ymax ymax],'color','r')
%filtrage
for i=1:L;
g(i)=Y(i)*heavi(abs(Y(i))-seuila(i,ymax));
end
figure(3)
i=1:L;
plot(i,abs(g(i)))
xlim([20 4000])
title('FFT de x(t) filtrée')
xlabel('Fréquence (Hz)')
ylabel('|Y(f)|');grid
 
%Transformée de Fourier inverse
figure(4)
iY=ifft(abs(g),L);
plot(t(1:L),real(iY(1:L)));
xlim([0 1])
title('Son produit par le diapason')
xlabel('Temps (secondes)')
 
%écoute du signal filtré
wavplay(real(iY),fs)
wavwrite(real(iY),fs,'diapaSon fin')
 
 
%Définition des fonctions
function s=seuila(n,ymax)
plage=[ymax 0;0 430;ymax 450];
 
if n>=plage(1,2)
s=plage(1,1);
end
if n>=plage(2,2)
s=plage(2,1);
end
if n>=plage(3,2)
s=plage(3,1);
end
 
function v=heavi(x)
if x<=0
v=0;
else
v=1;
end

Merci pour ta réponse, l'enregistrement marche, je suppose que l'écoute aussi mais le problème est que je n'entend rien, peut être que l'amplitude du signal est trop faible ? comment faire pour l'amplifier pour pouvoir entendre le son ? (si la question se posait, j'ai testé l'écoute du signal de départ, donc tout va bien pour ce qui est de mes écouteurs

)

De plus, est-ce normal que le signal après ifft soit symétrique ? (voir pièce jointe) car ca ne devrait pas être le cas pour une note de musique ...

Merci encore

**Jerome Briot** · 01/04/2009, 21h27

Je ne suis pas un spécialiste mais je crois savoir que le volume d'un fichier wav est donné (entre autre) par l'amplitude du signal qui y est enregistré (l'amplitude variant entre -1 et 1)

On procède donc généralement à une normalisation du signal :

Code :

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

1
2
3
4
5
x = ...
 
x = x-min(x);
x = x/max(x);
x = 2*x-1;

Quoique je ne suis pas très sûr de moi la...

**zugeur** · 01/04/2009, 21h58

Je ne comprend vraiment pas la forme du signal obtenu (pièce jointe du dernier message), c'est censé etre la fréquence 440 uniquement ou preque, donc un son pur sinusoïdal ... si quelqu'un voit une erreur quelconque dans le filtrage, ou une explication à la symétrie qui est apparu ca m'aiderai beaucoup.

x = ...

x = x-min(x);
x = x/max(x);
x = 2*x-1;

Dut, je ne comprend pas vraiment ou tu veux en venir, que représente x ? pourquoi le soumet tu à plusieurs "définitions" comme cela ?

Comment amplifier un signal ?

Merci d'avance

**phryte** · 02/04/2009, 16h39

Bonjour.

est-ce normal que le signal après ifft soit symétrique ?

C'est le phénomène de repliement.
Tu dois étudier la fft entre 0 et fs/2 Hz (théorème de Shannon)

Pour le gain je te propose :

Code :

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

1
2
3
4
%Ajout de 20 dB
k=exp(20*2.3026/10)
signal=real(k*Y);
wavwrite(signal,fe,'LA')

**phryte** · 03/04/2009, 13h34

Bonjour.

est-ce normal que le signal après ifft soit symétrique

Fais un filtre un peu plus sélectif :

Code :

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

1
2
3
plage=[ymax 0;0 439.5;ymax 440.5];
 
line([0 439.5 439.5 440.5 440.5 L],[ymax ymax 0 0 ymax ymax],'color','r')

**zugeur** · 03/04/2009, 16h31

%Ajout de 20 dB
k=exp(20*2.3026/10)
signal=real(k*Y);
wavwrite(signal,fe,'LA')

Bonjour. Je n'ai pas réussi à résoudre mon problème de symétrie, mais par contre l'amplification du signal a fonctionnée, j'ai rajouté le code

Code :

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

1
2
3
k=exp(50*2.3026/10)
signal=real(k*Y);
wavwrite(signal,fe,'LA')

J'ai remplacé 20 par 50 car l'amplification par 20dB ne suffisait pas.Cela signifi-t-il que j'ai augmententer l'amplitude de 50dB ? je dis ca par analogie à ce que tu avais marqué... Peux tu m'expliquer exactement à quoi correspondent les differentes oppérations définissant k stp ?

Merci

Si vous savez quoi faire pour la symétrie, n'hésitez pas.

**phryte** · 03/04/2009, 17h17

Le gain en dB entre deux puissances P ou deux amplitudes U est :

Code :

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

X(dB) = 10*log(P1/P0) = 10*log(U1^2/U0^2)

Si P1=k*P0 :

Code :

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

X=10*log(k)

or log(x)=ln(x)/ln(10)=ln(x)/2.3026 et 2.3026*X/10=ln(k) ou k=exp(2.3026/10*X).
Si je veux 50 dB de plus : k =exp(2.3026/10*50).

Mais il faut raisonner en amplitude donc tu fais pareil avec :
Si U1=k*U0,

Code :

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

U1^2 = k^2*U0^2

donc il faut prendre la racine carrée de k si tu fais le calcul que je t'ai indiqué.
http://fr.wikipedia.org/wiki/D%C3%A9cibel

Filtrage par fft-ifft

Signal

Discussions similaires

Partager

Partager