probleme avec la FFT

**Iven1** · 24/05/2009, 14h32

Salut tt le monde

je travail avec ce code fft:

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
#define PI 3.141592
#define SWAP(a,b) tempr=(a);(a)=(b);(b)=tempr
 
/* ******** FFT ********************************/
void FFT(float data[], int N, int isign)
{
 
    int n,mmax,m,j,istep,i;
    float wtemp,wr,wpr,wpi,wi,theta,tempr,tempi;
 
    n=N * 2;
 
 
    j=0;
    for (i=0;i<n/2;i+=2) {
        if (j > i) {
 
            SWAP(data[j],data[i]);
 
            SWAP(data[j+1],data[i+1]);
 
            if((j/2)<(n/4)){
 
                SWAP(data[(n-(i+2))],data[(n-(j+2))]);
 
                SWAP(data[(n-(i+2))+1],data[(n-(j+2))+1]);
            }
        }
        m=n/2;
        while (m >= 2 && j >= m) {
            j -= m;
            m = m/2;
        }
        j += m;
    }
 
    mmax=2;
    //external loop
    while (n > mmax)
    {
        istep = mmax<<  1;
        theta=isign*(6.2831853/mmax);
        wtemp=sin(0.5*theta);
        wpr = -2.0*wtemp*wtemp;
        wpi=sin(theta);
        wr=1.0;
        wi=0.0;
        //internal loops
        for (m=0;m<mmax;m+=2) {
            for (i= m;i<=n;i+=istep) {
                j=i+mmax;
                tempr=wr*data[j-1]-wi*data[j];
                tempi=wr*data[j]+wi*data[j-1];
                data[j-1]=data[i-1]-tempr;
                data[j]=data[i]-tempi;
                data[i-1] += tempr;
                data[i] += tempi;
            }
            wr=(wtemp=wr)*wpr-wi*wpi+wr;
            wi=wi*wpr+wtemp*wpi+wi;
        }
        mmax=istep;
    }
    }

Quand je l'appel par cette ligne:

Code :

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

FFT(Fft,t,1);

tq Fft un vecteur de données de F valeurs, et quand j'affiche les résultats je vois que certains case ne changent pas, elles concervent leur valeurs que celles de l'entrée.

- Le vecteur Fft ( allocation mémoire) contient des valeurs (float).
- Cette ligne me gène : SWAP(data[j+1],data[i+1]); je la comprend pas.
- F est de puissance de 2
- t = 1024
- je pense que c'est une fft reel mais quelque part je vois qu'elle n'est pas spéciale de la partie reel

merci d'avance.

**diogene** · 24/05/2009, 17h58

Cette ligne me gène : SWAP(data[j+1],data[i+1]); je la comprend pas.

Elle utilise la macro SWAP, définie plus haut et qui procède à l'échange des deux éléments du tableau. Cette partie est en fait pour classer les éléments du tableau au préalable dans l'ordre des index codés avec les bits inversés (pour permettre de faire la FFT sur place) .

- je pense que c'est une fft reel mais quelque part je vois qu'elle n'est pas spéciale de la partie reel

En regardant le code

Code :

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

1
2
3
4
5
6
                tempr=wr*data[j-1]-wi*data[j];
                tempi=wr*data[j]+wi*data[j-1];
                data[j-1]=data[i-1]-tempr;
                data[j]=data[i]-tempi;
                data[i-1] += tempr;
                data[i] += tempi;

on peut voir que ce n'est PAS une FFT qui porte sur les nombres réels. Les données de départ sont des complexes rangés dans le tableau de float de la façon suivante :
partie réelle point 0, partie imaginaire point 0, partie réelle point 1, partie imaginaire point 1,...

Il reste à savoir si N est le nombre d'éléments complexes ou le nombre d'éléments du tableau (soit deux fois plus) et j'ai la paresse d'étudier les détails du code pour le savoir.

**Iven1** · 24/05/2009, 18h28

Envoyé par diogene

En regardant le code

Code :

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

