fftn.c

来自「Audacity是一款用於錄音和編輯聲音的、免費的開放源碼軟體。它可以執行於Ma」· C语言 代码 · 共 1,047 行 · 第 1/2 页

      cd = sin(sd);
      cd = 2.0 * cd * cd;
      sd = sin(sd + sd);
      kk = 1;
      ii++;

      switch (factor [ii - 1]) {
       case 2:
     /* transform for factor of 2 (including rotation factor) */
     kspan /= 2;
     k1 = kspan + 2;
     do {
        do {
           k2 = kk + kspan;
           ak = Re [k2];
           bk = Im [k2];
           Re [k2] = Re [kk] - ak;
           Im [k2] = Im [kk] - bk;
           Re [kk] += ak;
           Im [kk] += bk;
           kk = k2 + kspan;
        } while (kk <= nn);
        kk -= nn;
     } while (kk <= jc);
     if (kk > kspan)
       goto Permute_Results_Label;		/* exit infinite loop */
     do {
        c1 = 1.0 - cd;
        s1 = sd;
        do {
           do {
          do {
             k2 = kk + kspan;
             ak = Re [kk] - Re [k2];
             bk = Im [kk] - Im [k2];
             Re [kk] += Re [k2];
             Im [kk] += Im [k2];
             Re [k2] = c1 * ak - s1 * bk;
             Im [k2] = s1 * ak + c1 * bk;
             kk = k2 + kspan;
          } while (kk < nt);
          k2 = kk - nt;
          c1 = -c1;
          kk = k1 - k2;
           } while (kk > k2);
           ak = c1 - (cd * c1 + sd * s1);
           s1 = sd * c1 - cd * s1 + s1;
           c1 = 2.0 - (ak * ak + s1 * s1);
           s1 *= c1;
           c1 *= ak;
           kk += jc;
        } while (kk < k2);
        k1 += inc + inc;
        kk = (k1 - kspan) / 2 + jc;
     } while (kk <= jc + jc);
     break;

       case 4:			/* transform for factor of 4 */
     ispan = kspan;
     kspan /= 4;

     do {
        c1 = 1.0;
        s1 = 0.0;
        do {
           do {
          k1 = kk + kspan;
          k2 = k1 + kspan;
          k3 = k2 + kspan;
          akp = Re [kk] + Re [k2];
          akm = Re [kk] - Re [k2];
          ajp = Re [k1] + Re [k3];
          ajm = Re [k1] - Re [k3];
          bkp = Im [kk] + Im [k2];
          bkm = Im [kk] - Im [k2];
          bjp = Im [k1] + Im [k3];
          bjm = Im [k1] - Im [k3];
          Re [kk] = akp + ajp;
          Im [kk] = bkp + bjp;
          ajp = akp - ajp;
          bjp = bkp - bjp;
          if (iSign < 0) {
             akp = akm + bjm;
             bkp = bkm - ajm;
             akm -= bjm;
             bkm += ajm;
          } else {
             akp = akm - bjm;
             bkp = bkm + ajm;
             akm += bjm;
             bkm -= ajm;
          }
          /* avoid useless multiplies */
          if (s1 == 0.0) {
             Re [k1] = akp;
             Re [k2] = ajp;
             Re [k3] = akm;
             Im [k1] = bkp;
             Im [k2] = bjp;
             Im [k3] = bkm;
          } else {
             Re [k1] = akp * c1 - bkp * s1;
             Re [k2] = ajp * c2 - bjp * s2;
             Re [k3] = akm * c3 - bkm * s3;
             Im [k1] = akp * s1 + bkp * c1;
             Im [k2] = ajp * s2 + bjp * c2;
             Im [k3] = akm * s3 + bkm * c3;
          }
          kk = k3 + kspan;
           } while (kk <= nt);

           c2 = c1 - (cd * c1 + sd * s1);
           s1 = sd * c1 - cd * s1 + s1;
           c1 = 2.0 - (c2 * c2 + s1 * s1);
           s1 *= c1;
           c1 *= c2;
           /* values of c2, c3, s2, s3 that will get used next time */
           c2 = c1 * c1 - s1 * s1;
           s2 = 2.0 * c1 * s1;
           c3 = c2 * c1 - s2 * s1;
           s3 = c2 * s1 + s2 * c1;
           kk = kk - nt + jc;
        } while (kk <= kspan);
        kk = kk - kspan + inc;
     } while (kk <= jc);
     if (kspan == jc)
       goto Permute_Results_Label;		/* exit infinite loop */
     break;

       default:
     /*  transform for odd factors */
#ifdef FFT_RADIX4
     fputs ("Error: " FFTRADIXS "(): compiled for radix 2/4 only\n", stderr);
     fft_free ();		/* free-up memory */
     return -1;
     break;
#else	/* FFT_RADIX4 */
     k = factor [ii - 1];
     ispan = kspan;
     kspan /= k;

     switch (k) {
      case 3:	/* transform for factor of 3 (optional code) */
        do {
           do {
          k1 = kk + kspan;
          k2 = k1 + kspan;
          ak = Re [kk];
          bk = Im [kk];
          aj = Re [k1] + Re [k2];
          bj = Im [k1] + Im [k2];
          Re [kk] = ak + aj;
          Im [kk] = bk + bj;
          ak -= 0.5 * aj;
          bk -= 0.5 * bj;
          aj = (Re [k1] - Re [k2]) * s60;
          bj = (Im [k1] - Im [k2]) * s60;
          Re [k1] = ak - bj;
          Re [k2] = ak + bj;
          Im [k1] = bk + aj;
          Im [k2] = bk - aj;
          kk = k2 + kspan;
           } while (kk < nn);
           kk -= nn;
        } while (kk <= kspan);
        break;

      case 5:	/*  transform for factor of 5 (optional code) */
        c2 = c72 * c72 - s72 * s72;
        s2 = 2.0 * c72 * s72;
        do {
           do {
          k1 = kk + kspan;
          k2 = k1 + kspan;
          k3 = k2 + kspan;
          k4 = k3 + kspan;
          akp = Re [k1] + Re [k4];
          akm = Re [k1] - Re [k4];
          bkp = Im [k1] + Im [k4];
          bkm = Im [k1] - Im [k4];
          ajp = Re [k2] + Re [k3];
          ajm = Re [k2] - Re [k3];
          bjp = Im [k2] + Im [k3];
          bjm = Im [k2] - Im [k3];
          aa = Re [kk];
          bb = Im [kk];
          Re [kk] = aa + akp + ajp;
          Im [kk] = bb + bkp + bjp;
          ak = akp * c72 + ajp * c2 + aa;
          bk = bkp * c72 + bjp * c2 + bb;
          aj = akm * s72 + ajm * s2;
          bj = bkm * s72 + bjm * s2;
          Re [k1] = ak - bj;
          Re [k4] = ak + bj;
          Im [k1] = bk + aj;
          Im [k4] = bk - aj;
          ak = akp * c2 + ajp * c72 + aa;
          bk = bkp * c2 + bjp * c72 + bb;
          aj = akm * s2 - ajm * s72;
          bj = bkm * s2 - bjm * s72;
          Re [k2] = ak - bj;
          Re [k3] = ak + bj;
          Im [k2] = bk + aj;
          Im [k3] = bk - aj;
          kk = k4 + kspan;
           } while (kk < nn);
           kk -= nn;
        } while (kk <= kspan);
        break;

      default:
        if (k != jf) {
           jf = k;
           s1 = pi2 / (double) k;
           c1 = cos(s1);
           s1 = sin(s1);
           if (jf > max_factors)
         goto Memory_Error_Label;
           Cos [jf - 1] = 1.0;
           Sin [jf - 1] = 0.0;
           j = 1;
           do {
          Cos [j - 1] = Cos [k - 1] * c1 + Sin [k - 1] * s1;
          Sin [j - 1] = Cos [k - 1] * s1 - Sin [k - 1] * c1;
          k--;
          Cos [k - 1] = Cos [j - 1];
          Sin [k - 1] = -Sin [j - 1];
          j++;
           } while (j < k);
        }
        do {
           do {
          k1 = kk;
          k2 = kk + ispan;
          ak = aa = Re [kk];
          bk = bb = Im [kk];
          j = 1;
          k1 += kspan;
          do {
             k2 -= kspan;
             j++;
             Rtmp [j - 1] = Re [k1] + Re [k2];
             ak += Rtmp [j - 1];
             Itmp [j - 1] = Im [k1] + Im [k2];
             bk += Itmp [j - 1];
             j++;
             Rtmp [j - 1] = Re [k1] - Re [k2];
             Itmp [j - 1] = Im [k1] - Im [k2];
             k1 += kspan;
          } while (k1 < k2);
          Re [kk] = ak;
          Im [kk] = bk;
          k1 = kk;
          k2 = kk + ispan;
          j = 1;
          do {
             k1 += kspan;
             k2 -= kspan;
             jj = j;
             ak = aa;
             bk = bb;
             aj = 0.0;
             bj = 0.0;
             k = 1;
             do {
            k++;
            ak += Rtmp [k - 1] * Cos [jj - 1];
            bk += Itmp [k - 1] * Cos [jj - 1];
            k++;
            aj += Rtmp [k - 1] * Sin [jj - 1];
            bj += Itmp [k - 1] * Sin [jj - 1];
            jj += j;
            if (jj > jf) {
               jj -= jf;
            }
             } while (k < jf);
             k = jf - j;
             Re [k1] = ak - bj;
             Im [k1] = bk + aj;
             Re [k2] = ak + bj;
             Im [k2] = bk - aj;
             j++;
          } while (j < k);
          kk += ispan;
           } while (kk <= nn);
           kk -= nn;
        } while (kk <= kspan);
        break;
     }
     /*  multiply by rotation factor (except for factors of 2 and 4) */
     if (ii == mfactor)
       goto Permute_Results_Label;		/* exit infinite loop */
     kk = jc + 1;
     do {
        c2 = 1.0 - cd;
        s1 = sd;
        do {
           c1 = c2;
           s2 = s1;
           kk += kspan;
           do {
          do {
             ak = Re [kk];
             Re [kk] = c2 * ak - s2 * Im [kk];
             Im [kk] = s2 * ak + c2 * Im [kk];
             kk += ispan;
          } while (kk <= nt);
          ak = s1 * s2;
          s2 = s1 * c2 + c1 * s2;
          c2 = c1 * c2 - ak;
          kk = kk - nt + kspan;
           } while (kk <= ispan);
           c2 = c1 - (cd * c1 + sd * s1);
           s1 += sd * c1 - cd * s1;
           c1 = 2.0 - (c2 * c2 + s1 * s1);
           s1 *= c1;
           c2 *= c1;
           kk = kk - ispan + jc;
        } while (kk <= kspan);
        kk = kk - kspan + jc + inc;
     } while (kk <= jc + jc);
     break;
#endif	/* FFT_RADIX4 */
      }
   }

