Probleme matlab

Résolu/Fermé
falla Messages postés 7 Date d'inscription mardi 4 mars 2008 Statut Membre Dernière intervention 6 mars 2008 - 4 mars 2008 à 14:07
 (abdou) - 11 sept. 2012 à 23:07
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.

9 réponses

Sacabouffe Messages postés 9427 Date d'inscription dimanche 19 août 2007 Statut Membre Dernière intervention 29 mai 2009 1 832
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
7
falla Messages postés 7 Date d'inscription mardi 4 mars 2008 Statut Membre Dernière intervention 6 mars 2008
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.
0
bj,
j'ai une petite problème avec les deux fonction ifft et fft. j'ai appliqué l' ifft a l'émission avec Zéros pading et au niveau de la réception j'ai appliqué l' fft mais je récupère une composante imaginaire 0.5954+0.0000i avec tous les bits ce qui cause des erreurs dans mon simulation:
function d=IFFT_ZP(X)
l1=length(X);
I=ifft(X);
l=length(I);
for n=128:128:l
if mod(n,128)==0
I=[I zeros(1,37)];
end
end
d=I;
end
***********************************
function s1=OLA_FFT(d)
m=length(d);
s=[];
n=1;
for u=1:m
if d(u)~=0
s(n)=d(u);
n=n+1;
end
end
s1=fft(s);
s1=real(s1);
end

Est ce que je peut obtenir une idée pour ce problème.
0
salut,Sacabouffe j'ai vue votre progrmme matlab et je suis intérréssé de la fft mais je travail dans matlab-simulink est ce que vous avez une idée ?
mercie.
(abdou)
0
Sacabouffe Messages postés 9427 Date d'inscription dimanche 19 août 2007 Statut Membre Dernière intervention 29 mai 2009 1 832
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
0
Char Snipeur Messages postés 9696 Date d'inscription vendredi 23 avril 2004 Statut Contributeur Dernière intervention 3 octobre 2023 1 297
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]')
0
falla Messages postés 7 Date d'inscription mardi 4 mars 2008 Statut Membre Dernière intervention 6 mars 2008
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 ?
0
falla Messages postés 7 Date d'inscription mardi 4 mars 2008 Statut Membre Dernière intervention 6 mars 2008
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.
0
Char Snipeur Messages postés 9696 Date d'inscription vendredi 23 avril 2004 Statut Contributeur Dernière intervention 3 octobre 2023 1 297
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.
0
Sacabouffe Messages postés 9427 Date d'inscription dimanche 19 août 2007 Statut Membre Dernière intervention 29 mai 2009 1 832
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
0
Sacabouffe Messages postés 9427 Date d'inscription dimanche 19 août 2007 Statut Membre Dernière intervention 29 mai 2009 1 832 > Sacabouffe Messages postés 9427 Date d'inscription dimanche 19 août 2007 Statut Membre Dernière intervention 29 mai 2009
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?
0
falla Messages postés 7 Date d'inscription mardi 4 mars 2008 Statut Membre Dernière intervention 6 mars 2008
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
0

Vous n’avez pas trouvé la réponse que vous recherchez ?

