Discussion :

[Iven1]

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]

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]

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]

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]

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é]

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