butter, filter et fft

Version imprimable

27/09/2012, 23h19
membreComplexe12

butter, filter et fft
Salut tous :zoubi:,
je voudrais enfin comprendre comment on filtre une certaine bande de fréquence d'un signal, disons par exemple aux alentours de 50Hz. J'ai regardé pas mal de discussions récentes et voici ce que j'ai identifié comme méthode :roll: :

1°) faire une "fft" et afficher le spectre. Supprimer les fréquences qui ne nous intéressent pas en mettant à 0 les coeff correspondant à ces fréquences.et faire une "ifft"
2°) la deuxième chose que j'ai identifié c'est un filtre coupe bande.
Pour faire ceci il semble y avoir deux solutions "filter" et "fir1"
3°) la troisième solution qui semble être possible est la fonction "butter"

Maintenant identifier les méthodes c'est bien mais mettre en place est une autre :cry: Voici ce que je fais :
1°) pour la fft j'ai un peu de mal à comprendre comment faire, voici ce que j'ai fais :
Code:

1 2 x=data(:,1); y=data(:,2); yFFT=fft(y,length(x));
le soucis est que jusqu'à présent je n'ai toujours pas compris comment déduire le vecteur de fréquences pour tracer mon spectre...
=> je me dis que la fréquence fondamentale est x(end)-x(1) ? mais comment déduire les autres fréquences ? (surtout que ça devrait plutôt être une données de sortie de fft non ?)

Ensuite, si je suppose que j'ai fais mon spectre. Je mets manuellement les coefficients de la variable yFFT qui ne m'interesse pas à 0 et je fais :
Code:

1 2 3 yTemporel=ifft(yFFT); plot(x,yTemporel)
et ça devrait fonctionner non ? :oops:

2°) Pour le filtre coupe bande j'ai une première méthode qui consiste à faire :
Code:

1 2 3 4 x=data(:,1); y=data(:,2); frequenceBasse=40; frequenceHaute=60; ordreFiltre=1; yFiltre = fir1(1,[2*pi*frequenceBasse 2*pi*frequenceHaute],'stop');
ceci ça devrait être pas mal, non ? :oops: mais par contre je n'ai pas compris à quoi ça sert d'avoir un filtre d'ordre supérieur à 1 ??? j'aurais tendance à dire que ça filtre de plus en plus si on augmente l'ordre car les pentes dans le diagramme de bode seront plus pentues ?

la deuxieme méthode consiste à utiliser directement "filter"
Code:

1 2 3 x=data(:,1); y=data(:,2); a=[];b=[] yFiltre = filter(b,a,y);
où "a" et "b" sont les coeff de la fonction de transfert. Dans l'idée ça me va mais je ne comprends pas quoi mettre pour a et b car sur matlab la fonction de transfert est avec des z^-1 alors que dans tous les cours que j'ai trouvé sur le net les filtres coupes bandes sont exprimé dans le domaine de Laplace ou Fourier... De plus, je ne vois pas comment choisir a et b afin que je coupe entre 40 et 60Hz :ouin:

3°) Pour le filtre de "butter" (que je n'ai pas vraiment compris dans le principe) la mise en place sous matlab semble comme "fir1" :
Code:

1 2 3 4 x=data(:,1); y=data(:,2); frequenceBasse=40; frequenceHaute=60; ordreFiltre=1; [z,p,k] = butter(1,[2*pi*frequenceBasse 2*pi*frequenceHaute],'stop');
le soucis est que je ne sais pas trop quoi faire à la fin de [z,p,k] :aie: il s'agit de poles et de gain mais comment reconstruire mon signal à partir de ceci
=> dans l'aide je ne comprends pas qu'es ce qu'ils font... ça me semble super compliqué par rapport aux autres filtre :roll: il y a un truc que je n'ai pas dû saisir 8O

Pourriez vous, s'il vous plait, m'expliquer les points que je n'ai pas saisi pour ces filtres ?

merci beaucoup :ccool:
28/09/2012, 12h50
phryte

Bonjour,
Il y a une quatrième solution qui est le REJECTEUR.
En fait c'est un coupe bande très sélectif.
28/09/2012, 13h46
Invité
Bonjour,
1. Le tracé du spectre est fait de a à z dans la documentation de la fonction fft (déjà dit dans une discussion à laquelle tu as participé :roll:)
2. Pour faire ceci il semble y avoir ~~deux~~ LA solutions fir1 + filter.
  Si tu te poses toujours la question pour les coefficients a revoie tes cours ;) (fir = Finite Impulse Response)
