Discussion :

[curlcie]

Bonjour à tous,

J'essaie depuis plusieurs jours de construire mon algorithme de la fft

mais j'obtiens des valeurs incohérentes : par exemple, quand je lui mets en entrée du cosinus échantillonné, j'obtiens en sortie des valeurs dépassant 1000 ! Je ne trouve vraiment pas où est mon erreur.

Si quelqu'un pouvait m'aider sur ce point, ça serait bien gentil de sa part!

Merci beaucoup!

ps : je me suis aidé de cette page : http://www.docstoc.com/docs/20472205...RETE-%28TFD%29

( la fonction attendre() n'est ni plus ni moins que _getch() )

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
#include "divers.h"
#include "principal.h"
#include "math.h"
 
// fonction qui réorganise les données :
int reverseBits(int x, int nbits)
{
int temp = 0;
for (int i = 0; i < nbits; i++)
{
	int v = 0;
	if (x & (1 << i))
	{
		int oppindex = abs(i - nbits + 1);
		int diff = oppindex - i;
		v = (1 << (i + diff));			
	}
	temp = temp + v;
}
return temp;
}
 
// fonction puissance avec des double 
double puissance_double (double a, int b)
{
	double c=1;
	if (b==0) {return 1;}
	else
	{
		for (int i=1; i<=b; i++)
		{
			c=c*a ;
		}
		return c;
	}
}
 
 
 
// fonction puissance avec des int :
int puissance_int (int a, int b)
{
	int c=1;
	if (b==0) {return 1;}
	else
	{
		for (int i=1; i<=b; i++)
		{
			c=c*a ;
		}
		return c;
	}
}
 
void main ( void ) {
// Lit le fichier de données et l'enregistre dans le tableau E. 
// Partie réelle dans la première colonne et partie imaginaire dans la seconde.
 
fstream F ;
double point ;
int N = 0 ;
double E[512][2] ;
 
F.open("donnees",ios::in) ;
 
while(!F.eof())
{
	F>>point ;
	if (!F.fail())
	{
		E[N][0] = point ;
		E[N][1] = 0 ;
		N = N+1 ;
	}
 
}
F.close() ;
 
// déclaration des variables :
 
int n = 9 ;
int k, p, i, M, m, d, l ;
double S[512][2], Wb[1][2], W[1][2] ;
double Pi = 3.141592 ;
 
// reorganisation des données :
for ( i=0 ; i<N ; i++)
{
	S[i][0] = E[reverseBits(i, n)][0] ;
}
for (i=0 ; i<N ; i++)
{
	E[i][0] = S[i][0] ; 
	E[i][1] = 0.0 ;
}
 
for (i=0 ; i<N ; i++)
{
	S[i][0] = 0.0 ; 
	S[i][1] = 0.0 ; 
}
 
// transformée :
 
for (k=1; k<=n;k++) // pour chaque étage
{
	p = puissance_int(2, k-1) ;   //taille des papillons
	Wb[0][0] = cos(Pi/(2*p)) ;   // coefficient de base
	Wb[0][1] = sin(Pi/(2*p)) ;   
	M = puissance_int(2,n-k) ;   //  nombre de paquets de papillons
 
	for (m=1 ; m<=M ; m++) // pour chaque paquet
	{
		W[0][0] = 1.0 ;
		W[0][1] = 0.0 ;
		d = (m-1)*(N/M) ;  // indice de début de paquet
 
 
		for (i=1; i<=p; i++)  // pour chaque papillon
		{
			l = d+i-1;  // premier point du papillon
			E[l+p][0] = E[l+p][0]*W[0][0] - E[l+p][1]*W[0][1] ;
			E[l+p][1] = E[l+p][1]*W[0][0] + E[l+p][0]*W[0][1] ;
 
			S[l][0] = E[l][0] + E[l+p][0] ;
			S[l][1] = E[l][1] + E[l+p][1] ;
 
			S[l+p][0] = E[l][0] - E[l+p][0] ;
			S[l+p][1] = E[l][1] - E[l+p][1] ;
 
			W[0][0] = W[0][0]*Wb[0][0] - W[0][1]*Wb[0][1] ;
			W[0][1] = W[0][1]*Wb[0][0] + W[0][0]*Wb[0][1] ;
		}
	}
	for (int h=0 ; h<N ; h++)
	{	E[h][0] = S[h][0] ;
		E[h][1] = S[h][1] ;	}
}
 
for (i=0 ; i<=4 ; i++)
{cout<<S[i][0]<<endl;}
attendre();
 
 
// calcul du module :
for (i=0 ; i<N ; i++)
{
	S[i][0] = puissance_double(S[i][0],2) + puissance_double(S[i][1],2) ;
}
 
 
// Écriture dans le fichier :
fstream a ;
a.open("spectre",ios::out) ;
 
for (i=0;i<N;i++)
{
	a<<S[i][0]<<endl ;
}
 
a.close() ;
 
}

[diogene]

Je n'ai pas regardé le détail (il y aurait beaucoup d'améliorations à faire) mais il doit y avoir une erreur dans le calcul du papillon ici:

