Bonsoir.
J'essaie d'effectuer une comparaison entre un algorithme basique de transformée de Fourier discrète et la fonction FFT de MATLAB. L'objectif étant de comparer les temps d'exécution respectifs pour une utilisation dans un programme avec de très nombreuses TF à réaliser. Pour cela je pars d'une gaussienne centrée autour de 0 et je lance les 2 méthodes. Ensuite j'effectue une transformée de Fourier inverse avec les 2 méthodes. Logiquement je dois revenir à la gaussienne d'origine. Cela fonctionne parfaitement pour l'amplitude de la gaussienne dans les 2 cas. Mais en ce qui concerne la phase j'ai une erreur avec la FFT (qui n’existe pas avec la TFD) sur la phase qui devrait être nulle mais qui ne l'est pas tout à fait.
Existe-t-il une personne avertie qui pourrait me dire où je commets une erreur sur l'exploitation de la FFT et éventuellement me donner une solution pour la corriger ?
Merci d'avance car cela commence à me prendre beaucoup de temps et d'énergie.
Noel.
Voici le code que j'utilise :
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
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
 
clc;
clear all 
close all                                                                                           
 
%% Génération de la gaussienne temporelle
dt=4.17e-12;
Ndt=1000; % nombre d'échantillons temporels sur la gaussienne
tmax=0e-9; % décalage de la gaussienne sur l'axe du temps
lmh=1e-9;  % largeur à mi hauteur de la Gaussienne 
alpha=4*log(2)/(lmh)^2; %Facteur de rapidité                                                                        
E=1;       % Amplitude max de la gaussienne
t_G(:,1)=(-(Ndt-1)*dt/2:dt:(Ndt-1)*dt/2)';   % Axe de temps de la Gaussienne générée                                                                                                                                                                      
% t_G(:,1)=(0:dt:Ndt*dt)';   % Axe de temps de la Gaussienne générée                                                                                                                                                                      
Gauss(:,1)=E.*exp(-alpha*((t_G(:,1))).^2);
[pks locs]=findpeaks(Gauss);
Nechmaxgauss=locs; %recherche de l'échantillon où se trouve le max de la guassienne centrée
Gauss(:,1)=E.*exp(-alpha*((t_G(:,1)-tmax)).^2);
[pks locs]=findpeaks(Gauss);
Nechtmax=locs; %recherche de l'échantillon où se trouve le max de la gaussienne non centrée
Gauss(:,1)=Gauss(:,1)./max(Gauss(:,1)); %normalisation de la Gaussienne  
X0(:,2)=Gauss(:,1); 
X0(:,1)=t_G(:,1);
 
%% initialisation constantes avant TF
fmin=-1.50e9;
fmax=1.5e9;
nf=501;
df=(fmax-fmin)/(nf-1);
Ndf=round(1/(dt*df));
Ndfmin=round(fmin/df);
Ndfmax=round(fmax/df);
freq = (Ndfmin:Ndfmax)*df;
 
%%   Transformée de Fourier discrète de base
tic;
i=1;   % compteur pas fréquentiel
j=1;   % compteur pas temporel
while (i < nf+1)
    AC = 0;
    AS = 0;
    j=1;
    while (j < Ndt+1)
        U=1;
        wt=2*pi*freq(i)*X0(j,1);
        if j==1
           U=0.5;
        end 
        if j==Ndt
           U=0.5;
        end
        if abs(wt)<=1e-02
           C1=(1-wt*wt/2)+(wt^4/24)+(wt^6/720);     % Développement limité du cos autour de 0
           S1=(wt -(wt^3/6))+(wt^5/120)-(wt^5/5040);
           AC = AC + X0(j,2)*C1*U;
           AS = AS - X0(j,2)*S1*U;
        else
           AC = AC + X0(j,2)*cos((wt))*U;
           AS = AS - X0(j,2)*sin((wt))*U;
        end
        j=j+1;
    end
 
    reel(i)=AC*dt;
    img(i)=AS*dt;
 
    i=i+1;
end
 
Arg_TF_X1(:,1)=atan(img(:)./reel(:)).*180/pi;
Mod_TF_X1(:,1)=(reel(:).^2+img(:).^2).^0.5;
toc;
 
 
%% FFT
% fftw('planner','patient');
tic;
 
TF_X0(:,1)=fft(X0(:,2),Ndf)*dt; % calcul de la fft avec le même pas que la TF classique (ici df)
TF_X0(:,1)=fftshift(TF_X0(:,1)); % centrage de la fft sur f=0
% % Arg_TF_X0(:,1)=atan(imag(TF_X0(:,1))./real(TF_X0(:,1))).*180/pi; % calcul de l'argument de la FFT en degrés
% % Mod_TF_X0(:,1)=abs(TF_X0(:,1)); % calcul du module de la FFT
freqfft = (0:(Ndf-1))*(1/(Ndf)); % calcul évolution de la fréquence de la fft pour recentrer la phase lorsque des temps et/ou des fréquences négatifs sont présents
Phifft_decal=2*pi*freqfft'*Nechmaxgauss; % calcul de la phase de la fft pour des signaux comportant des temps et/ou des fréquences négatifs
Phasefft_shift=complex(cos(Phifft_decal),sin(Phifft_decal)); % décalage de la phase
% Arg_Phase_shift=atan(imag(Phasefft_shift)./real(Phasefft_shift)).*180/pi; % calcul de l'argument de la TF en degrés
TF_X0(:,1)=TF_X0(:,1).*Phasefft_shift;
TF_BP_X0(:,1)=TF_X0((Ndf/2+(Ndfmin)+1:Ndf/2+Ndfmax+1),1); % Limitation de la FFT aux fréquences comprises entre fmin et fmax comme pour la transformée classique
Arg_TF_BP_X0(:,1)=atan(imag(TF_BP_X0(:,1))./real(TF_BP_X0(:,1))).*180/pi; % calcul de l'argument de la TF en degrés
Mod_TF_BP_X0(:,1)=abs(TF_BP_X0(:,1)); % calcul du module de la TF
toc;
MethodFFT=fftw('planner');
 
 
%% FFT inverse
% Arg_TF_BP_X0(:,1)=Arg_TF_X1(:,1)
% t_G(:,1)=(-(Ndt-1)*dt/2:dt:(Ndt-1)*dt/2)';   % Axe de temps de la Gaussienne générée
tic;
tmin=-2.082915e-9;
% tmin=0;
tmax=2.082915e-9;
nt=1000;
dt=(tmax-tmin)/(nt-1);
% Ndf=round(1/(dt*df));
Ndtmin=round(tmin/df);
Ndtmax=round(tmax/df);
temp= (Ndtmin:Ndtmax)*dt;
tempfft=(1:(Ndf))*(1/(Ndf)); 
 