Posez votre question
Char Snipeur Messages postés 9696 Date d'inscription vendredi 23 avril 2004 Statut Contributeur Dernière intervention 3 octobre 2023 1 297
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.
0
Sacabouffe Messages postés 9427 Date d'inscription dimanche 19 août 2007 Statut Membre Dernière intervention 29 mai 2009 1 832
7 mars 2008 à 14:51
Ah... ben j'en sais rien... chaque fois que j'ai vu un bruit ajouté au données, c'était fait comme ça... Mais j'ai jamais vraiment réussi à savoir si physiquement c'était correct... Du coup j'utilisais le même :-D
Quant au raisonnement, ben... là non plus je sais pas trop... mais numériquement ça marchait super bien. En tronquant la TF on enlevait toutes les hautes fréquences et du coup même en mettant un bruit super fort (avec le modèle dont j'ai parlé) les résultats étaient quasiment toujours les mêmes.
Pour le code que j'ai écrit, je pense qu'il y a un petit bug. J'essaie de voir ça plus tard. Il y a un décalage de 25Hz dans le spectre j'ai l'impression.

falla, ce que tu vois avec le code que j'ai écrit c'est justement la fréquence 50Hz, si tu veux voir plus d'harmonique, modifie le fc. Mais comme je le disais, il y a encore un truc qui cloche.

A plus
0
Sacabouffe Messages postés 9427 Date d'inscription dimanche 19 août 2007 Statut Membre Dernière intervention 29 mai 2009 1 832 > Sacabouffe Messages postés 9427 Date d'inscription dimanche 19 août 2007 Statut Membre Dernière intervention 29 mai 2009
10 mars 2008 à 01:20
Bon j'ai corrigé 2 ou 3 trucs...
On oublie le fftshift, vu que le signal de départ est pas symétrique par rapport à 0, j'arrive pas à l'utiliser correctement et... ben... j'ai pas envie de chercher :-D
J'ai rajouté un truc parce qu'en fait, la ifft du spectre tronqué a une composante imaginaire qu'il faut prendre en compte, on peut pas l'enlever comme j'ai fait... Par contre ça marche seulement pour prendre la fondamentale (fréquence de coupure fc=50), faut bidouiller si on veut prendre plus d'harmoniques, toujours pour prendre en compte de manière correcte les composantes imaginaires... mais pour l'instant j'arrive pas à le faire de manière correcte... ça me gave!!!
Char Snipeur!!! Une idée s'il te plaît!!!
Fs = 819200;
L = 16384;
T = 1/Fs;
t = (0:L-1)*T;
f = linspace(0,(L-1)*Fs/L,L);
fc = 50; %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(1:Nmaxplot),abs(y(1:Nmaxplot)));
title('Amplitude Spectrum of y(t)');
xlabel('Frequency (Hz)');
ylabel('|Y(f)|') ;

%filtrage des harmoniques
Yfilter = linspace(0,(L-1)*Fs/L,L);
Yfilter = (Yfilter<=fc);
Yfilter = Yfilter.*y;
%Yfilter = fftshift(Yfilter);
sig_filter=Fs*ifft(Yfilter);
figure(3);
% la fondamentale c'est la composante réelle plus la composante imaginaire
% déphasée de pi/2
% en espérant que je raconte pas de conneries... je revérifierai les calculs...
nper = find(t==1/(4*50),1);
plot(t,real(sig_filter)+circshift(-imag(sig_filter),[0 4097])); %
title('Filtered signal');
xlabel('Temps (s)');

Bon voilà... En espérant que Char Snipeur aura une idée lumineuse... :-D
0
Sacabouffe Messages postés 9427 Date d'inscription dimanche 19 août 2007 Statut Membre Dernière intervention 29 mai 2009 1 832 > Sacabouffe Messages postés 9427 Date d'inscription dimanche 19 août 2007 Statut Membre Dernière intervention 29 mai 2009
10 mars 2008 à 16:55
Non, non, non et non!!! C'est tout pourri ce que j'ai écrit!!!
D'une manière ou d'une autre il faut l'utiliser le fftshift sinon il y a des harmoniques qui restent à la fin du spectre... il faut recentrer tout ça... pas le choix...

Le décalage de 25Hz que j'avais avant... ben c'est parce qu'il y avait un nombre pair de fréquence pour le spectre, donc quand on recentre avec fftshift et ben... c'est pas possible de recentrer en 0 puisque les fréquences du spectre sont un truc du genre ...-75 -25 25 75...
Et sinon... le fftshift sous Matlab est pas une application involutive (je viens de le découvrir...). Donc obligé de faire un circshift aussi. Maintenant je pense que c'est bon...
Fs = 819150;
L = 16383;
T = 1/Fs;
t = (0:L-1)*T;
f = linspace(-(L-1)*Fs/2/L,(L-1)*Fs/2/L,L);
nharm = 1; % nombre d'harmoniques voulues
fc = 50*(2*nharm-1); % 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+1)/2:Nmaxplot),abs(y((L+1)/2:Nmaxplot)));
title('Amplitude Spectrum of y(t)');
xlabel('Frequency (Hz)');
ylabel('|Y(f)|') ;

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

