Bonjour,
Je cherche à retrouver temporellement la réponse impulsionnelle d'un système en connaissant les accélérations d'entrée et de sortie de mon système. Je cherche donc à passer par Fourier avec la fonction fft de Matlab, diviser la transformée de ma sortie par celle de mon entrée, retourner en temporel grace à la transformée inverse ifft et donc en théorie avoir ma réponse. Or quand je teste sur des cas simples je n'obtiens pas les bons résultats, je ne sais pas quelle commande j'oublie.
mon code :
Code:function [hf] = fftinv(t,ae,as,fsyst) %on cherche à déconvoluer un signal de sorte à retrouver la rép %impulsionnelle d'un systeme %t est le vecteur temps %ae est l'accélération en entrée du systeme %as est l'accélération mesurée en sortie %fsyst est la fréquence d'excitation du systeme dt=t(3)-t(2); fech=1/dt; %paramètres physiques du systeme m=10; %masse ksi=0.05; %taux d'amortissement voulu omega=2*pi*fsyst; k=m*omega*omega; c=2*ksi*m*omega; %amortissement omegaD=omega*sqrt(1-ksi*ksi); %lorsque je cherche à obtenir la sortie de mon systeme en connaissant la réponse impulsionnelle et l'entrée, je devais multiplier le résultat de ma convolution pour obtenir ma vraie accélération (Duhamel), donc je divise par la meme chose de sorte à ce que as corresponde au résultat de mon produit de convolution. as=m*as*omegaD/(omega*omega*dt); %transfo de fourier de ae et as Nt=length(ae); Nfft=2^nextpow2(Nt); freq=fech/2*linspace(0,1,Nfft/2); %calcul du spectre aef=fft(ae,Nfft); asf=fft(as,Nfft); figure; plot(freq, abs(aef(1:Nfft/2,1))); hold on plot(freq, abs(asf(1:Nfft/2,1)),'color','red'); hold on; xlabel('frequence (Hz)'); % jusque là, lorsque j'entre un sinus tq f0=5Hz, j'ai bien un seul pic à 5Hz donc ok hf=asf./aef; figure; plot( freq,abs(hf(1:Nfft/2,1))); % dans ce spectre, un pic à 20Hz apparait et je ne sais pas d'ou il vient .. h=real(ifft(hf)); h=h(1:Nt); %on ne prend que les Nt premiers points, qui ont la vraie info temporelle figure; plot(t,h) %là aussi il y a probleme
J'entre t de 0 à 4 avec dt=0.004, ae=10*sin(2*pi*5*t) jusqu'à 2sec puis 0 ensuite et ma réponse impulsionnelle h est définie par h(i+1)=exp(-ksi*omega*t(i+1))*(sin(omegaD* t(i+1)). Je calcule as par convolution entre ae et h et je multiplie par dt*w²/(m*wd). Je cherche donc à retrouver h en utilisant le programme ci-dessus. Le résultat est juste jusqu'à 2sec, puis devient n'importe quoi ensuite (il reprend comme au début au lieu de continuer de décroitre vers 0).
Je cherche ensuite à tester sur de vrais signaux sismiques et là j'obtiens une "bande" très proche de 0 donc également fausse.
Je suppose qu'il s'agit d'une ligne manquante/fausse mais je ne vois pas où. Si qqn a des idées, merci d'avance!
-----