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
| #include <stdio.h>
#include <math.h>
#include "complex.h"
#include "tfr1d.h"
// Cette unite gere tous les calculs lies a la TFR1D.
// Par rapport au TP initial, TOUTES les procedures sont maintenant parametrees
// dynamiquement, a l'execution, et non plus statiquement a la compilation.
// C'est le role du parametre "int n" terminal des procedures, qui remplace
// partout la constante 'N' du sujet.
// NOTE IMPORTANTE : AUCUNE ALLOCATION DE MEMOIRE N'EST FAITE.
// En consequence, l'utilisateur doit allouer LUI-MEME l'espace necessaire
// pour stocker le spectre calcule.
// L'utilisation du type "complexe" permet de definir des tableaux a une seule
// dimension pour l'utilisateur, meme si, en fait, ce sont des matrices.
// Afin d'accelerer le plus possible les calculs, il a ete fait un usage
// intensif des decalages binaires. Ceci a permis, entre autres, d'eliminer
// totalement tout calcul explicite d'une quelconque puissance de 2 lors
// de l'implementation de l'algorithme "papillon". De meme, le parametre
// initial 'P' a ete totalement supprime explicitement, meme si de nombreuses
// variables jouent son role ('j' dans la boucle la plus externe de la fonction
// "tfr_hard", par exemple).
// Conversion d'une tableau reel en tableau complexe.
void convert ( float signal[], complex spectre[], int n) {
int i ;
for (i=0;i<n;i++)
spectre[i]=c_algb(signal[i],0);
}
// Affiche un tableau de complexes.
void display_c ( complex spectre[], int n ) {
int i ;
for (i=0;i<n;i++)
printf("(%f,%f) ",spectre[i].Real,spectre[i].Imag);
printf("\n\n");
}
// Affiche les MODULES d'un tableau de complexes.
void display_m ( complex spectre[], int n ) {
int i ;
for (i=0;i<n;i++)
printf("%f ",c_mod(spectre[i]));
printf("\n\n");
}
// Procedure demandee dans le sujet.
void affiche ( complex spectre[], int n ) {
int i ;
for (i=0;i<n;i++)
printf("|(%f,%f)|=%f\n",spectre[i].Real,spectre[i].Imag,c_mod(spectre[i]));
printf("\n\n");
}
/*
La methode ci-dessous, bien que plus economique en temps machine, perd de la
precision a la 6eme decimale. Le probleme pourrait etre regle en passant les
calculs internes en "double" au lieu de "float", afin de garder la precision.
void init_w ( complex w[], int n , int sens ) {
complex p ;
int i ;
p=c_trig(1,-2*sens*M_PI/n) ;
n>>=1;
w[0]=c_algb(1,0);
for (i=1;i<n;i++)
w[i]=c_mul(w[i-1],p);
}
*/
// Solution "lente", mais legerement plus precise.
// Au lieu d'un seul calcul sin/cos, on en a maintenant n/2...
void init_w ( complex w[], int n , int sens ) {
int i ;
float t ;
n>>=1; // On divise n par 2, cela allege tous les calculs.
for (i=0;i<n;i++) {
t=-i*sens*M_PI;
w[i]=c_trig(1,t/n); // On a DEJA divise n par 2 !
}
}
// Calcule le miroir d'un chiffre i binaire sur n bits.
int shuffle ( int i, int n ) {
int k,r ;
// En fait, via decalages, on "pousse" le chiffre dans un buffer,
// mais dans l'ordre inverse.
for (k=n,r=0;k!=1;k>>=1,i>>=1)
r=(r << 1) | (i & 1) ;
return r ;
}
// Une fois la procedure "shuffle" ecrite, il n'y a plus de problemes.
void init_spectre ( complex signal[], complex spectre[] , long n ) {
long i ;
for (i=0;i<n;i++)
spectre[i]=signal[shuffle(i,n)];
}
// Calcul reel de TFR/TFR-1, n'initialise RIEN, par contre, et considere que
// ses parametres SONT initialises. Cette separation est necessaire pour
// permettre la reutilisation du code par TFR2D.
void tfr_hard ( complex spectre[], complex w[] , int n , int sens ) {
int i,j,k,x,y,z,t ;
complex a ;
// Notations :
// i : Nombre d'elements par bloc.
// j : Nombre de blocs.
// k : Moitie de i, par soucis de rapidite.
for (i=2,j=n>>1;j;i<<=1,j>>=1) {
k=i>>1;
// Phase d'additions/soustractions.
// Notations :
// x : Numero du bloc en cours
// y : Indice courant (global)
for (x=0,y=0;x<j;x++) {
// On effectue 'k' operations doubles.
for (z=0;z<k;z++,y++) {
// Il faut "sauver" au moins une des valeurs initiales.
a=spectre[y];
spectre[y]=c_add(a,spectre[y+k]);
spectre[y+k]=c_sub(a,spectre[y+k]);
}
y+=k; // On "saute" le deuxieme demi-bloc, il est deja traite.
// Phase de multiplications.
// On multiplie sur les blocs IMPAIRS, il suffit donc de tester
// le bit de poids faible de x pour savoir si une multiplication
// est necessaire.
// Notations :
// z : Indice courant (global)
// t : Puissance de Wn utilis‚e
if (x & 1)
// On incremente t en fonction du nombre de blocs, mais
// en divisant d'abord ce dernier par 2.
for (z=y-i,t=0;z<y;z++,t+=j>>1)
spectre[z]=c_mul(spectre[z],w[t]);
}
}
if (sens==TFRDI) {
a=c_algb(n,0);
for (i=0;i<n;i++)
spectre[i]=c_div(spectre[i],a);
}
}
// Calcule la TFR d'un signal discret.
void tfr1D ( complex signal[], complex spectre[], int n , int sens ) {
// Le tableau W est local.
complex w[n/2] ;
// Initialisation des puissances de Wn et du spectre.
init_w(w,n,sens);
init_spectre(signal,spectre,n);
// Appel du "vrai" calcul de TFR. Noter le parametre "w" supplementaire.
tfr_hard(spectre,w,n,sens);
} |
Transformée de fourier
Envoyé par rstaros
Répondre avec citation
Partager