Partager vos
astuces Déco
Posez votre question Signaler

Probleme matlab [Résolu]

falla 7Messages postés 4 mars 2008Date d'inscription - Dernière réponse le 7 janv. 2010 à 15:08
Bonjour,
j'ai un probleme de programmation sur matlab qui est un outil nouveau pour moi. il se trouve que j'ai un signal sinus et je doit faire l'extraction des harmoniques. j'ai réussi à faire la fft et de visualiser le fondamental et les harmoniques. donc il me reste à faire le filtrage de ces harmoniques(filtre passe bas).mon probleme c que je c pas comment programmer ca sur mathlab. apres ce travail je doit faire programme de la methode recurssive de la fft pour le comparer avec la fft. merci de votre comprehension je suis vraiment dans la merde c une petite partie de mon sujet de stage.
Lire la suite 

Probleme matlab »

Suggestions
Plus
29 réponses
Réponse
+4
moins plus
Sacabouffe 9604Messages postés 19 août 2007Date d'inscription 6 mars 2008 à 14:42
Bon... on va arrêter le carnage ici parce que c'est édifiant...

1) Un filtre passe-bas est un filtre qui laisse passer les basses fréquences et qui atténue les hautes fréquences, c'est-à-dire les fréquences supérieures à la fréquence de coupure. (Wikipédia)

2) Dans la mesure où un bruit est généralement un signal haute fréquence, filtrer un signal bruité en ne gardant que les hautes fréquences revient à garder seulement le bruit... d'une utilité phénoménale!!!

3) La fft Matlab est une transformée de Fourier discrète "standard" comme beaucoup d'autres (librairie fftw Fortran, etc...), par conséquent ça donne un spectre qui est shifté. Pour avoir correctement le spectre il faut utiliser la fonction fftshift de Matlab (ou faire à la main... mais bon... vu qu'il y a une fonction, autant en profiter...)

