Discussion :

[takout]

Bonjour à tous,

Dans le but de modéliser une courbe fermée, j'utilise la transformée de fourier discrète. J'aimerai ré-échantillonner cette courbe à l'aide de la librairie C++ FFTW3. Pour se faire, je détermine les coefficients de fourrier de la courbe et je calcule la transformée de fourier inverse en souhaitant obtenir en sortie un tableau de la taille de l'échantillon voulu. Cependant, la fonction de la transformée inverse que j'ai implémenté ne peut pas utiliser un tableau de taille différente que la courbe.
Comment puis-je faire une interpolation en utilisant cette librairie ?

Code C :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
 
#include <complex>
#include <fftw3.h>
 
void transformeeFourierComplex(complex<double> *curve, complex<double> *descriptor, int n)
{
 
  fftw_plan plan;
 
  plan = fftw_plan_dft_1d(n, reinterpret_cast<fftw_complex*>(curve), reinterpret_cast<fftw_complex*>(descriptor), FFTW_FORWARD, FFTW_ESTIMATE);
  fftw_execute(plan);
 
   /*liberation de la memoire*/
  fftw_destroy_plan(plan);
 
  return;
}
 
void transformeeFourierInv(complex<double> *descriptor, complex<double> *curve, int n)
{
 
  fftw_plan plan;
 
  plan = fftw_plan_dft_1d(n, reinterpret_cast<fftw_complex*>(descriptor), reinterpret_cast<fftw_complex*>(curve), FFTW_BACKWARD, FFTW_ESTIMATE);
  fftw_execute(plan);
 
  fftw_destroy_plan(plan);
}

[Ehouarn]

Bonjour,

Puisque le but est d'interpoler, je ne vois pas pourquoi passer par la transformée de Fourier inverse.
A partir du moment où tu as les coefficients de Fourier, tu peux réaliser ton interpolation en n'importe quel point, et donc construire ton signal ré-échantillonné.

[takout]

Hum je n'ai pas trop compris.
Je pensais que pour une interpolation il fallait déterminer les coefficients de Fourier puis passer par la transformée inverse pour trouver le points interpolé de la courbe.
Comment fait-tu uniquement avec les descripteurs de Fourier pour faire l'interpolation ?

J'ai recodé la transformée inverse mais mon souci c'est que mon implémentation est trop lente, c'est pour cela que je veux utiliser FFTW3.

[pseudocode]

Cependant, la fonction de la transformée inverse que j'ai implémenté ne peut pas utiliser un tableau de taille différente que la courbe.
Comment puis-je faire une interpolation en utilisant cette librairie ?

Il te faut augmenter la taille de ton tableau en insérant des zéros pour les fréquences hautes (pour lesquelles on n'a pas d'information dans le signal original).

google: fftw3 interpolation zero padding

[takout]

Merci Pseudocode,

Pour un sur-échantillonnage le zero-padding fonctionne mais tu fais comment pour un sous-échantillonnage ?

[pseudocode]

Merci Pseudocode,

Pour un sur-échantillonnage le zero-padding fonctionne mais tu fais comment pour un sous-échantillonnage ?

Réciproquement, tu raccourcis le tableau en supprimant les fréquences hautes avant de faire la transformée inverse.

Pour du sous-échantillonnage, je ne suis pas certain que passer par FFT/iFFT soit réellement plus performant que de faire une convolution dans le domaine temporel. Tout dépend de la taille du signal et du filtre utilisé.

[Ehouarn]

Bonjour,

Comment fait-tu uniquement avec les descripteurs de Fourier pour faire l'interpolation ?

