Discussion :

[toto_madit]

Bonjour,
je souhaiterais utiliser la librairie FFTW3 (la fameuse) pour traiter une série de profils temporels de grande taille (plus de 5000 points par profils et plus de 20k profils par traitement en double real complexes). Le déroulement du traitement pour chaque profil est plutôt simple: un opération de convolution effectuée dans l'espace de Fourier sous la forme d'un produit, puis FFT inverse du résultat. Il s'avère que je n'arrive pas à effectuer la transformée inverse (message d'erreur). Le problème viendrait de la définissions du plan pour effectuer la TF inverse. Mon programme fonctionne correctement tant que je n’essaie pas d'appliquer la TF inverse.

Est ce que quelqu'un pourrait, pour que je puisse m'en inspirer, me montrer un exemple de traitement FFT qui comprendrait:

une FFT direct d'un signal et sa transformée inverse dans le même programme (FORTRAN 90)?

par avance merci!

[toto_madit]

Bon, je crois que ça n'inspire pas grand monde, si je trouve la solution je la mettrai en message ..
@bientôt

[toto_madit]

Bonsoir tout le monde,

suite à mon message, je poste ici la solution que j'ai pu trouver
il s'agit d'un code fortran90 qui utilise la fftw3 pour calculer la transformé de fourier de plusieurs signaux types (nul, porte, chirp, sinusoïde) avec la possibilité de bruiter le signal. Tout est codé en Hard (code à modifier) pour choisir quel signal on veux traiter. Je suis intéressé par d'éventuels feedback (erreurs?).

Pour rappel, la TF de la librairie fftw3 se déroule en 3 étapes:
1- détermination d'un plan de calcul: call dfftw_plan_dft_1d
2- le calcule de la TF direct ou inverse: call dfftw_execute
3- la suppression du plan de calcul: dfftw_destroy_plan
ce code permet en outre de gérer automatiquement le zéro-padding pour ajuster la longueur du signal à 2**n le plus proche (supérieur).

Une remarque, dans les exemples courants donnés pour l'utilisation de cette librairie, la définition du plan de calcul se fait APRES l'initiation du tableau d'entré (ici IN). Or j'ai remarqué que dans cet ordre , la variable IN est systématiquement réinitialisé à 0. Dans le code présent, j'initialise mon tableau APRES la définition du plan de calcul pour que ça fonctionne. Est-ce normal?

Une autre question: lors de l'affichage du résultat à l'écran (l 160: "abs(out(n-i+1))/n"), je dois inverser l'ordre de la variable "out" qui contient la TF pour que les fréquences soient bien calées... est-ce aussi normal?

Bonne lecture et bon tests

Code :

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
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
program test_fftw3
implicit none
 
!notes								!
!code testant la fftw3 sur plusieurs signaux type (nul, rectangulaire, sinusoïdale, chirp)!
!pas de garanties ni explicite ni implicite!
!code libre GNU!
 
!compilation dans répertoire courant:
!#>gfortran fftw_test.f90 -lfftw3 -lm -o fftw_test
!#>./fftw_test
 
integer :: n, i
integer :: FFTW_FORWARD=1, FFTW_ESTIMATE=0, FFTW_BACKWARD=-1
integer*8 :: plan,plan2
logical	:: period,rect,chirp,bruit,nul
real*8	:: f0,f1,PI=3.141592,duree,freq_ech,t,x,y,fact,A0,A1
double complex, allocatable :: in(:),in2(:),out(:),dsp(:)
 
!choix du signal (rect=1/0 pour rectangulaire, period=1/0 pour sinusoïdale, bruit=1/0 pour signal bruité)
 
!signal nul
 
nul=.FALSE.
 
!rectangulaire------
 
rect=.FALSE.
 
!sinusoïdale--------
 
period=.TRUE.
 
!fréquences d'oscillation (Hz)
 
if(period)then
	A0=1.	!amplitude 0
	A1=0.5	!amplitude 1
	f0=10. 	!fréquence 0
	f1=15.	!fréquence 1
endif
 
!chirp
 
chirp=.FALSE.
 
!fréquence centrale
 
if(chirp)then
	f0=100.
endif
 
!bruit-------
 
bruit=.TRUE.
 
!amplitude bruit
 
if(bruit)then
	fact=0.2
	x=-1.
	y=1.
endif
 
! durée de l'enregistrement, ici 10 secondes.
 
duree=1.
 
! fréquence d'échantillonnage: nombre de points par unité de temps (Hz)
 
freq_ech=10000. 
 
!longueur du profil en nombre de points
 
n=floor(duree*freq_ech) 
 
!allocation mémoire in, in2, out
 
allocate(in(n)) !signal d'entré
allocate(in2(n)) !signal d'entré reconstruit à partir de la fft
allocate(out(n)) ! fft
allocate(dsp(n)) ! densité de puissance spectrale
 
!définition du plan de calcul pour la transformé directe et inverse
 
call dfftw_plan_dft_1d(plan,n,in,out,FFTW_FORWARD,FFTW_ESTIMATE)
call dfftw_plan_dft_1d(plan2,n,out,in2,FFTW_BACKWARD,FFTW_ESTIMATE)
 
!------------------------------
!---initialisation du tableau in
!------------------------------
 
do i=1,n
 
	!calcul du temps écoulé
 
	t=i/freq_ech
 
	!signal nul
 
	if (nul) then
		in(i)=cmplx(0.,0.)
	endif
 
	!signal rectangulaire
 
	if(rect) then
		if(t .ge. duree/10. .and. t .le. 2*duree/10 ) then
			in(i) = cmplx(1,0)
		else
			in(i) = cmplx(0,0)
		endif
	endif
 
	!signal périodique exp(i2PI*f0*t)+exp(i2PI*f1*t)
 
	if(period) then
    	in(i) = A0*cmplx(cos(2*PI*f0*t),sin(2*PI*f0*t))
    	in(i) = in(i)+A1*cmplx(cos(2*PI*f1*t),sin(2*PI*f1*t))
    endif
 
    !signal chirp exp(i2PI*f0*t**2)
 
	if(chirp) then
    	in(i) = cmplx(cos(2*PI*f0*t*t),sin(2*PI*f0*t*t))
    endif
 
    !introduction de bruit aléatoire
 
    if(bruit) then
    	in(i) = in(i)+fact*((rand(0)*(y-x))+x)
    endif
 
enddo
 
!--------------------------
!calcul de la transformé de fourier
!utilisation de la librairie libre fftw3
!--------------------------
 
! Forward Fourier transform
 
call dfftw_execute(plan)
call dfftw_destroy_plan(plan)
 
! Inverse Fourier transform
 
call dfftw_execute(plan2)
call dfftw_destroy_plan(plan2)
 
!!---------------------------
 
!calcul de la dsp (non normalisée)
 
dsp=out*conjg(out)
 
! sortie écran: fréquence, temps, abs(fft), real(DSP), real(signal), real(signal reconstruit)
do i=1,n
  write(*,*) i*freq_ech/n, i/freq_ech, abs(out(n-i+1))/n, real(dsp(n-i+1)/n**2), real(in(i)), real(in2(i))/n
enddo
 
deallocate(in,in2,out,dsp)
 
end program test_fftw3

[toto_madit]

Bonsoir,
après un petit moment d'absence, me voilà de retour. Est ce que quelqu'un peut apporter une précision ou une réponse concernant les deux questions abordées précédemment? Toute aide serait la bienvenue.
cordialement.

[toto_madit]

bonsoir,
Je viens de lire ici ce qui semble être la réponse à ma première question. A savoir, l'utilisation de la fonction qui permet d'initialiser le plan de calcul écrase les tableaux d'entré et de sortie et les initialise à 0. Reste à comprendre le deuxième point qui reste un mystère pour moi. Cette librairie ne me parait pas facile à utiliser surtout pour les adeptes du Fortran. Quelques retours d'utilisateurs? (en espérant que je ne sois pas tout seul à chercher à utiliser FFTW3 en Fortran...).
+