Répondre à la discussion
Affichage des résultats 1 à 4 sur 4

probleme fft/ifft matlab



  1. #1
    blabla1010

    probleme fft/ifft matlab


    ------

    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!

    -----
    Dernière modification par Antoane ; 04/07/2015 à 19h40. Motif: Ajout balises code

  2. Publicité
  3. 📣 Nouveau projet éditorial de Futura
    🔥🧠 Le Mag Futura est lancé, découvrez notre 1er magazine papier

    Une belle revue de plus de 200 pages et 4 dossiers scientifiques pour tout comprendre à la science qui fera le futur. Nous avons besoin de vous 🙏 pour nous aider à le lancer...

    👉 Je découvre le projet

    Quatre questions à explorer en 2022 :
    → Quels mystères nous cache encore la Lune 🌙 ?
    → Pourra-t-on bientôt tout guérir grâce aux gènes 👩‍⚕️?
    → Comment nourrir le monde sans le détruire 🌍 ?
    → L’intelligence artificielle peut-elle devenir vraiment intelligente 🤖 ?
  4. #2
    phuphus

    Re : probleme fft/ifft matlab

    Bonsoir,

    déjà, il n'est pas normal de ne pas avoir directement la bonne as à partir de ae et de la RI. Si c'est le cas, c'est que la RI n'est pas OK.

    Tu as un signal réel, dont la FFT est symétrique conjuguée par rapport à la fréquence de Nyquist. Avant de faire la FFt inverse de ta fonction de transfert en fréquentiel, il faut reconstruire tous les résultats au delà de la fréquence de Nyquist via cette symétrie conjuguée. Evite aussi le zero padding en te mettant sur la prochaine puissance de 2 : c'est source de problèmes si tu ne manipules pas bien ces concepts.

  5. #3
    phuphus

    Re : probleme fft/ifft matlab

    Bonjour,

    autre chose : si tu veux obtenir ta RI à partir des deux signaux ae et as, il faut t'assurer que tout point de as est dû aux données présentes dans ae. Classiquement, cela se fait avec un système au repos au début et à la fin de chaque signal. A défaut, tu pourras minimiser l'influence des parties de as n'étant pas dues à ae (le début), et réciproquement aux parties de ae n'influençant pas as (la fin) en prenant ae et as très longs devant la durée significative de la RI.

  6. #4
    blabla1010

    Re : probleme fft/ifft matlab

    Bonjour,
    Merci de votre réponse.

    Le fait que je multiplie le résultat de ma convolution de ae par ma RI est du à Duhamel je crois (pour un cas parfait, la réponse du système s'écrit sous la forme d'une intégrale multipliée par un coeff devant, c'est ce coeff que j'ai rajouté).

    Comment reconstruire les résultats via la symétrie conjuguée sous Matlab ? J'avais déjà essayé de faire hf=hf*conj(hf)/Nfft mais ça n'apportait rien.

  7. A voir en vidéo sur Futura

Discussions similaires

  1. Matlab IFFT
    Par maestro88 dans le forum Électronique
    Réponses: 4
    Dernier message: 02/07/2010, 17h43
  2. ifft sur Matlab PSD
    Par pamath dans le forum Mathématiques du supérieur
    Réponses: 0
    Dernier message: 28/03/2010, 12h18
  3. FFT et IFFT sur musique
    Par segatasan dans le forum Mathématiques du supérieur
    Réponses: 3
    Dernier message: 03/06/2009, 12h51
  4. [matlab] et ifft : comment ça marche exactement ?
    Par bratisla dans le forum Logiciel - Software - Open Source
    Réponses: 2
    Dernier message: 07/05/2009, 12h43
  5. phase x = fftshift(ifft(X)) avec MATLAB
    Par oliAudio dans le forum Logiciel - Software - Open Source
    Réponses: 0
    Dernier message: 24/06/2008, 11h50