4) falla, ce n'est qu'un avis personnel mais je pense qu'un bruit blanc (rand) où un bruit gaussien (randn) modéliseraient mieux le bruit que l'on a sur un signal (si c'est bien ça que tu voulais faire)... 
Fs = 819050;
L = 16384;
T = L/(L-1)*1/Fs;
t = (0:L-1)*T;
f = linspace(-Fs/2,Fs/2,L);
fc = 100; %fréquence de coupure
freqplot = 1050; % fréquence max pour la courbe des harmoniques
Nmaxplot = find(f>=freqplot,1);

x = 55.2*sin(2*pi*50*t-2*pi*6/360)+...
42.78*sin(2*3*pi*50*t+2*pi*158/360)+23.69*sin(2*5*pi*50*t-2*pi*41/360)+...
7.82*sin(2*7*pi*50*t+2*pi*105/360)+3.68*sin(2*9*pi*50*t+ 2*pi*175/360)+...
3.91*sin(2*11*pi*50*t-2*pi*58/360)+1.84*sin(2*13*pi*50*t+ 2*pi*69/360)+...
1.38*sin(2*15*pi*50*t+2*pi*146/360)+1.15*sin(2*17*pi*50*t-2*pi*91/360)+...
0.69*sin(2*19*pi*50*t)+0.92*sin(2*21*pi*50*t+2*pi*104/360);

figure(1);
plot(t,x);
title('Signal Corrupted with Zero-Mean Random Noise');
xlabel('temps (s)');

y=1/Fs*fft(x);
y=fftshift(y);
figure(2);
%Plot single-sided amplitude spectrum
plot(f(L/2:Nmaxplot),abs(y(L/2:Nmaxplot)));
title('Amplitude Spectrum of y(t)');
xlabel('Frequency (Hz)');
ylabel('|Y(f)|') ;

%filtrage des harmoniques
Yfilter = linspace(0,Fs,L);
Yfilter = (abs(Yfilter-Fs/2)<=fc);
Yfilter = Yfilter.*y;
Yfilter = fftshift(Yfilter);
sig_filter=Fs*ifft(Yfilter);
figure(3);
plot(t,real(sig_filter))
title('Filtered signal');
xlabel('Temps (s)');

Dis-moi s'il y a encore des bugs...
Ciao

Ajouter un commentaire - Site web
falla - 6 mars 2008 à 17:48
ca marche mais je ne comprend pas ta fft car je ne vois pas du tout le fondamental (50hz) et les harmoniques. sinon ca marche a mon avis je vais le voir encors de prés ce soir et je vais vous répondre.
merci tu me redonne du courage.
Ajouter un commentaire
Réponse
+0
moins plus
Char Snipeur 4885Messages postés 23 avril 2004Date d'inscription 9 février 2012Dernière intervention 5 mars 2008 à 22:35
Voila un code qui fonctionne à peu près.
Mais décidemment, je ne vois pas pourquoi tu divise par L ta fft; c'est pour ça qu'a la fin pour obtenir un signal d'amplitude similaire à l'original je multiplie par L.
Le résultat ressemble à quelque chose, mais je ne suis pas sur que ça soit juste, pour le filtrage, il conviens peut être aussi de filtré les spectre "miroir" (ce qui est supérieur à NFFT/2)
Fs = 819050; % frequence d'echantillonnage
T = 1/Fs; % periode d'echantillonnage
L = 16384; % longueur du signal
t = (0:L-1)*T; % temps

x = 55.2*sin(2*pi*50*t-2*pi*6/360)+ 42.78*sin(2*3*pi*50*t-2*pi*158/360)+23.69*sin(2*5*pi*50*t-2*pi*41/360)+7.82*sin(-2*7*pi*50*t+ 2*pi*105/360)+3.68*sin(2*9*pi*50*t+ 2*pi*175/360)+3.91*sin(2*11*pi*50*t- 2*pi*58/360)+1.84*sin(2*13*pi*50*t+ 2*pi*69/360)+1.38*sin(2*15*pi*50*t+ 2*pi*146/360)+1.15*sin(2*17*pi*50*t- 2*pi*91/360)+0.69*sin(2*19*pi*50*t)+0.92*sin(2*21*pi*50*t+ 2*pi*104/360);
%x=sin(800000*t)+sin(2000*t);
x=x/(max(x)-min(x));
x=x-mean(x);

y = x;
plot(Fs*t(1:16300),y(1:16300))
title('Signal Corrupted with Zero-Mean Random Noise')
xlabel('temps (millisecondes)')
NFFT = 2^nextpow2(L); %nombre d'echatillon
Y = fft(y,NFFT)/L;
f = Fs/2*linspace(0,1,NFFT/2);

%Plot single-sided amplitude spectrum.
plot(f,2*abs(Y(1:NFFT/2)))
title('Single-Sided Amplitude Spectrum of y(t)')
xlabel('Frequency (Hz)')
ylabel('|Y(f)|')
pause

%filtrage des harmoniques
Yfilter = Y;
m = floor(size(Y));
%Yfilter(3)=Y(3);
[M_Y,im]=max(abs(Y));
Yfilter(1:2*im)=zeros(2*im,1);
sig_filter=ifft(Yfilter);
plot(f,2*abs(sig_filter(1:NFFT/2)))
title('Single-Sided1 Amplitude1 Spectrum1 of Yfilter')
xlabel('Temps en millisecondes')
pause

%reconversion en fréquentiel
Y1 = fft(sig_filter,NFFT/2)/L;
f = Fs/2*linspace(0,1,NFFT/2);
plot(f,2*abs(Y1(1:NFFT/2)))
title('Single-Sided2 Amplitude2 Spectrum2 of Yfilter')
xlabel('Frequency2 (Hz)')
pause

%comparaison du signal original et du signal filtré (filtre passe haut)
plot(t,x,t,abs(sig_filter)*L)
title('Superposition du signal et du signal filtré')
xlabel('temps [s]')

Ajouter un commentaire - Site web
falla - 6 mars 2008 à 14:35
la si tu veux ton programme tourne bien mais dans mes objectifs le signal filter doit etre un sinus et doit avoir une amplitude de 55.2.la j'essaye de trouver ca ce qui n'est pas facile.si tu a des idées par rapport à ca ?
falla - 6 mars 2008 à 17:36
dit moi si par hazrad tu n'a pas l'algo du filtre de butterworth passe bas peut etre ca va régler le probleme.
Ajouter un commentaire
Réponse
+0
moins plus
Char Snipeur 4885Messages postés 23 avril 2004Date d'inscription 9 février 2012Dernière intervention 6 mars 2008 à 15:35
salut.
1) En effet, erreur de ma part dans le message 2. dsl.
2) pour la fréquence du bruit, je ne suis pas d'accord. Si tu as un bruit "blanc", par définition, il est à toute les fréquences. Je pense que le mieux est de savoir dans quelle gamme est le signal et de faire un filtre passe bande (mais il est vrai que les basses fréquences sont moins gênantes).
3) pour ma culture, c'est quoi une fftshifté ?
Si tu ne veux conservé que ton "fondamentale", essai ma petite manip : 
[M_Y,im]=max(abs(Y));
Yfilter(im)=Y(im);
peut être te donnera elle un meilleur résultat.
Et je réitère (3ème fois) pourquoi divise tu par L ? Ou alors remultiplie à la fin.
Il est vrai que l'amplitude de la fft semble dépendre du nombre de point, mais c'est aussi valable pour l'inverse. Donc si tu diminue l'amplitude du spectre, tu diminuera l'amplitude du signal en fesant la transformé inverse.
Voilà ce que je peux te dire.

