| 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.
Salutation ! avant je croyais, maintenant je suis fixé.Jésus Christ
Char Snipeur | 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
Some folks are born made to wave the flag, ooo, they’re red, white and blue.
And when the band plays "Hail to the Chief", ooo, they point the cannon at you, y’all! | 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
Some folks are born made to wave the flag, ooo, they’re red, white and blue.
And when the band plays "Hail to the Chief", ooo, they point the cannon at you, y’all! | 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
Some folks are born made to wave the flag, ooo, they’re red, white and blue.
And when the band plays "Hail to the Chief", ooo, they point the cannon at you, y’all! | 22 falla, le 11 mar 2008 à 13:54:33Merci 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 | Salut
Je regarde ça...
A plus
Some folks are born made to wave the flag, ooo, they’re red, white and blue.
And when the band plays "Hail to the Chief", ooo, they point the cannon at you, y’all! |
|
|
|
|
| 23 falla, le 12 mar 2008 à 12:57:04Il 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 | Salut
Désolé j'y comprends que dalle.
C'est du Fortran ça?
Ça ressemble à du Basic des vieux CPC... :-D
Désolé... à plus.
Some folks are born made to wave the flag, ooo, they’re red, white and blue.
And when the band plays "Hail to the Chief", ooo, they point the cannon at you, y’all! |
|
|
Communication rs232 d'un gps et un pc avec matlab un code