evoltemp=ifftshift(ifft((TF_BP_X0),Ndf));
[pks locs]=findpeaks(Mod_TF_BP_X0);
NechmaxTFgauss=locs; % numéro échantillon constituant l'origine de la phase si différent de 1
Phitemp_decal=-2*pi*(tempfft')*(NechmaxTFgauss-1); % calcul de la phase de la fft inverse
Phasetemp_shift=complex(cos(Phitemp_decal),sin(Phitemp_decal)); % décalage de la phase
Arg_Phase_shift=atan(imag(Phasetemp_shift)./real(Phasetemp_shift)).*180/pi; % calcul de l'argument de la TF en degrés
evoltemp=evoltemp.*Phasetemp_shift;
% Arg_evoltemp=atan(imag(evoltemp)./real(evoltemp)).*180/pi; % calcul de l'argument de la TF en degrés
% Arg_evoltemp_bp(:,1)=Arg_evoltemp((-Ndt/2+(Ndf/2)+1:Ndt/2+Ndf/2),1); % Limitation de la transformée de Fourier inverse aux temps compris entre tmin et tmax
evoltemp_bp(:,1)=evoltemp((-Ndt/2+(Ndf/2)+1:Ndt/2+Ndf/2),1); % Limitation de la transformée de Fourier inverse aux temps compris entre tmin et tmax
Arg_evoltemp_bp=0+atan(imag(evoltemp_bp)./real(evoltemp_bp)).*180/pi; % calcul de l'argument de la TF en degrés
evoltemp_bp=abs(evoltemp_bp)/dt;
toc;
 
%%   Transformée de Fourier discrète inverse de base
tic;
 
fr(:,1)=freq;                                                                                                                                               
data(:,1)=reel+i*img;
 
% ----- Paramètres frequentiels -------------
nbre_F=length(fr);
Df=abs(fr(2,1)-fr(1,1));                                                                    
nbre_T=nt;
% ----- Création du vecteur de temps --------
clear time;
j=1;                                                                                            
while (j < nbre_T+1)
    time(j,1)=((j-1)*(tmax-tmin)/(nbre_T-1))+tmin;                                              
    j=j+1;
end
% ----------
 
k=1;                                                                                            
j=1;                                                                                            
%h_wait=waitbar(0,'Transformée de Fourier inverse en cours...');
while (j < nbre_T+1)
    Amplitude=0;                                                                               
    k=1;
    while (k < nbre_F+1)
        U=1;
        var=2*pi*fr(k,1)*time(j,1);
        if k==1
            U=0.5;
        end 
        if k==nbre_F
            U=0.5;
        end
        if abs(var)<=1e-02
            C1=(1-var*var/2)+(var^4/24)+(var^6/720);                                           
            S1=(var-(var^3/6))+(var^5/120)-(var^5/5040);                                       
            Amplitude=Amplitude+data(k,1)*(C1+i*S1)*U;
        else
            Amplitude=Amplitude+data(k,1)*(cos(var)+i*sin(var))*U;
        end
        k=k+1;
    end
    Partie_reelle=real(Amplitude)*Df;
    Partie_imaginaire=imag(Amplitude)*Df;
    phi_T(j,1)=atan(Partie_imaginaire./Partie_reelle);
    argH_T(j,1)=phi_T(j,1)*180/pi;                                                                                                                                    
    mod_signal(j,1)=sqrt(Partie_reelle*Partie_reelle+Partie_imaginaire*Partie_imaginaire);     
    j=j+1;
    %waitbar(j/(nbre_T+1),h_wait);
end
 
Arg_GaussTF(:,1)=argH_T;
Mod_GaussTF(:,1)=mod_signal;
toc;
 
%% figures
figure(1)
subplot(211)
plot(freq,(Arg_TF_BP_X0(:,1)),'b');
hold on;
plot(freq,(Arg_TF_X1(:,1)),'r');
 
subplot(212)
plot(freq,(Mod_TF_BP_X0(:,1)),'b');
hold on;
plot(freq,(Mod_TF_X1(:,1)),'r');
 
figure(2)
plot(t_G,evoltemp_bp,'b');
hold on;
plot(t_G,Gauss,'r');
hold on;
plot(t_G,Mod_GaussTF,'g')
 
figure(3)
plot(Arg_evoltemp_bp,'b');
hold on;
plot(Arg_GaussTF,'r');
% plot(Gauss)