3. Idem : butter + filter
28/09/2012, 14h29
membreComplexe12
merci beaucoup d'avoir pris le temps de répondre !!!

Citation:

Envoyé par phryte

Bonjour,
Il y a une quatrième solution qui est le REJECTEUR.
En fait c'est un coupe bande très sélectif.

merci Phrytre. En fait, dans ma tête coupe-bande c'est équivalent à réjecteur... ce n'est pas le cas ?

Citation:

Envoyé par Winjerome

Le tracé du spectre est fait de a à z dans la documentation de la fonction fft (déjà dit dans une discussion à laquelle tu as participé :roll:)

oui j'ai vu mais c'est pas le tracé qui m'intéresse mais comprendre ce qui est fait. En fait ni moi ni l'auteur de la discussion où j'ai participé nous n'avons compris pourquoi on définie comme ceci le vecteur de fréquence...

http://www.developpez.net/forums/d12...ur-frequences/

Citation:

Envoyé par Winjerome

Pour faire ceci il semble y avoir ~~deux~~ LA solutions fir1 + filter.
Si tu te poses toujours la question pour les coefficients a revoie tes cours ;) (fir = Finite Impulse Response)

=> fir1 me donne les coefficient de ma fonction de transfert
et ensuite il me suffit de les utiliser dans filter ? :ccool:
=> par contre dans "fir1" il y a un seul argument de sortie et dans "filter" il y a deux argument de sortie... du coup ça veut dire que les coefficient "a" de filter doivent tous être égaux à 1 pour ce type de filtre ?
(je n'ai jamais eu de cours sur le filtrage)

Citation:

Envoyé par Winjerome

Idem : butter + filter

Pour Butter j'ai compris du coup, merci.

je fais :
Code:

1 2 3 [b,a]=butter(n,Wn,'stop') yFiltre=filter(b,a,y)
et ça doit donner quelque chose de OK :zoubi:

ps: une dernière question : plus un filtre à un ordre élevé plus il filtre precisement une gamme de frequence ?

en fait il ne me reste qu'un petit soucis (car si j'ai bien compris a=1 pour FIR1)

=> comment définir le vecteur de fréquence X pour la fft ?
=> je sais qu'il y a un lien avec la fréquence d'echantillonage mais je ne vois pas trop lequel...

la frequence d'echantillonage va me donner une frequence discrete qui est le minimum que je peux avoir dans mon signal discret. Ensuite, pour tout les autres points se seront des multiples de cet fréquence ? et au final lorsque j'aurais fait "Fe*nbPoint" j'aurais le fondamental de mon signal ?

=> par contre, pourquoi dans l'exemple de matlab il y a une sorte de normalisation dans l'utilisation de la FFT ? (une division par L)

bref, pourriez vous m'éclaircir ce probleme d'echantillonage / frequence ?

merci d'avance

En fait mon gros probleme ce trouve dans ces lignes sur l'exemple de l'aide de matlab :
Code:

1 2 3 NFFT = 2^nextpow2(L); % Next power of 2 from length of y Y = fft(y,NFFT)/L; f = Fs/2*linspace(0,1,NFFT/2+1);
mais je ne comprends pas trop bien ceci aussi :
Code:

1 2 3 4 5 Fs = 1000; % Sampling frequency T = 1/Fs; % Sample time L = 1000; % Length of signal t = (0:L-1)*T; % Time vector plot(Fs*t(1:50),y(1:50))
les 4 première lignes ça va mais la dernière non :(
=> en fait dans les 4premieres lignes voici ce que l'on fait :
1°) définition frequence echantillonage
2°) definition periode echantillonage
3°) nombre de points que l'on veut
4°) on fait le vecteur temps qui correspond au nb de point à cette frequence
5°) ici "Fs*t(1:50)" je ne comprends pas ce que l'on fait...? si on aurait voulu tracer y=f(t) on aurait dû faire :