Code :

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

1
2
E[l+p][0] = E[l+p][0]*W[0][0] - E[l+p][1]*W[0][1] ;
E[l+p][1] = E[l+p][1]*W[0][0] + E[l+p][0]*W[0][1] ;

La deuxième expression utilise E[l+p][0] qui a été modifié par le calcul de la première expression au lieu d'utiliser la valeur d'origine.

[curlcie]

Oui merci, je n'avais pas vu!

Grâce à ta remarque j'ai pu trouver une deuxième erreur un peu plus bas, avec les nouveaux W[][].

Mais, même en modifiant, j'obtiens des valeurs farfelues!

Je commence à me demander si c'est pas trop compliqué pour moi...

Quelqu'un voit-il une autre erreur?

Merci

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
// fonction qui réorganise les données :
int reverseBits(int x, int nbits)
{
int temp = 0;
for (int i = 0; i < nbits; i++)
{
    int v = 0;
    if (x & (1 << i))
    {
        int oppindex = abs(i - nbits + 1);
        int diff = oppindex - i;
        v = (1 << (i + diff));            
    }
    temp = temp + v;
}
return temp;
}
 
// fonction puissance avec des double 
double puissance_double (double a, int b)
{
    double c=1;
    if (b==0) {return 1;}
    else
    {
        for (int i=1; i<=b; i++)
        {
            c=c*a ;
        }
        return c;
    }
}
 
 
 
// fonction puissance avec des int :
int puissance_int (int a, int b)
{
    int c=1;
    if (b==0) {return 1;}
    else
    {
        for (int i=1; i<=b; i++)
        {
            c=c*a ;
        }
        return c;
    }
}
 
void main ( void ) {
// Lit le fichier de données et l'enregistre dans le tableau E. 
// Partie réelle dans la première colonne et partie imaginaire dans la seconde.
 
fstream F ;
double point ;
int N = 0 ;
double E[32][2] ;
 
F.open("donnees",ios::in) ;
 
while(!F.eof())
{
    F>>point ;
    if (!F.fail())
    {
        E[N][0] = point ;
        E[N][1] = 0 ;
        N = N+1 ;
    }
 
}
F.close() ;
 
// déclaration des variables :
N=32;
int n = 5 ;
int k, p, i, M, m, d, l ;
double S[32][2], Wb[1][2], W[1][2] ;
double Pi = 3.141592 ;
 
// reorganisation des données :
for ( i=0 ; i<N ; i++)
{
    S[i][0] = E[reverseBits(i, n)][0] ;
}
for (i=0 ; i<N ; i++)
{
    E[i][0] = S[i][0] ; 
    E[i][1] = 0.0 ;
}
 
for (i=0 ; i<N ; i++)
{
    S[i][0] = 0.0 ; 
    S[i][1] = 0.0 ; 
}
 
// transformée :
 
for (k=1; k<=n;k++) // pour chaque étage
{
    p = puissance_int(2, k-1) ;   //taille des papillons
    Wb[0][0] = cos(Pi/(2*p)) ;   // coefficient de base
    Wb[0][1] = sin(Pi/(2*p)) ;   
    M = puissance_int(2,n-k) ;   //  nombre de paquets de papillons
 
    for (m=1 ; m<=M ; m++) // pour chaque paquet
    {
        W[0][0] = 1.0 ;
        W[0][1] = 0.0 ;
        d = (m-1)*(N/M) ;  // indice de début de paquet
 
 
        for (i=1; i<=p; i++)  // pour chaque papillon
        {
            l = d+i-1 ;  // premier point du papillon
 
            double temp[1][2] ;
 
            temp[0][0] = E[l+p][0] ;
            E[l+p][0] = E[l+p][0]*W[0][0] - E[l+p][1]*W[0][1] ;
            E[l+p][1] = E[l+p][1]*W[0][0] + temp[0][0]*W[0][1] ;
 
            S[l][0] = E[l][0] + E[l+p][0] ;
            S[l][1] = E[l][1] + E[l+p][1] ;
 
            S[l+p][0] = E[l][0] - E[l+p][0] ;
            S[l+p][1] = E[l][1] - E[l+p][1] ;
 
            temp[0][0] = W[0][0] ;
            W[0][0] = W[0][0]*Wb[0][0] - W[0][1]*Wb[0][1] ;
            W[0][1] = W[0][1]*Wb[0][0] + temp[0][0]*Wb[0][1] ;
        }
    }
    for (int h=0 ; h<N ; h++)
    {    E[h][0] = S[h][0] ;
        E[h][1] = S[h][1] ;    }
}
 
// calcul du module :
for (i=0 ; i<N ; i++)
{
    S[i][0] = puissance_double(S[i][0],2) + puissance_double(S[i][1],2) ;
}
 
 
// Écriture dans le fichier :
fstream a ;
a.open("spectre",ios::out) ;
 
for (i=0;i<N;i++)
{
    a<<S[i][0]<<endl ;
}
 
a.close() ;
 
}

