C: Recettes numériques (FFT)

Question

Cette question s'adresse à tous les fans de recettes numériques ou toute personne qui comprend bien FFT.

Quelqu'un peut-il expliquer pourquoi la composante réelle est calculée par -2 * (sin (thêta / 2)) ^ 2? Je ne peux pas sembler envelopper ma tête autour d'elle. Je l'ai vu d'autres exemples tels que http://www.dspdimension.com/admin / DFT-a-pied / tutoriel qui prend simplement cos (thêta) comme réel et -sin (thêta) comme imaginaire. J'ai aussi vu ici en base http://www.dspguide.com/ch12/3. htm qu'il énumère comme cos (thêta) comme réel et -sin (thêta) comme imaginaire. Je peux penser à un peu plus de ressources qui prennent simplement le cos et -sin comme réelle et imaginaire.

cos(theta) = 1-2*(sin(theta/2))^2

si l'identité trigonométrique ci-dessus est vrai, alors pourquoi ne pas ce folllow?

theta=isign*(6.28318530717959/mmax);
wtemp=sin(0.5*theta);
wpr = -2.0*wtemp*wtemp;
wpi=sin(theta);

Je présume que la recette numérique doit utiliser une identité trigonométrique? Je ne peux pas sembler figurer dehors et le livre n'explique pas du tout.

code trouvé ici: http: //ronispc.chem.mcgill .ca / Ronis / chem593 / sinfft.c.html

#define SWAP(a,b) tempr=(a);(a)=(b);(b)=tempr

void four1(double *data,unsigned long nn,int isign)
{
        unsigned long n,mmax,m,j,istep,i;
        double wtemp,wr,wpr,wpi,wi,theta;
        double tempr,tempi;

        n=nn << 1;
        j=1;
        for (i=1;i<n;i+=2) {
                if (j > i) {
                        SWAP(data[j],data[i]);
                        SWAP(data[j+1],data[i+1]);
                }
                m=n >> 1;
                while (m >= 2 && j > m) {
                        j -= m;
                        m >>= 1;
                }
                j += m;
        }
        mmax=2;
        while (n > mmax) {
                istep=mmax << 1;
                theta=isign*(6.28318530717959/mmax);
                wtemp=sin(0.5*theta);
                wpr = -2.0*wtemp*wtemp;
                wpi=sin(theta);
                wr=1.0;
                wi=0.0;
                for (m=1;m<mmax;m+=2) {
                        for (i=m;i<=n;i+=istep) {
                                j=i+mmax;
                                tempr=wr*data[j]-wi*data[j+1];
                                tempi=wr*data[j+1]+wi*data[j];
                                data[j]=data[i]-tempr;
                                data[j+1]=data[i+1]-tempi;
                                data[i] += tempr;
                                data[i+1] += tempi;
                        }
                        wr=(wtemp=wr)*wpr-wi*wpi+wr;
                        wi=wi*wpr+wtemp*wpi+wi;
                }
                mmax=istep;
        }
}
#undef SWAP

Answer 1

Départ:

• cos (A + B) = cos (A) cos (b) - sin (a) sin (B)
• sin (a + b) = sin (A) cos (b) + cos (A) sin (B)
• cos (2A) = 1 - 2 sin ² (A)
• e ^{i θ} = cos (θ) + i sin (θ)

e ^{i (φ + δ)}

= cos (φ + δ) + i sin (φ + δ)

= cos (φ) cos (δ) - sin (φ) sin (δ) + i [sin (φ) cos (A) + cos (φ) sin (δ)]

= cos (φ) [1 - 2 sin ² (δ / 2)] + i sin (φ) [1 - 2 sin ² (δ / 2 )] + i sin (δ) [i * sin (φ) + cos (φ)]

= [cos (φ) + i sin (φ)] [1 - 2 sin ² (δ / 2)] + [cos (φ) + i sin (φ)] i sin (δ)

= e ^{i φ} + e ^{i φ} [- 2 sin ² (δ / 2) + i sin (δ)]