Ajouter un commentaire - Site web
Sacabouffe - 6 mars 2008 à 15:58
Salut Char Snipeur

1)Arf... oui... ce coup-ci c'est moi qui me suis planté, t'as tout à fait raison... :-D
Toutes mes excuses...

2) Oh ben là honnêtement je m'y connais pas trop mais oui... ça me paraît logique...

3) shiftée, c'est parce que la fréquence nulle est pas au centre du spectre. Si tu refais tous les calculs en utilisant la formule de la fft discrète tu verras qu'il y a des petits trucs à changer, tu peux pas l'utiliser comme ça directement pour les fonctions (sauf exception). La fonction fftshift de Matlab remet tout correctement

4) Pour les multiplications, divisions, c'est pareil... faut refaire les calculs, c'est un peu chiant...

Pour ton petit truc pour garder le fondamental, ça marche quand le fondamental a la plus grande amplitude, non?
Sinon ça te donne plutôt le mode dominant...
Mais pour falla je pense que ça colle du coup.

A plus
Sacabouffe - 6 mars 2008 à 16:45
Finalement je pense que sais pourquoi je pensais ça pour le bruit blanc...
Quand on rajoute un bruit blanc sur des données, on le fait avec un truc du type données=données+bruit_blanc_normalisé*données, non?
Dans ce cas, ça se manifeste le plus dans les hautes fréquences... je me plante encore?
falla - 6 mars 2008 à 17:34
parceque si je divise par L ma fft aura de bonne amplitude pour pour chaque harmonique.par contre si je ne sivise pas par L le fondamental aurra une amplitude supérieur à 1300.ce qui n'est pas normale car l'amplitude du signal pour le fondamental est de 55.2.sinon je sais que c dure ce truc mais je vais toujours continuer jusqu'a le réussir car sinon je vais me retrouver a la porte. sinon j'ai essayé ton example mais j'ai toujours mon soucis qui est mon signal sinus qui correspond au 50Hz.
merci encors pour l'aide
Ajouter un commentaire
Réponse
+0
moins plus
Char Snipeur 4885Messages postés 23 avril 2004Date d'inscription 9 février 2012Dernière intervention 7 mars 2008 à 08:55
Ok, je comprends ton raisonnement, mais je doute qu'il soit juste.
Comme le fesait remarquer Sacabouffe, la transformée de Fourier d'un sinus est une fonction de Dirac, c'est à dire qu'elle tend vers l'infini pour la fréquence du sinus, et ce quel que soit son amplitude. Je ne suis pas sur qu'il faille chercher un rapport direct en amplitude du spectre et amplitude de l'onde.
Toi tu fait un signalà partir de fonction réellement périodique, alors que sous MATLAB tu ne donne qu'un échantillon.
Si tu fait la transformée de Fourier de sin(w.t) F=TF(sin(w.t)); alors F(w)->+oo
Par contre, si tu as un signal qui vaut zéro sauf entre 0 et T (w=2*pi/T)
Alors F(w)~0.5, si le signal dure entre 0 et 2T, F(w)~1
Tu obtiens une amplitude de spectre fonction du nombre de périodes effectuées sans rapport avec l'amplitude devant le sinus qui est ici de 1 !
à mon avis, si tu veux remoner à l'amplitude de chaque sinusoïde , il vaut mieux diviser par le nombre de cycle effectué durant ton temps de mesure.
Enfin, après je peux me tromper, je ne suis pas expert en traitement du signal.

Quant au bruit blanc... je ne sais pas trop. Entre comment on le défini et comment il est dans la réalité...
En cour, il me semble que nous ajoutions simplement le bruit, sans le multiplier. Après dans ta façon de faire, pour moi ce n'est pas évident que ça intervienne plus dans les hautes fréquences.

Ajouter un commentaire - Site web
Afficher les 9 commentaires
Sacabouffe - 15 mars 2008 à 18:15
Salut
Désolé j'y comprends que dalle.
C'est du Fortran ça?
Ça ressemble à du Basic des vieux CPC... :-D
Désolé... à plus.
thamina84 - 23 nov. 2009 à 10:10
bjr sacabouffe,