Et bien si tu connais les coefficients de Fourier A_i, B_i , alors la fonction d'interpolation, valide sur tout l'intervalle est de la forme:
F(x) = Somme_i A_i cos(...k_i...x) + B_i sin(...k_i...x)
(la forme exacte dépend de l'intervalle considéré [-pi;pi] [0:2pi] et/ou si la représentation est sous forme réelle ou complexe)
Qu'il suffit alors d'appliquer pour les valeurs de x auxquelles tu souhaites ré-échantillonner.
Non?

[membreComplexe12]

bonjour tous,

je suis désolé de l'incruster dans la discussion mais j'ai une question
car je n'ai jamais vu d'interpolation par FFT, comment ça marche ?

-> tu fais une FFT sur ta courbe expérimentale afin d'avoir les amplitudes des différentes fréquences.

-> ensuite tu considère que les "pics" qui ont une amplitudes suffisante sont les plus représentatifs et tu les conserves (comment savoir le nombre de pics à conserver?)

Ce qui me gène dans cette histoire c'est que si tu as pas trop de points expérimentaux alors du va louper pas mal de frequences dans ton signal et donc au final ton résultat sera surement moins bon qu'avec une interpolation de type polynomiale ou spline ?

en gros quel est l'intéret de ce type d'interpolation ?

[pseudocode]

bonjour tous,

je suis désolé de l'incruster dans la discussion mais j'ai une question
car je n'ai jamais vu d'interpolation par FFT, comment ça marche ?

On considère que les données de départ sont le résultat de l'échantillonnage d'un signal continu réel (qui est inconnu).

Avec la transformée de Fourier discrète, on construit une fonction qui représente le signal continu.

Dans le domaine temporel, cette fonction f(t) est une somme de cosinus/sinus. On peut donc l'évaluer pour n'importe quelle valeur de "t", et donc faire un nouvel échantillonnage.

Cette technique s'appelle le resampling.

[membreComplexe12]

merci pour cette précision pseudocode! j'essaierai de faire ce genre d'interpolation un de ces quatre pour voir ce que ça donne

A+

[Ehouarn]

Bonjour,

Ce qui me gène dans cette histoire c'est que si tu as pas trop de points expérimentaux alors du va louper pas mal de frequences dans ton signal et donc au final ton résultat sera surement moins bon qu'avec une interpolation de type polynomiale ou spline ?

en gros quel est l'intéret de ce type d'interpolation ?

Il faut effectivement se poser la question de la précision des valeurs des points sur lesquels on base toute interpolation et du possible niveau de bruit associé.

Ceci dit, si on souhaite une interpolation aussi précise que possible, alors une interpolation via une FFT est bien meilleure (convergence exponentielle vers la solution en fonction du nombre de points utilisés) qu'une interpolation de type spline ou polynomiale (typiquement d'une convergence d'ordre 2 ou 3).

[membreComplexe12]

Bonjour,

Bonjour et merci d'avoir pris le temps de repondre, c'est très gentil.

(convergence exponentielle vers la solution en fonction du nombre de points utilisés) qu'une interpolation de type spline ou polynomiale (typiquement d'une convergence d'ordre 2 ou 3).

OK, je comprends tout à fait l'intérêt du coup.

Par contre, comment peut on montrer/demontrer ceci ? ça m'intéresse de vraiment comprendre en profondeur ce genre de chose.

merci

[pseudocode]

Par contre, comment peut on montrer/demontrer ceci ? ça m'intéresse de vraiment comprendre en profondeur ce genre de chose.

Tu peux trouver un début d'explication sur la page de wikipedia qui s'appelle "Convergence of Fourier series".

Pour plus d'information, il faut se pencher sur le vaste domaine de l'analyse harmonique.

[membreComplexe12]

Tu peux trouver un début d'explication sur la page de wikipedia qui s'appelle "Convergence of Fourier series".
Pour plus d'information, il faut se pencher sur le vaste domaine de l'analyse harmonique.

merci beaucoup !!! très intéressant.

entre temps, je viens de faire des tests avec la fonction interpft de matlab :

test numéro 1:

Code MATLAB :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
10
11
12
13
14
clear all;close all;clc;
%definition générales
facteur=100;
dt=1;
%signal echantillonné
t1=0:dt:10;
y1=sin(2*pi*5*t1)+sin(2*pi*10*t1);
Fe=length(t1)/(t1(end)-t1(1))
%reconstruction interpolation
t=0:dt/facteur:10;
y=interpft(y1,length(t));
%plot
figure (1)
plot(t1,y1,'*b',t,y,'r')

voici le résultat en PJ n°1 "interpolation1.jpg"

ici j'ai une fréquence échantillonnage de Fe=1.1Hz or comme mon signal à ses fréquences supérieures à Fshanon=2.2Hz alors
c'est normal que mon interpolation ne soit pas bonne, c'est bien cela ? :pff:

test numéro 2:

Code MATLAB :

Sélectionner tout - Visualiser dans une fenêtre à part

y1=sin(2*pi*0.05*t1)+sin(2*pi*0.1*t1);

cette fois j'ai remplacé la fonction par une fonction de fréquence plus basse et voici le résultat en PJ n°2 "interpolation2.jpg"

ici je suis un peu étonné que l'interpolation ne soit pas meilleure car j'ai une fréquence échantillonnage qui est 2*supérieur à la plus haute fréquence
de mon signal. Je m'attendais à une interpolation quasi parfaite.... ?

qu'es ce qui fait que ce ne soit pas meilleur ? que dois je faire pour que ce soit meilleur ?
Je pense qu'il faut augmenté de nombre de points expérimentaux et donc Fe mais ce n'est pas toujours possible...

[pseudocode]

Je pense qu'il faut augmenté de nombre de points expérimentaux et donc Fe mais ce n'est pas toujours possible...

Tu peux augmenter le nombre de points expérimentaux en ajoutant des zéros.

Ces points sont bien sur "fictifs", mais ils permettent de calculer plus de coefficients de Fourier , et donc d'avoir une meilleure précision.

[membreComplexe12]

merci beaucoup pseudo code, je crois avoir compris à présent

en fait avec ce "bourrage de zeros" ça devient une super méthode l'interpolation par FFT !

[membreComplexe12]

merci pseudo code mais un dernier blocage "pratique":

je viens d'essayer de faire du zero padding temporel avant de faire l'interpolation FFT et voici mon code MATLAB :

Code MATLAB :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
clear all;close all;clc;
%definition générales
facteur=100;
dt=1;
%signal echantillonné
t1=0:dt:10;
y1=sin(2*pi*5*t1)+sin(2*pi*10*t1);
Fe=length(t1)/(t1(end)-t1(1))
%reconstruction interpolation
t=0:dt/facteur:10;
%zero padding temporel
y1padding=[y1 zeros(1,10)];
FeAvecPadding=length(y1padding)/(t1(end)-t1(1))
y=interpft(y1padding,length(t));
%plot
figure (1)
plot(t1,y1,'*b',t,y,'r')

la fréquence d'echantillonage a été multipliée par 2, je m'attendais donc à un résultat meilleur que l'interpolation précédente mais ce n'est pas le cas...
j'obtient des oscillations

pourriez vous me dire comment corriger ceci dans mon code et m'expliquer pourquoi c'est faux ? (je pensais vraiment que c'était correct cette fois...)
merci beaucoup !!!

[Ehouarn]

Bonjour,

Par contre, comment peut on montrer/demontrer ceci ? ça m'intéresse de vraiment comprendre en profondeur ce genre de chose.

C'est à la base des méthodes spectrales (de résolution numérique d'équations différentielles), donc tu trouveras ces infos dans les livres/cours sur le sujet comme par exemple celui-ci

Pour ce qui est des exemples que tu montres (je n'utilises pas matlab; je n'ai pas regardé en détail), c'est plus que louche. Ton interpolation n'en est visiblement pas une puisque tu ne retrouves même pas les "points de contrôle" du signal initial...

[pseudocode]

la fréquence d'echantillonage a été multipliée par 2, je m'attendais donc à un résultat meilleur que l'interpolation précédente mais ce n'est pas le cas...
j'obtient des oscillations

pourriez vous me dire comment corriger ceci dans mon code et m'expliquer pourquoi c'est faux ? (je pensais vraiment que c'était correct cette fois...)
merci beaucoup !!!

1. Comme je l'ai dit post #4, pour faire une interpolation dans le domaine temporel il faut faire du zero-padding dans le domaine fréquentiel. Et réciproquement.

2. Ca me semble curieux que tous tes points expérimentaux soient dans le range 10E-14. Je ne connais pas trop Matlab, mais est-ce que par hasard tu n'aurais pas utilisé que des valeurs entières pour "t1". Auquel cas l'expression "y1=sin(2*pi*5*t1)+sin(2*pi*10*t1);" vaudrait tout le temps zéro.

[membreComplexe12]

Bonjour tous et merci de votre aide !!!

C'est à la base des méthodes spectrales (de résolution numérique d'équations différentielles), donc tu trouveras ces infos dans les livres/cours sur le sujet comme par exemple celui-ci

merci beaucoup pour le lien !

Pour ce qui est des exemples que tu montres (je n'utilises pas matlab; je n'ai pas regardé en détail), c'est plus que louche. Ton interpolation n'en est visiblement pas une puisque tu ne retrouves même pas les "points de contrôle" du signal initial...

une interpolation ne passe forcement par les points de contrôle... donc ça me parait pas illogique (exemple: une interpolation de moindre carré)

1. Comme je l'ai dit post #4, pour faire une interpolation dans le domaine temporel il faut faire du zero-padding dans le domaine fréquentiel. Et réciproquement.

Je pense avoir compris ceci. En faisant ceci on va conserver les fréquences du signal de départ mais on va avoir plus de point dans le temporel car on a fait du zeroPadding frequenciel....

Ici ce n'est pas vraiment ça ma question.

J'ai ce signal :

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

yExp=sin(2*pi*25*tExp)+sin(2*pi*120*tExp);

et une frequence d'echantillonage d'environ 10Hz car mon vecteur temps et défini comme ceci :

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
dt=0.1;
tExp=0:dt:10;

j'ai donc deux frequences 25 et 120Hz mais si ma frequence d'echantillonage
est de 10 Hz alors mon interpolation FFT ne sera pas capable de représenter les frequences au dessus de 5Hz si je m'en tiens au theoreme de Shanon ?
tu coup, je me dis que si je veux augmenter la frequence d'echantillonage alors il faut que je fasse du zero padding temporelle ? (avant de faire de l'interpolation FFT qui est fait du zeroPadding frequentiel)

2. Ca me semble curieux que tous tes points expérimentaux soient dans le range 10E-14. Je ne connais pas trop Matlab, mais est-ce que par hasard tu n'aurais pas utilisé que des valeurs entières pour "t1". Auquel cas l'expression "y1=sin(2*pi*5*t1)+sin(2*pi*10*t1);" vaudrait tout le temps zéro.

oui exactement, j'ai fais cette erreur dans mon message précédent