[curlcie]

C'est bon!!! j'ai réussi!!

je vous montre le code pour ceux que ça interesse

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
#include "divers.h"
#include "principal.h"
#include "math.h"
 
// fonction qui réorganise les données :
int reverseBits(int x, int nbits)
{
int temp = 0;
for (int i = 0; i < nbits; i++)
{
    int v = 0;
    if (x & (1 << i))
    {
        int oppindex = abs(i - nbits + 1);
        int diff = oppindex - i;
        v = (1 << (i + diff));            
    }
    temp = temp + v;
}
return temp;
}
 
// fonction puissance avec des double 
double puissance_double (double a, int b)
{
    double c=1;
    if (b==0) {return 1;}
    else
    {
        for (int i=1; i<=b; i++)
        {
            c=c*a ;
        }
        return c;
    }
}
 
 
 
// fonction puissance avec des int :
int puissance_int (int a, int b)
{
    int c=1;
    if (b==0) {return 1;}
    else
    {
        for (int i=1; i<=b; i++)
        {
            c=c*a ;
        }
        return c;
    }
}
 
void main ( void ) {
// Lit le fichier de données et l'enregistre dans le tableau E. 
// Partie réelle dans la première colonne et partie imaginaire dans la seconde.
 
fstream F ;
double point ;
int N = 0 ;
double E[32][2] ;
 
F.open("donnees",ios::in) ;
 
while(!F.eof())
{
    F>>point ;
    if (!F.fail())
    {
        E[N][0] = point ;
        E[N][1] = 0 ;
        N = N+1 ;
    }
 
}
F.close() ;
 
// déclaration des variables :
N=32;
int n = 5 ;
int k, p, i, M, m, d, l ;
double S[32][2], Wb[1][2], W[1][2] ;
double Pi = 3.1415926535897932384626433832795 ;
 
// reorganisation des données :
for ( i=0 ; i<N ; i++)
{
    S[i][0] = E[reverseBits(i, n)][0] ;
}
for (i=0 ; i<N ; i++)
{
    E[i][0] = S[i][0] ; 
    E[i][1] = 0.0 ;
}
 
for (i=0 ; i<N ; i++)
{
    S[i][0] = 0.0 ; 
    S[i][1] = 0.0 ; 
}
 
// transformée :
 
for (k=1; k<=n;k++) // pour chaque étage
{
    p = puissance_int(2, k-1) ;   //taille des papillons
    M = puissance_int(2,n-k) ;   //  nombre de paquets de papillons
 
    for (m=1 ; m<=M ; m++) // pour chaque paquet
    {
        d = (m-1)*2*p ;  // indice de début de paquet
 
 
        for (i=1; i<=p; i++)  // pour chaque papillon
        {
            l = d+i-1 ;  // premier point du papillon
 
            double temp[1][2] ;
 
            temp[0][0] = E[l+p][0]*cos(Pi*l/(puissance_int(2,k-1))) + E[l+p][1]*sin(Pi*l/(puissance_int(2,k-1))) ;
            temp[0][1] = E[l+p][1]*cos(Pi*l/(puissance_int(2,k-1))) - E[l+p][0]*sin(Pi*l/(puissance_int(2,k-1))) ;
 
            S[l][0] = E[l][0] + temp[0][0] ;
            S[l][1] = E[l][1] + temp[0][1] ;
 
            temp[0][0] = E[l+p][0]*cos(Pi*(l+p)/(puissance_int(2,k-1))) + E[l+p][1]*sin(Pi*(l+p)/(puissance_int(2,k-1))) ;
            temp[0][1] = E[l+p][1]*cos(Pi*(l+p)/(puissance_int(2,k-1))) - E[l+p][0]*sin(Pi*(l+p)/(puissance_int(2,k-1))) ;
 
            S[l+p][0] = E[l][0] + temp[0][0] ;
            S[l+p][1] = E[l][1] + temp[0][1] ;
 
        }
    }
    for (int h=0 ; h<N ; h++)
    {    E[h][0] = S[h][0] ;
        E[h][1] = S[h][1] ;    }
}
 
for (i=0 ; i<N ; i++)
{cout<<S[i][0]<<"  + i*   "<<S[i][1]<<endl;}
attendre();
 
 
// calcul du module au carré :
for (i=0 ; i<N ; i++)
{
    S[i][0] = puissance_double(S[i][0],2) + puissance_double(S[i][1],2) ;
}
 
 
// Écriture dans le fichier :
fstream a ;
a.open("spectre",ios::out) ;
 
for (i=0;i<N;i++)
{
    a<<S[i][0]<<endl ;
}
 
a.close() ;
 
}