Modifier : C'était beaucoup de formatage inutile de ma part. Il est en fait plus simple:

y ^{(a + b)} = y ^a x y ^b pour tout y. Donc:

e ^{i (φ + δ)}

= e ^{i φ} e ^{i δ}

= e ^{i φ} [cos (δ) + i sin (δ)]

= e ^{i φ} [1 - 2 sin ² (δ / 2) + i sin (δ)]

Answer 2

Une forme de l'identité de demi-angle de cosinus est:

cos(theta) = 1 - 2*(sin(theta/2)^2)

Je ne sais pas si cela répond à votre question.

Answer 3

La raison est la précision numérique. Si vous examinez de près le code suivant:

wtemp=sin(0.5*theta);
wpr = -2.0*wtemp*wtemp;
wpi=sin(theta);

et

wr=(wtemp=wr)*wpr-wi*wpi+wr;
wi=wi*wpr+wtemp*wpi+wi;

ils sont conçus pour travailler ensemble pour obtenir le résultat attendu correct. serait mis en œuvre une approche naïve comme suit:

wpr = cos(theta);
wpi = sin(theta);

et

wtemp = wr;
wr =wr*wpr - wi*wpi;
wi =wi*wpr + wtemp*wpi;

et avec une précision infinie serait fonctionnellement équivalent.

Toutefois, lorsque theta est proche de zéro (par exemple un grand nombre de points d'échantillonnage ou grand nn), cos(theta) devient problématique, car pour les petits angles cos(theta) approches 1 et a une pente approchant 0. Et à un cos(theta) == 1 angle. Mon expérience montre que, pour float cos(2*PI/N) == 1 exactement N >= 25736 pour flotteur (à savoir de précision 32 bits). Une FFT point 25736 est concevable. Donc, pour éviter ce problème wpr est calculé comme cos(theta)-1 en utilisant la formule de demi-angle de la trigonométrie. Il est calculé avec sin qui est très précis pour les petits angles, donc pour les petits angles à la fois wpr et wpi sont petites et précises. Ceci est alors annulée dans le code de mise à jour en ajoutant le dos à dos 1 après la multiplication complexe. Exprimé mathématiquement, nous obtenons:

w_p = cos(theta) - 1    + i*sin(theta) 
    = -2*sin^2(theta/2) + i*sin(theta)

et la règle de mise à jour est:

w = w * (w_p + 1) 
  = w + w*w_p

Answer 4

Je ne sais pas FFT bien que je l'ai fait une, mais sa fait longtemps.

Je regardais les identités trig http: //www.sosmath. com / TRIG / trig5 / trig5 / trig5.html

et de la première identité "produit-To-Somme" pour le péché (u) * sin (v) nous avons

^ sin 2 (theta / 2) = (1/2) [cos (zéro) - cos (theta)] = 0,5 - 0,5 cos (theta)

Est-ce que cette aide?

Answer 5

Qu'est-ce que la confusion est que NR utilise la version Math / Physique de la FFT qui fait tourner le tripoter facteurs dans le sens opposé qui définissent la FFT EEs. Ainsi, l'avant NR est l'inverse EE et inversement. Notez que dans NR vers l'avant a un exposant positif au lieu de l'exposant négatif EE. La méthode de EE transforme le temps de fréquence où la version Math et Physique joue avec un moment angulaire.

Answer 6

Ils utilisent des identités trigonométriques pour minimiser le nombre de fonctions circulaires dont ils ont besoin pour calculer. De nombreuses implémentations rapides il suffit de regarder jusqu'à ces fonctions. Si vous voulez vraiment savoir que vous devez travailler simplement les détails par la boucle au-dessus de dérouler et de faire des substitutions variables appropriées .... oui fastidieux.

BTW, la mise en œuvre de NR est connu pour être très lent.

Paul

Answer 7

Ok, donc ici est l'identité trigonométrique. La raison ne est pas la moitié cos (thêta) identité trigonométrique est à cause de la dépendance à l'égard euler et nombres imaginaires. Le calcul est toujours au-delà de moi ...

_{(source: librow.com )}