Code:

plot(t,y)

tout simplement ... ? pourquoi prendre que les 50premieres valeurs ? et pourquoi multiplier par la frequence Fs ???

30/09/2012, 20h04
thomMonteillet

ce genre de problème correspond exactement à ce que je ne comprends pas moi aussi.
En tout cas j'ai pu grâce à votre discussion sur ces différentes fonctions résoudre mon problème de filtrage.
Merci et bon courage pour la fft que je n'ai pas comprise moi aussi
30/09/2012, 20h32
Invité

Un lien qui sera plus parlant que quelques mots...

Citation:

Envoyé par 21did21

(je n'ai jamais eu de cours sur le filtrage)

Les cours sur le traitement du signal (et de ressources sur la FFT) sont plus que nombreux sur la toile, Google saura te conseiller ;)
30/09/2012, 22h13
membreComplexe12
merci beaucoup Jerôme :D, je vais regarder ceci. En effet, on dirait que c ce que je cherche :ccool:
Concernant Butter :
j'ai bien saisi à présent, voici le type de code que j'utilise :
Code:

1 2 3 [b,a]=butter(5,[0.7 0.9],'stop'); yfilter=filter(b,a,y);
par contre il y a un truc qui m'échappe : c'est quoi ces histoire de normalisation ?
autrement dit, si je veux filtrer un signal entre 200 et 300 Hz que doit je rentrer à la place de mon 0.7 et 0.9 ? (je normalise par rapport à quoi ?)
j'ai cru comprendre qu'il faillait faire un truc dans ce genre :

Code:

[b,a]=butter(5,[fbas/fEchan/2 fhaut/fEchan/2],'stop');

avec :
fbas : frequence basse de ma bande de coupure
fhaut : frequence haute de ma bande de coupure
fEchan : frequence d'echantillonage

=> par contre juste pour vérif : la frequence d'echantillonage et la fréquence d'espacement entre mes points du vecteur initial ? donc ça doit etre

Code:

1/(x(2)-x(1))

?
=> où la fréquence d'echantillonage n'est pas quelque chose donnée par notre signal mais quelque chose que l'on impose ?

(mais à ce moment là comment la fonction reconnait que c'est l'echantillonage qui à changé et non les fréquences que je souhaite couper ?)

j'ai regardé dans la FAQ mais j'ai l'impression qu'un truc m'echappe...
Concernant la FFT :

Code:

NFFT = 2^nextpow2(L);

cette ligne etait vraiment bloquante pour moi mais en lisant le lien que tu m'as envoyé j'ai cru comprendre que c'était une façon de faire qui permet d'optimiser les temps de calcul ?

je peux faire comme ceci à la place pour definir mon vecteur frequences ?

Code:

f=0 : Fech/Nfft : Fech/2

et la normalisation c'est quoi son intérêt ?
=> on peut bien remplacer :

Code:

Y = fft(y,NFFT)/L;

par :

Code:

Y = fft(y,NFFT);

et

Code:

f = Fs/2*linspace(0,1,NFFT/2+1);

par

Code:

f=0 : Fech/Nfft : Fech/2

je ne vois pas pourquoi ça ne marcherai pas...

je viens de faire un bout de code qui est la somme de deux sinus de fréquence différente.

ensuite, j'essai de filtrer un des deux sinus mais je n'y arrive pas :
Code:

1 2 3 4 5 6 7 8 9 10 clear all;close all;clc; %% signal origine t = 0:0.01:10; y = sin(2*pi*0.3*t)+sin(2*pi*2*t); %% signal traité fs = 1000; Wn = 2*[1 3]/fs; [b]=fir1(1,Wn,'stop'); yfilter=filter(b,1,y); figure (1);hold on;plot(t,y);plot(t,yfilter,'r');