/*  permute the results to normal order---done in two stages */
/*  permutation for square factors of n */
Permute_Results_Label:
   Perm [0] = ns;
   if (kt) {
      k = kt + kt + 1;
      if (mfactor < k)
    k--;
      j = 1;
      Perm [k] = jc;
      do {
     Perm [j] = Perm [j - 1] / factor [j - 1];
     Perm [k - 1] = Perm [k] * factor [j - 1];
     j++;
     k--;
      } while (j < k);
      k3 = Perm [k];
      kspan = Perm [1];
      kk = jc + 1;
      k2 = kspan + 1;
      j = 1;
      if (nPass != nTotal) {
/*  permutation for multivariate transform */
Permute_Multi_Label:
     do {
        do {
           k = kk + jc;
           do {
          /* swap Re [kk] <> Re [k2], Im [kk] <> Im [k2] */
          ak = Re [kk]; Re [kk] = Re [k2]; Re [k2] = ak;
          bk = Im [kk]; Im [kk] = Im [k2]; Im [k2] = bk;
          kk += inc;
          k2 += inc;
           } while (kk < k);
           kk += ns - jc;
           k2 += ns - jc;
        } while (kk < nt);
        k2 = k2 - nt + kspan;
        kk = kk - nt + jc;
     } while (k2 < ns);
     do {
        do {
           k2 -= Perm [j - 1];
           j++;
           k2 = Perm [j] + k2;
        } while (k2 > Perm [j - 1]);
        j = 1;
        do {
           if (kk < k2)
         goto Permute_Multi_Label;
           kk += jc;
           k2 += kspan;
        } while (k2 < ns);
     } while (kk < ns);
      } else {
/*  permutation for single-variate transform (optional code) */
Permute_Single_Label:
     do {
        /* swap Re [kk] <> Re [k2], Im [kk] <> Im [k2] */
        ak = Re [kk]; Re [kk] = Re [k2]; Re [k2] = ak;
        bk = Im [kk]; Im [kk] = Im [k2]; Im [k2] = bk;
        kk += inc;
        k2 += kspan;
     } while (k2 < ns);
     do {
        do {
           k2 -= Perm [j - 1];
           j++;
           k2 = Perm [j] + k2;
        } while (k2 > Perm [j - 1]);
        j = 1;
        do {
           if (kk < k2)
         goto Permute_Single_Label;
           kk += inc;
           k2 += kspan;
        } while (k2 < ns);
     } while (kk < ns);
      }
      jc = k3;
   }

   if ((kt << 1) + 1 >= mfactor)
     return 0;
   ispan = Perm [kt];
   /* permutation for square-free factors of n */
   j = mfactor - kt;
   factor [j] = 1;
   do {
      factor [j - 1] *= factor [j];
      j--;
   } while (j != kt);
   kt++;
   nn = factor [kt - 1] - 1;
   if (nn > max_perm)
     goto Memory_Error_Label;
   j = jj = 0;
   for (;;) {
      k = kt + 1;
      k2 = factor [kt - 1];
      kk = factor [k - 1];
      j++;
      if (j > nn)
    break;				/* exit infinite loop */
      jj += kk;
      while (jj >= k2) {
     jj -= k2;
     k2 = kk;
     k++;
     kk = factor [k - 1];
     jj += kk;
      }
      Perm [j - 1] = jj;
   }
   /*  determine the permutation cycles of length greater than 1 */
   j = 0;
   for (;;) {
      do {
     j++;
     kk = Perm [j - 1];
      } while (kk < 0);
      if (kk != j) {
     do {
        k = kk;
        kk = Perm [k - 1];
        Perm [k - 1] = -kk;
     } while (kk != j);
     k3 = kk;
      } else {
     Perm [j - 1] = -j;
     if (j == nn)
       break;		/* exit infinite loop */
      }
   }
   max_factors *= inc;
   /*  reorder a and b, following the permutation cycles */
   for (;;) {
      j = k3 + 1;
      nt -= ispan;
      ii = nt - inc + 1;
      if (nt < 0)
    break;			/* exit infinite loop */
      do {
     do {
        j--;
     } while (Perm [j - 1] < 0);
     jj = jc;
     do {
        kspan = jj;
        if (jj > max_factors) {
           kspan = max_factors;
        }
        jj -= kspan;
        k = Perm [j - 1];
        kk = jc * k + ii + jj;
        k1 = kk + kspan;
        k2 = 0;
        do {
           k2++;
           Rtmp [k2 - 1] = Re [k1];
           Itmp [k2 - 1] = Im [k1];
           k1 -= inc;
        } while (k1 != kk);
        do {
           k1 = kk + kspan;
           k2 = k1 - jc * (k + Perm [k - 1]);
           k = -Perm [k - 1];
           do {
          Re [k1] = Re [k2];
          Im [k1] = Im [k2];
          k1 -= inc;
          k2 -= inc;
           } while (k1 != kk);
           kk = k2;
        } while (k != j);
        k1 = kk + kspan;
        k2 = 0;
        do {
           k2++;
           Re [k1] = Rtmp [k2 - 1];
           Im [k1] = Itmp [k2 - 1];
           k1 -= inc;
        } while (k1 != kk);
     } while (jj);
      } while (j != 1);
   }
   return 0;			/* exit point here */

   /* alloc or other problem, do some clean-up */
Memory_Error_Label:
   fputs ("Error: " FFTRADIXS "() - insufficient memory.\n", stderr);
   fft_free ();			/* free-up memory */
   return -1;
}
#endif	/* _FFTN_C */
/* ---------------------- end-of-file (c source) ---------------------- */