1
2
3
4
5
6
                tempr=wr*data[j-1]-wi*data[j];
                tempi=wr*data[j]+wi*data[j-1];
                data[j-1]=data[i-1]-tempr;
                data[j]=data[i]-tempi;
                data[i-1] += tempr;
                data[i] += tempi;

on peut voir que ce n'est PAS une FFT qui porte sur les nombres réels. Les données de départ sont des complexes rangés dans le tableau de float de la façon suivante :
partie réelle point 0, partie imaginaire point 0, partie réelle point 1, partie imaginaire point 1,...

Il reste à savoir si N est le nombre d'éléments complexes ou le nombre d'éléments du tableau (soit deux fois plus) et j'ai la paresse d'étudier les détails du code pour le savoir.

Merci infiniment
J'ai attendu quelqu'un me confirme ce qui est en marron me confirme que cette fft traite aussi la partie imaginaire que je ne l'ai pas besoin. Quoi je doit faire pour la transformer en fft traite que de la partie réelle, si j'ai deux tableaux l'un pour la parie réelle et l'autre pour la partie imaginaire, le probleme ne reste pas posé; j'annule l'imaginaire et ça se fait.

Je vais concentrer sur le code surtout pour le N, je vois pas la boucle qui traite tout le tableau, est ce que elle le traite par bloc de 1024 ou commment!
Merci une autre fois, diogene,et j'aimerais avoir encore d'aide.

**diogene** · 24/05/2009, 19h44