g besoin d aide sur la résolution des système masse ressort amorti (second ordre) et la visualisation de la réponse sous matlab
(par exemple oscillateur de duffing )ou un exemple simple
g désespérément besoin d'aide parce que ce ne sont que des exemples simples par rapport à ce que j'ai à faire :( je travaille sur la commande d'un système masse ressort à un seul degré de liberté que je dois visualiser la réponse temporelle )
merci :D
thamina84 › Sacabouffe - 7 janv. 2010 à 15:08
bonjour,

j'ai un programme matlab sur lequel je travaille et je troue des difficultés.
j'ai un systéme discret qui s'écrit sous la forme:

x[n+1] = Ax[n] + Bu[n]
y[n] = Cx[n] + Du[n]
u étant un vecteur de perturbation sinusoïdale
dt=0.05;
Tmp=[0:dt:20];
n_pas=length(Tmp);
u=0.5*sin(2*pi*F*(Tmp+dt))
j'arrive à représenter la réponse temporelle par
[y,x]=dlsim(A,P,C,D,u(:,1:n_pas));

mais je veux représenter la réponse fréquentielle si je varie la fréquence de perturbation F de u
comment faire?
j'ai vraiment besoin d'aide

meeerci
Ajouter un commentaire
Réponse
-1
moins plus
Sacabouffe 9604Messages postés 19 août 2007Date d'inscription 4 mars 2008 à 16:36
Salut falla
1) Un signal sinusoïdal n'a pas d'harmonique
2) La TF d'un sinus, c'est un Dirac, il y a rien à filtrer...
3) Si donc... ton signal n'est pas un signal sinusoïdal... mais simplement périodique...
Quel type de filtre passe-bas tu veux? Quelle fréquence de coupure, quelle décroissance en fonction de la fréquence, etc, etc?
A plus tard

Ajouter un commentaire - Site web
Ajouter un commentaire
Réponse
-1
moins plus
Char Snipeur 4885Messages postés 23 avril 2004Date d'inscription 9 février 2012Dernière intervention 4 mars 2008 à 16:44
Salut.
un filtre passe bas consiste à supprimer toutes les fréquences basses.
donc un truc dans le genre :
spectre=fft(signal);
i=spectre.length();
spectre_filtré=[zeros(1,i/2) , spectre(i/2:end) ];
signal_filtré=ifft(spectre_filtré);
plot(signal_filtré)

Ajouter un commentaire - Site web
falla - 5 mars 2008 à 14:56
bonjour je suis tres content de votre reponse car ca ma donné du jus.sinon voici mon signal et la transformée de fourier ,le filtrage puis l'inverse et la conversion en fréquentiel pour voir si j'ai seulement mon fondamental. mais la il se trouve que apres filtrage je constate que l'amplitude du fondamental est trés faible donc je me dit qu'il a tout filtré donc je me dit qu'il y a quelque chose qui bloque.

Fs = 819050; % frequence d'echantillonnage
T = 1/Fs; % periode d'echantillonnage
L = 16380; % logueur du signal
t = (0:L-1)*T; % temps

x = 55.2*sin(2*pi*50*t-2*pi*6/360)+ 42.78*sin(2*3*pi*50*t-2*pi*158/360)+23.69*sin(2*5*pi*50*t-2*pi*41/360)+7.82*sin(2*7*pi*50*t+ 2*pi*105/360)+3.68*sin(2*9*pi*50*t+ 2*pi*175/360)+3.91*sin(2*11*pi*50*t- 2*pi*58/360)+1.84*sin(2*13*pi*50*t+ 2*pi*69/360)+1.38*sin(2*15*pi*50*t+ 2*pi*146/360)+1.15*sin(2*17*pi*50*t- 2*pi*91/360)+0.69*sin(2*19*pi*50*t)+0.92*sin(2*21*pi*50*t+ 2*pi*104/360)

y = x;
plot(Fs*t(1:16300),y(1:16300))
title('Signal Corrupted with Zero-Mean Random Noise')
xlabel('temps (millisecondes)')
NFFT = 2^nextpow2(L); %nombre d'echatillon
Y = fft(y,NFFT)/L;
f = Fs/2*linspace(0,1,NFFT/2);

%Plot single-sided amplitude spectrum.
plot(f,2*abs(Y(1:NFFT/2)))
title('Single-Sided Amplitude Spectrum of y(t)')
xlabel('Frequency (Hz)')
ylabel('|Y(f)|')