A plus
0
falla > Sacabouffe Messages postés 9427 Date d'inscription dimanche 19 août 2007 Statut Membre Dernière intervention 29 mai 2009
11 mars 2008 à 13:54
merci pour ta réponse je l'ai essayé mais je vais y faire quelque réglages. La il se trouve que sujet à un peu changé de forme car je doit totalement écrire l'algorithme de FFT ou de la FFT récursive pour l'adapté à réseau électrique qui débite sur une charge polluante. donc l'objectif est de traduire l'algoritme de la FFt mathématique en matlab qui me permet d'afficher tous le 20ms la FFT et de sélection les harmoniques à compenser afin de faire la IFFT de ces harmoniques et de les réinjecter au réseau. Ce programme doit pouvoir filtré jusqu'au rang 50.Donc pour le paramétre je vais prendre N=256, Fe=12800, le calcul se fait tous les 20 ms.
voici un programme en fortran qui calcul la fft mais bon je ne sait pas si je pourrai le traduire en matlab et puis l'améliorer pour qu'il puisse faire la sélection de mes harmonique à compenser.

1000 'THE FAST FOURIER TRANSFORM
1010 'Upon entry, N% contains the number of points in the DFT, REX[ ] and
1020 'IMX[ ] contain the real and imaginary parts of the input. Upon return,
1030 'REX[ ] and IMX[ ] contain the DFT output. All signals run from 0 to N%-1.
1040 '
1050 PI = 3.14159265 'Set constants
1060 NM1% = N%-1
1070 ND2% = N%/2
1080 M% = CINT(LOG(N%)/LOG(2))
1090 J% = ND2%
1100 '
1110 FOR I% = 1 TO N%-2 'Bit reversal sorting
1120 IF I% >= J% THEN GOTO 1190
1130 TR = REX[J%]
1140 TI = IMX[J%]
1150 REX[J%] = REX[I%]
1160 IMX[J%] = IMX[I%]
1170 REX[I%] = TR
1180 IMX[I%] = TI
1190 K% = ND2%
1200 IF K% > J% THEN GOTO 1240
1210 J% = J%-K%
1220 K% = K%/2
1230 GOTO 1200
1240 J% = J%+K%
1250 NEXT I%
1260 '
1270 FOR L% = 1 TO M% 'Loop for each stage
1280 LE% = CINT(2^L%)
1290 LE2% = LE%/2
1300 UR = 1
1310 UI = 0
1320 SR = COS(PI/LE2%) 'Calculate sine & cosine values
1330 SI = -SIN(PI/LE2%)
1340 FOR J% = 1 TO LE2% 'Loop for each sub DFT
1350 JM1% = J%-1
1360 FOR I% = JM1% TO NM1% STEP LE% 'Loop for each butterfly
1370 IP% = I%+LE2%
1380 TR = REX[IP%]*UR - IMX[IP%]*UI 'Butterfly calculation
1390 TI = REX[IP%]*UI + IMX[IP%]*UR
1400 REX[IP%] = REX[I%]-TR
1410 IMX[IP%] = IMX[I%]-TI
1420 REX[I%] = REX[I%]+TR
1430 IMX[I%] = IMX[I%]+TI
1440 NEXT I%
1450 TR = UR
1460 UR = TR*SR - UI*SI
1470 UI = TR*SI + UI*SR
1480 NEXT J%
1490 NEXT L%
1500 '
1510 RETURN

merci encore
0
il se trouve que sujet à un peu changé de forme car je doit totalement écrire l'algorithme de FFT ou de la FFT récursive pour l'adapté à réseau électrique qui débite sur une charge polluante. donc l'objectif est de traduire l'algoritme de la FFt mathématique en matlab qui me permet d'afficher tous le 20ms la FFT et de sélection les harmoniques à compenser afin de faire la IFFT de ces harmoniques et de les réinjecter au réseau. Ce programme doit pouvoir filtré jusqu'au rang 50.Donc pour le paramétre je vais prendre N=256, Fe=12800, le calcul se fait tous les 20 ms.
voici un programme en fortran qui calcul la fft mais bon je ne sait pas si je pourrai le traduire en matlab et puis l'améliorer pour qu'il puisse faire la sélection de mes harmonique à compenser.