Il y a deux formes habituellement utilisées pour les données complexes
- Alterner dans un même tableau partie réelle et partie imaginaire (c'est ton cas)
- Utiliser un tableau pour les parties réelles et un autre pour les parties imaginaires

Lorsque le signal est réel, on peut se dispenser de la partie imaginaire (puisque elle est connue comme nulle). Toutefois, la Transformée de Fourier Discrète (TFD) d'un signal réel est un signal complexe. Il comporte donc à priori deux fois plus de données (float) que les données de départ. On utilise alors la propriété que la TFD d'un signal réel a une partie réelle symétrique et une partie imaginaire antisymétrique pour ne stocker que la moitié des données. Naturellement, le programme de calcul est différent de celui, plus général, qui calcule la TFD complexe d'un signal complexe.

Si, c'est ton cas, tu as l'algorithme général, il faut mettre à 0 explicitement les parties imaginaires. Donc tu construit ton tableau de la façon suivante :
partie réelle point 0, 0, partie réelle point 1, 0,...

est ce que elle le traite par bloc de 1024 ou commment!

Non, il traite la TFD pour N égal à une puissance de 2, pas par bloc de 1024

**Iven1** · 24/05/2009, 20h51

diogene;
J'ai essayé de bien comprendre le code fft, j'arrive à comprendre le N; moi j'ai un tableau (Fft) de (F=f*1024) valeurs je vais tout les traiter pas seulement 1024 points, alors j'ai integré cette boucle:

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
 
 
    int  d=0;
      for( d=0;d<f;d++){
 
     for( i=0;i<1024;i++){
        Tableau[2*i]=Fft[(1024*d)+i];
        Tableau[2*i+1]=0;}
 
 
 
        FFT(Tableau,1024,1);
 
        for ( k=0;k<1024;k++){
            Fft[(1024*d)+k]=Tableau[2*k];}
 
 
       }

Je vois maintenant un changement de resultat et je vais verifier si celle que je veut.
Merci bien diogene

Invité · 24/05/2009, 23h01

Bonjour,

Ton code fft ressemble furieusement à celui de Numerical Recipes in C (jusqu'aux noms des variables...). Si tu as ce livre, c'est à la page 507 de la deuxième édition (sinon je crois qu'on peut le trouver en ligne).

A priori, ce code traite bien des nombres complexes, saisis sous la forme indiquée par Diogène. Maintenant, dans Num Rec, ils ont l'étrange habitude de n'avoir que des tableaux en base 1, ton code a été modifié pour fonctionner avec des tableaux en base zéro... Oui mais...

Code :

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

1
2
3
4
5
 
if((j/2)<(n/4)){
      SWAP(data[(n-(i+2))],data[(n-(j+2))]);
      SWAP(data[(n-(i+2))+1],data[(n-(j+2))+1]);
}

Est un ajout au code original, qui semble tirer parti d'une symétrie de ta fonction d'origine... Il faudrait regarder ce que c'est.

Code :

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

1
2
3
4
5
6
7
8
9
10
11
12
13
        for (m=0;m<mmax;m+=2) {
            for (i= m;i<=n;i+=istep) {
                j=i+mmax;
                tempr=wr*data[j-1]-wi*data[j];
                tempi=wr*data[j]+wi*data[j-1];
                data[j-1]=data[i-1]-tempr;
                data[j]=data[i]-tempi;
                data[i-1] += tempr;
                data[i] += tempi;
            }
            wr=(wtemp=wr)*wpr-wi*wpi+wr;
            wi=wi*wpr+wtemp*wpi+wi;
        }

Cette boucle me semble fautive : lors de la première itération, tu fais référence à data[i-1] alors que i vaut zéro... Ce ne va pas bien se passer.

Pour porter le code de NumRec dans des tableaux à base zéro, il te faudrait :

- soit avoir ci dessus

Code :

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

for (m=1;m<=mmax;m+=2) {

soit remplacer dans toute la boucle data[j-1] par data[j], data[j] par data[j+1] et pareil pour les i.

Je n'ai pas testé, mais je pense que ca devrait aller...

Francois

**Iven1** · 24/05/2009, 23h27

Merci Francois pour ta réponse.

J'ai trouvé ce code dans ce forum et ça me fait plaisir si tu poste le code original d'apres le livre ( Numerical Recipes in C).

pour ce code:

Code :

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

1
2
3
4
if((j/2)<(n/4)){
      SWAP(data[(n-(i+2))],data[(n-(j+2))]);
      SWAP(data[(n-(i+2))+1],data[(n-(j+2))+1]);
}

Je ne sais pas vraiment à quoi sert.

Invité · 24/05/2009, 23h42

Le passage correspondant de Numerical Recipes (2nd Ed) est disponible en ligne là

http://www.fizyka.umk.pl/nrbook/c12-2.pdf

(l'ensemble du livre est là http://www.fizyka.umk.pl/nrbook/bookcpdf.html)

Maintenant, il y a une nouvelle version de ce livre, qui date de 2007. Sur google books (http://books.google.com/books?id=1aA...um=4#PPA613,M1), on peut en voir un bout, et on y trouve ce code avec des data[j-1]. Je soupconne donc que ton code est une adaptation incorrecte du code de la Troisième édition.

Si tu l'as trouvé sur Developpez, il devait y avoir une référence, non?

Je vais certainement acheter ce livre (je te conseille de faire de même, programmer des calculs sans l'avoir lu, c'est pas sérieux). Je regarderai quand je l'aurai reçu (fin de la semaine, je pense, le temps que monsieur amazon fasse son boulot).

Francois

**Iven1** · 24/05/2009, 23h49

Je te remercie beaucoup Monsieur Francois

Invité · 15/06/2009, 20h09

Comme promis, j'ai jeté un coup d'oeil sur Numerical Recipes 3 (reçu aujourd'hui).

C'est très probablement la source principale d'inspiration de ce code (ou alors l'inventeur de cette autre fonction a choisi exactement les mêmes noms de variables, déclarées dans le même ordre).

Il y a clairement une erreur de recopie dans la boucle centrale :

Code :

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

1
2
for (m=0;m<mmax;m+=2) {
            for (i= m;i<=n;i+=istep) {

devrait être

Code :

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

1
2
for (m=1;m<mmax;m+=2) {
            for (i= m;i<=n;i+=istep) {

Le reste semble bon. Pour la première partie, c'est un ajout du copiste, qui modifie la boucle initiale de "bit reversal" et dépliant deux passages dans la boucle en un seul. Je ne pense pas que cela nuise, mais il ne parait pas clair que ca gagne du temps non plus...

Donc, voila, en changeant la valeur de départ de 0 en 1, on doit retrouver, à peu de choses près, le code de Numerical Recipes, qui est certainement bon.

Francois

probleme avec la FFT

C

Discussions similaires

Partager

Partager