%filtrage des harmoniques
Yfilter = 0*Y;
m = floor(size(Y));
Yfilter(3)=Y(3);
sig_filter=ifft(Yfilter);
plot(f,2*abs(sig_filter(1:NFFT/2)))
title('Single-Sided1 Amplitude1 Spectrum1 of Yfilter')
xlabel('Temps en millisecondes')

%reconversion en fréquentiel
Y1 = fft(sig_filter,NFFT/2)/L;
f = Fs/2*linspace(0,1,NFFT/2);
plot(f,2*abs(Y1(1:NFFT/2)))
title('Single-Sided2 Amplitude2 Spectrum2 of Yfilter')
xlabel('Frequency2 (Hz)')
Ajouter un commentaire
Réponse
-1
moins plus
Char Snipeur 4885Messages postés 23 avril 2004Date d'inscription 9 février 2012Dernière intervention 5 mars 2008 à 21:53
Salut.
J'ai récupérer ton code, j'ai un peu de mal à comprendre, il faudrait que je me remette aux fft.
Je ne sais pas exactement pourquoi ton signal final est quasi nul, mais il faut dire que ton filtrage est sacrément violant !
Yfilter est un vecteur nul avec une seule valeur non nulle...
je pense que ton filtrage n'est pas bon.

Ajouter un commentaire - Site web
falla - 6 mars 2008 à 08:51
pour le signal j'ai fait une erreur de signe que j'ai rectifier donc la j'obtient bien mon signal pollué mais j'ai toujours le meme probléme.
merci pour votre compréhension.

x = 55.2*sin(2*pi*50*t-2*pi*6/360)+ 42.78*sin(2*3*pi*50*t+2*pi*158/360)+23.69*sin(2*5*pi*50*t-2*pi*41/360)+7.82*sin(2*7*pi*50*­t+ 2*pi*105/360)+3.68*sin(2*9*pi*50*t+ 2*pi*175/360)+3.91*sin(2*11*pi*50*t- 2*pi*58/360)+1.84*sin(2*13*pi*50*t+ 2*pi*69/360)+1.38*sin(2*15*pi*50*t+ 2*pi*146/360)+1.15*sin(2*17*pi*50*t- 2*pi*91/360)+0.69*sin(2*19*pi*50*t)+0.92*sin(2*21*pi*50*t+ 2*pi*104/360)
Char Snipeur - 6 mars 2008 à 10:26
Tu as essayé ce que je t'ai envoyé ?
Je pense vraiment que ton filtrage est trop violent. Ou alors, multiplie ton Yfilter par L, voire par plus.
Il faut garder à l'esprit qu'une fft n'est pas parfaite, car tout signal fini dans le temps comporte une incertitude sur la fréquence (c'est d'ailleurs une bonne façon d'illustrer l'incertitude d'Heisenberg au passage) Donc, tant que ton signal n'aura pas une durée infini tu n'aurai pas un vrai spectre discret, mais il y aura une dispersion : filtre moins fort !
c'est sur que tes valeurs sont faible, come je te l'ai dit, Y(3) est faible, tu divise par L et tu ne prend qu'une valeur pour ta transformée inverse.
Ajouter un commentaire
Réponse
-1
moins plus
hahad 6 juin 2008 à 13:56
je suis en stage et j'ai un problemme au niveau de matlab
en faite, je dois faire un programme sur matlab qui calacule la transformé de fourier inverse d'une fonction qlcq sauf qui doit la faire sur une bande comprise entre 40Mhz et 18 GHz .

merci de me repondre je suis vraiment bloqué .

 Ajouter un commentaire
Sacabouffe - 6 juin 2008 à 14:21
Salut
ifft et tu choisis l'intervalle
A plus
Ajouter un commentaire
Réponse
-2
moins plus
maryam 8 mars 2008 à 10:08
bonjour
si vous permettez je rencontre n grand probléme avec matlab, je ne suis qu'une débutante , mais à la moindre des choses il ne m'affiche pas le résultat d'un histogramme
c'est quoi le truc???

 Ajouter un commentaire
Sacabouffe - 9 mars 2008 à 01:37
Salut
Fais ton propre sujet et explique mieux les choses.
A plus
Ajouter un commentaire
Posez votre question
Ce document intitulé « probleme matlab » issu de CommentCaMarche (www.commentcamarche.net) est mis à disposition sous les termes de la licence Creative Commons. Vous pouvez copier, modifier des copies de cette page, dans les conditions fixées par la licence, tant que cette note apparaît clairement.
Dossier à la une
5 extensions si vous voulez revenir à l'ancien Facebook