1000 'THE FAST FOURIER TRANSFORM
1010 'Upon entry, N% contains the number of points in the DFT, REX[ ] and
1020 'IMX[ ] contain the real and imaginary parts of the input. Upon return,
1030 'REX[ ] and IMX[ ] contain the DFT output. All signals run from 0 to N%-1.
1040 '
1050 PI = 3.14159265 'Set constants
1060 NM1% = N%-1
1070 ND2% = N%/2
1080 M% = CINT(LOG(N%)/LOG(2))
1090 J% = ND2%
1100 '
1110 FOR I% = 1 TO N%-2 'Bit reversal sorting
1120 IF I% >= J% THEN GOTO 1190
1130 TR = REX[J%]
1140 TI = IMX[J%]
1150 REX[J%] = REX[I%]
1160 IMX[J%] = IMX[I%]
1170 REX[I%] = TR
1180 IMX[I%] = TI
1190 K% = ND2%
1200 IF K% > J% THEN GOTO 1240
1210 J% = J%-K%
1220 K% = K%/2
1230 GOTO 1200
1240 J% = J%+K%
1250 NEXT I%
1260 '
1270 FOR L% = 1 TO M% 'Loop for each stage
1280 LE% = CINT(2^L%)
1290 LE2% = LE%/2
1300 UR = 1
1310 UI = 0
1320 SR = COS(PI/LE2%) 'Calculate sine & cosine values
1330 SI = -SIN(PI/LE2%)
1340 FOR J% = 1 TO LE2% 'Loop for each sub DFT
1350 JM1% = J%-1
1360 FOR I% = JM1% TO NM1% STEP LE% 'Loop for each butterfly
1370 IP% = I%+LE2%
1380 TR = REX[IP%]*UR - IMX[IP%]*UI 'Butterfly calculation
1390 TI = REX[IP%]*UI + IMX[IP%]*UR
1400 REX[IP%] = REX[I%]-TR
1410 IMX[IP%] = IMX[I%]-TI
1420 REX[I%] = REX[I%]+TR
1430 IMX[I%] = IMX[I%]+TI
1440 NEXT I%
1450 TR = UR
1460 UR = TR*SR - UI*SI
1470 UI = TR*SI + UI*SR
1480 NEXT J%
1490 NEXT L%
1500 '
1510 RETURN
0
Char Snipeur Messages postés 9696 Date d'inscription vendredi 23 avril 2004 Statut Contributeur Dernière intervention 3 octobre 2023 1 297
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é)
-1
falla Messages postés 7 Date d'inscription mardi 4 mars 2008 Statut Membre Dernière intervention 6 mars 2008
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)')
0
Char Snipeur Messages postés 9696 Date d'inscription vendredi 23 avril 2004 Statut Contributeur Dernière intervention 3 octobre 2023 1 297
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.
-1
falla Messages postés 7 Date d'inscription mardi 4 mars 2008 Statut Membre Dernière intervention 6 mars 2008
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)
0
Char Snipeur Messages postés 9696 Date d'inscription vendredi 23 avril 2004 Statut Contributeur Dernière intervention 3 octobre 2023 1 297 > falla Messages postés 7 Date d'inscription mardi 4 mars 2008 Statut Membre Dernière intervention 6 mars 2008
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.
0
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é .
-1
Sacabouffe Messages postés 9427 Date d'inscription dimanche 19 août 2007 Statut Membre Dernière intervention 29 mai 2009 1 832
6 juin 2008 à 14:21
Salut
ifft et tu choisis l'intervalle
A plus
0
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???
-2
Sacabouffe Messages postés 9427 Date d'inscription dimanche 19 août 2007 Statut Membre Dernière intervention 29 mai 2009 1 832
9 mars 2008 à 01:37
Salut
Fais ton propre sujet et explique mieux les choses.
A plus
0