[beast/devel: 1/6] BSE: Make gslfft results compatible with fftw.



commit 6e1d04dc25f6f9d30b6b0c0159a14c27fff6b4dc
Author: Stefan Westerfeld <stefan space twc de>
Date:   Wed Oct 13 11:20:35 2010 +0200

    BSE: Make gslfft results compatible with fftw.
    
    The Numerical Recipies FFT and FFTW behave differently because a minus
    sign in the formula is different (so Numerical Recipies analysis performs
    the same computation like FFTW synthesis and vice versa). This change
    makes gslfft and FFTW results equivalent. Also, while FFTW never scales
    results, gslfft used to scale results during backward transform. To get
    the same results, scaling is now only performed optionally in gslfft.
    
    Finally, real variants (like fftar/fftsr) have been changed to produce the
    same results like real FFTW.

 bse/gsl-fftconf.sh |  150 ++++++++++++++++++++++++++++++++++++++--------------
 bse/gsl-fftgen.pl  |  104 +++++++++++++++++++++++++++---------
 bse/gslfft.h       |   85 ++++++++++++++++++++++++------
 3 files changed, 257 insertions(+), 82 deletions(-)
---
diff --git a/bse/gsl-fftconf.sh b/bse/gsl-fftconf.sh
index 109c66e..f4f97bf 100755
--- a/bse/gsl-fftconf.sh
+++ b/bse/gsl-fftconf.sh
@@ -29,7 +29,7 @@ $MKFFT $OPTIONS 0
 #
 # generate small fft sizes, seperating stage 2
 #
-for dir in --analysis --synthesis ; do
+for dir in --analysis --synthesis --synthesis-scale ; do
   DOPT="$OPTIONS --skip-macros $dir"
   echo "Generating FFT functions: $dir" >&2
   $MKFFT $DOPT	   2 F				# standalone fft2
@@ -68,31 +68,44 @@ cat <<__EOF
 
 /* public FFT frontends and generic handling of huge fft sizes */
 
+typedef enum {
+  BIG_ANALYSIS = 1,
+  BIG_SYNTHESIS = 2,
+  BIG_SYNTHESIS_SCALE = 3
+} FFTMode;
 
 static void
 gsl_power2_fftc_big (const unsigned int n_values,
 		     const $IEEE_TYPE  *rivalues_in,
 		     $IEEE_TYPE        *rivalues,
-                     const int          esign)
+                     FFTMode            mode)
 {
+  const int esign = (mode == BIG_ANALYSIS) ? -1 : 1;
   const unsigned int n_values2 = n_values << 1;
   double theta = esign < 0 ? -3.1415926535897932384626433832795029 : 3.1415926535897932384626433832795029;
   unsigned int i, block_size = 8192 << 1;
   double last_sin;
 
-  if (esign > 0)
+  if (mode == BIG_SYNTHESIS)
     {
       if (rivalues_in)
-        bitreverse_fft2analysis (n_values, rivalues_in, rivalues);
+        bitreverse_fft2synthesis (n_values, rivalues_in, rivalues);
       for (i = 0; i < n_values; i += 8192)
-        gsl_power2_fft8192analysis_skip2 (rivalues + (i << 1), rivalues + (i << 1));
+        gsl_power2_fft8192synthesis_skip2 (rivalues + (i << 1), rivalues + (i << 1));
     }
-  else
+  else if (mode == BIG_SYNTHESIS_SCALE)
     {
       if (rivalues_in)
-        bitreverse_fft2synthesis (n_values, rivalues_in, rivalues);
+        bitreverse_fft2synthesis_scale (n_values, rivalues_in, rivalues);
       for (i = 0; i < n_values; i += 8192)
-        gsl_power2_fft8192synthesis_skip2 (rivalues + (i << 1), rivalues + (i << 1));
+        gsl_power2_fft8192synthesis_scale_skip2 (rivalues + (i << 1), rivalues + (i << 1));
+    }
+  else if (mode == BIG_ANALYSIS)
+    {
+      if (rivalues_in)
+        bitreverse_fft2analysis (n_values, rivalues_in, rivalues);
+      for (i = 0; i < n_values; i += 8192)
+        gsl_power2_fft8192analysis_skip2 (rivalues + (i << 1), rivalues + (i << 1));
     }
   theta *= (double) 1.0 / 8192.;
   last_sin = sin (theta);
@@ -231,7 +244,7 @@ gsl_power2_fftac (const unsigned int n_values,
       case 2048: gsl_power2_fft2048analysis (rivalues_in, rivalues_out);	break;
       case 4096: gsl_power2_fft4096analysis (rivalues_in, rivalues_out);	break;
       case 8192: gsl_power2_fft8192analysis (rivalues_in, rivalues_out);	break;
-      default:	 gsl_power2_fftc_big (n_values, rivalues_in, rivalues_out, +1);
+      default:	 gsl_power2_fftc_big (n_values, rivalues_in, rivalues_out, BIG_ANALYSIS);
     }
 }
 void
@@ -257,7 +270,34 @@ gsl_power2_fftsc (const unsigned int n_values,
       case 2048: gsl_power2_fft2048synthesis (rivalues_in, rivalues_out);	break;
       case 4096: gsl_power2_fft4096synthesis (rivalues_in, rivalues_out);	break;
       case 8192: gsl_power2_fft8192synthesis (rivalues_in, rivalues_out);	break;
-      default:	 gsl_power2_fftc_big (n_values, rivalues_in, rivalues_out, -1);
+      default:	 gsl_power2_fftc_big (n_values, rivalues_in, rivalues_out, BIG_SYNTHESIS);
+    }
+  /* { unsigned int i; for (i = 0; i < n_values << 1; i++) rivalues_out[i] *= (double) n_values; } */
+}
+void
+gsl_power2_fftsc_scale (const unsigned int n_values,
+                        const $IEEE_TYPE  *rivalues_in,
+                        $IEEE_TYPE        *rivalues_out)
+{
+  g_return_if_fail ((n_values & (n_values - 1)) == 0 && n_values >= 1);
+
+  switch (n_values)
+    {
+      case    1: rivalues_out[0] = rivalues_in[0], rivalues_out[1] = rivalues_in[1]; break;
+      case    2: gsl_power2_fft2synthesis_scale (rivalues_in, rivalues_out);	break;
+      case    4: gsl_power2_fft4synthesis_scale (rivalues_in, rivalues_out);	break;
+      case    8: gsl_power2_fft8synthesis_scale (rivalues_in, rivalues_out);	break;
+      case   16: gsl_power2_fft16synthesis_scale (rivalues_in, rivalues_out);	break;
+      case   32: gsl_power2_fft32synthesis_scale (rivalues_in, rivalues_out);	break;
+      case   64: gsl_power2_fft64synthesis_scale (rivalues_in, rivalues_out);	break;
+      case  128: gsl_power2_fft128synthesis_scale (rivalues_in, rivalues_out);	break;
+      case  256: gsl_power2_fft256synthesis_scale (rivalues_in, rivalues_out);	break;
+      case  512: gsl_power2_fft512synthesis_scale (rivalues_in, rivalues_out);	break;
+      case 1024: gsl_power2_fft1024synthesis_scale (rivalues_in, rivalues_out);	break;
+      case 2048: gsl_power2_fft2048synthesis_scale (rivalues_in, rivalues_out);	break;
+      case 4096: gsl_power2_fft4096synthesis_scale (rivalues_in, rivalues_out);	break;
+      case 8192: gsl_power2_fft8192synthesis_scale (rivalues_in, rivalues_out);	break;
+      default:	 gsl_power2_fftc_big (n_values, rivalues_in, rivalues_out, BIG_SYNTHESIS_SCALE);
     }
   /* { unsigned int i; for (i = 0; i < n_values << 1; i++) rivalues_out[i] *= (double) n_values; } */
 }
@@ -280,7 +320,7 @@ gsl_power2_fftar (const unsigned int n_values,
 
   g_return_if_fail ((n_values & (n_values - 1)) == 0 && n_values >= 2);
 
-  gsl_power2_fftac (n_cvalues, r_values_in, rivalues_out);
+  gsl_power2_fftsc (n_cvalues, r_values_in, rivalues_out);
   theta = 3.1415926535897932384626433832795029;
   theta /= (double) n_cvalues;
 
@@ -309,12 +349,12 @@ gsl_power2_fftar (const unsigned int n_values,
       H1im = F2im + F1re;
       H1re = F1im - F2re;
       H2re = F2re - F1im;
-      H2im = H1im - FEim;
+      H2im = FEim - H1im;
       H1re += FEre;
       H2re += FEre;
       H1im += FEim;
       rivalues_out[i] = H1re;
-      rivalues_out[i + 1] = H1im;
+      rivalues_out[i + 1] = -H1im;
       rivalues_out[r] = H2re;
       rivalues_out[r + 1] = H2im;
       WMULTIPLY (Wre, Wim, Dre, Dim);
@@ -322,11 +362,14 @@ gsl_power2_fftar (const unsigned int n_values,
   Dre = rivalues_out[0];
   rivalues_out[0] = Dre + rivalues_out[1];
   rivalues_out[1] = Dre - rivalues_out[1];
+  rivalues_out[n_cvalues + 1] = -rivalues_out[n_cvalues + 1];
 }
-void
-gsl_power2_fftsr (const unsigned int n_values,
-                  const double      *rivalues_in,
-                  double            *r_values_out)
+
+static inline void
+gsl_power2_fftsr_impl (const unsigned int n_values,
+                       const double      *rivalues_in,
+                       double            *r_values_out,
+                       gboolean           need_scale)
 {
   unsigned int n_cvalues = n_values >> 1;
   double Dre, Dim, Wre, Wim, theta, scale;
@@ -360,9 +403,8 @@ gsl_power2_fftsr (const unsigned int n_values,
       ri |= j;
 
       FOre = -FOre;
-      FOim = -FOim;
       FEre *= 0.5;
-      FEim *= 0.5;
+      FEim *= -0.5;
       F2re = FOre * Wim;
       F2im = FOim * Wim;
       F1re = FOre * Wre;
@@ -391,9 +433,16 @@ gsl_power2_fftsr (const unsigned int n_values,
   if (n_values < 4)
     return;
   r_values_out[2] = rivalues_in[i];
-  r_values_out[2 + 1] = rivalues_in[i + 1];
-  scale = n_cvalues;
-  scale = 1.0 / scale;
+  r_values_out[2 + 1] = -rivalues_in[i + 1];
+  if (need_scale)
+    {
+      scale = n_cvalues;
+      scale = 1 / scale;
+    }
+  else
+    {
+      scale = 2;
+    }
   for (i = 0; i < n_values; i += 4)
     BUTTERFLY_10scale (r_values_out[i], r_values_out[i + 1],
                        r_values_out[i + 2], r_values_out[i + 3],
@@ -403,21 +452,39 @@ gsl_power2_fftsr (const unsigned int n_values,
   switch (n_cvalues)
     {
       case    2: break;
-      case    4: gsl_power2_fft4synthesis_skip2 (NULL, r_values_out);	 break;
-      case    8: gsl_power2_fft8synthesis_skip2 (NULL, r_values_out);	 break;
-      case   16: gsl_power2_fft16synthesis_skip2 (NULL, r_values_out);	 break;
-      case   32: gsl_power2_fft32synthesis_skip2 (NULL, r_values_out);	 break;
-      case   64: gsl_power2_fft64synthesis_skip2 (NULL, r_values_out);	 break;
-      case  128: gsl_power2_fft128synthesis_skip2 (NULL, r_values_out);	 break;
-      case  256: gsl_power2_fft256synthesis_skip2 (NULL, r_values_out);	 break;
-      case  512: gsl_power2_fft512synthesis_skip2 (NULL, r_values_out);	 break;
-      case 1024: gsl_power2_fft1024synthesis_skip2 (NULL, r_values_out); break;
-      case 2048: gsl_power2_fft2048synthesis_skip2 (NULL, r_values_out); break;
-      case 4096: gsl_power2_fft4096synthesis_skip2 (NULL, r_values_out); break;
-      case 8192: gsl_power2_fft8192synthesis_skip2 (NULL, r_values_out); break;
-      default:	 gsl_power2_fftc_big (n_cvalues, NULL, r_values_out, -1);
+      case    4: gsl_power2_fft4analysis_skip2 (NULL, r_values_out);	 break;
+      case    8: gsl_power2_fft8analysis_skip2 (NULL, r_values_out);	 break;
+      case   16: gsl_power2_fft16analysis_skip2 (NULL, r_values_out);	 break;
+      case   32: gsl_power2_fft32analysis_skip2 (NULL, r_values_out);	 break;
+      case   64: gsl_power2_fft64analysis_skip2 (NULL, r_values_out);	 break;
+      case  128: gsl_power2_fft128analysis_skip2 (NULL, r_values_out);	 break;
+      case  256: gsl_power2_fft256analysis_skip2 (NULL, r_values_out);	 break;
+      case  512: gsl_power2_fft512analysis_skip2 (NULL, r_values_out);	 break;
+      case 1024: gsl_power2_fft1024analysis_skip2 (NULL, r_values_out); break;
+      case 2048: gsl_power2_fft2048analysis_skip2 (NULL, r_values_out); break;
+      case 4096: gsl_power2_fft4096analysis_skip2 (NULL, r_values_out); break;
+      case 8192: gsl_power2_fft8192analysis_skip2 (NULL, r_values_out); break;
+      default:	 gsl_power2_fftc_big (n_cvalues, NULL, r_values_out, BIG_ANALYSIS);
     }
 }
+
+void
+gsl_power2_fftsr (const unsigned int n_values,
+                  const double      *rivalues_in,
+                  double            *r_values_out)
+{
+  gsl_power2_fftsr_impl (n_values, rivalues_in, r_values_out, FALSE);
+}
+
+void
+gsl_power2_fftsr_scale (const unsigned int n_values,
+                       const double      *rivalues_in,
+                       double            *r_values_out)
+{
+  gsl_power2_fftsr_impl (n_values, rivalues_in, r_values_out, TRUE);
+}
+
+
 void
 gsl_power2_fftar_simple (const unsigned int n_values,
                          const float       *real_values,
@@ -442,10 +509,14 @@ gsl_power2_fftar_simple (const unsigned int n_values,
   complex_values[n_values + 1] = 0.0;
   g_free (rv);
 }
+__EOF
+for FUNC_NAME in fftsr fftsr_scale
+do
+cat << __EOF
 void
-gsl_power2_fftsr_simple (const unsigned int n_values,
-                         const float       *complex_values,
-                         float             *real_values)
+gsl_power2_${FUNC_NAME}_simple (const unsigned int n_values,
+                                const float       *complex_values,
+                                float             *real_values)
 {
   double *cv, *rv;
   guint i;
@@ -458,10 +529,11 @@ gsl_power2_fftsr_simple (const unsigned int n_values,
   while (i--)
     cv[i] = complex_values[i];
   cv[1] = complex_values[n_values];
-  gsl_power2_fftsr (n_values, cv, rv);
+  gsl_power2_$FUNC_NAME (n_values, cv, rv);
   i = n_values;
   while (i--)
     real_values[i] = rv[i];
   g_free (cv);
 }
 __EOF
+done
diff --git a/bse/gsl-fftgen.pl b/bse/gsl-fftgen.pl
index 46a8bd6..ef9a198 100755
--- a/bse/gsl-fftgen.pl
+++ b/bse/gsl-fftgen.pl
@@ -32,7 +32,9 @@ my $min_split = 2048;
 my $min_compress = 512;
 my $fft2loop = 0;
 my $single_stage = 0;
-my $negate_sign = 0;
+my $negate_sign = 1;
+my $func_name = "analysis";
+my $scale = 0;
 my $Wtest = 0;
 
 #
@@ -42,8 +44,9 @@ while ($_ = $ARGV[0], defined $_ && /^-/) {
     shift;
     last if (/^--$/);
     if (/^--fft2loop$/) { $fft2loop = 1 }
-    elsif (/^--analysis$/) { $negate_sign = 0 }
-    elsif (/^--synthesis$/) { $negate_sign = 1 }
+    elsif (/^--analysis$/) { $negate_sign = 1; $func_name = "analysis"; $scale = 0; }
+    elsif (/^--synthesis$/) { $negate_sign = 0; $func_name = "synthesis"; $scale = 0; }
+    elsif (/^--synthesis-scale$/) { $negate_sign = 0; $func_name = "synthesis_scale"; $scale = 1; }
     elsif (/^--Wtest$/) { $Wtest = 1 }
     elsif (/^--double$/) { $ieee_type = "double" }
     elsif (/^--skip-macros$/) { $gen_macros = 0 }
@@ -52,6 +55,7 @@ while ($_ = $ARGV[0], defined $_ && /^-/) {
     elsif (/^--min-compress$/) { $min_compress = shift }
     elsif (/^--single-stage$/) { $single_stage = shift }
 }
+
 # parse arguments
 my @arguments = 0;
 if (defined $ARGV[0]) {
@@ -217,18 +221,18 @@ sub unroll_stage {
 		my $offset2 = $offset1 + $n_points;
 		my $offset1r = bitreverse ($fft_size, $offset1 >> 1) << 1;
 		my $offset2r = bitreverse ($fft_size, $offset2 >> 1) << 1;
-		my $scale = !$negate_sign ? "__1, __0" : sprintf "1.0 / (%s) %u", $tmp_ieee_type, $fft_size;
+		my $scale_args = ! $scale ? "__1, __0" : sprintf "1.0 / (%s) %u", $tmp_ieee_type, $fft_size;
 		printf($indent."BUTTERFLY_10%s (X[%s], X[%s + 1],\n".
 		       $indent."                X[%s], X[%s + 1],\n".
 		       $indent."                Y[%s], Y[%s + 1],\n".
 		       $indent."                Y[%s], Y[%s + 1],\n".
 		       $indent."                %s);\n",
-		       $negate_sign ? "scale" : "",
+		       $scale ? "scale" : "",
 		       $offset1r, $offset1r,
 		       $offset2r, $offset2r,
 		       $offset1, $offset1,
 		       $offset2, $offset2,
-		       $scale);
+		       $scale_args);
 	    }
 	}
     } elsif ($unroll_outer) {
@@ -325,7 +329,7 @@ sub bitreverse_fft2 {
     # mul_result = gsl_complex (c1.re * c2.re - c1.im * c2.im, c1.re * c2.im + c1.im * c2.re);
     printf "
 static inline void
-bitreverse_fft2analysis (const uint         n,
+bitreverse_fft2synthesis (const unsigned int n,
                          const %-6s        *X,
                          %-6s              *Y)
 {
@@ -370,25 +374,26 @@ bitreverse_fft2analysis (const uint         n,
     }
 }
 static inline void
-bitreverse_fft2synthesis (const uint         n,
-                          const %-6s        *X,
-                          %-6s              *Y)
+bitreverse_fft2synthesis_scale (const unsigned int n,
+                               const %-6s        *X,
+                               %-6s              *Y)
 {
   const uint n2 = n >> 1, n1 = n + n2, max = n >> 2;
   uint i, r;
   %s scale = n;
 
   scale = 1.0 / scale;
+
   BUTTERFLY_10scale (X[0], X[1],
 		     X[n], X[n + 1],
 		     Y[0], Y[1],
 		     Y[2], Y[3],
 		     scale);
   BUTTERFLY_10scale (X[n2], X[n2 + 1],
-		     X[n1], X[n1 + 1],
-		     Y[4], Y[5],
-		     Y[6], Y[7],
-		     scale);
+		X[n1], X[n1 + 1],
+		Y[4], Y[5],
+		Y[6], Y[7],
+		scale);
   for (i = 1, r = 0; i < max; i++)
     {
       uint k, j = n >> 1;
@@ -403,20 +408,65 @@ bitreverse_fft2synthesis (const uint         n,
       k = r >> 1;
       j = i << 3;
       BUTTERFLY_10scale (X[k], X[k + 1],
-			 X[k + n], X[k + n + 1],
-			 Y[j], Y[j + 1],
-			 Y[j + 2], Y[j + 3],
-			 scale);
+		         X[k + n], X[k + n + 1],
+		         Y[j], Y[j + 1],
+		         Y[j + 2], Y[j + 3],
+		         scale);
       k += n2;
       j += 4;
       BUTTERFLY_10scale (X[k], X[k + 1],
-			 X[k + n], X[k + n + 1],
-			 Y[j], Y[j + 1],
-			 Y[j + 2], Y[j + 3],
-			 scale);
+		         X[k + n], X[k + n + 1],
+		         Y[j], Y[j + 1],
+		         Y[j + 2], Y[j + 3],
+		         scale);
+    }
+}
+static inline void
+bitreverse_fft2analysis (const unsigned int n,
+                          const %-6s        *X,
+                          %-6s              *Y)
+{
+  const unsigned int n2 = n >> 1, n1 = n + n2, max = n >> 2;
+  unsigned int i, r;
+
+  BUTTERFLY_10 (X[0], X[1],
+		X[n], X[n + 1],
+		Y[0], Y[1],
+		Y[2], Y[3],
+		__0, __1);
+  BUTTERFLY_10 (X[n2], X[n2 + 1],
+		X[n1], X[n1 + 1],
+		Y[4], Y[5],
+		Y[6], Y[7],
+		__0, __1);
+  for (i = 1, r = 0; i < max; i++)
+    {
+      unsigned int k, j = n >> 1;
+
+      while (r >= j)
+	{
+	  r -= j;
+	  j >>= 1;
+	}
+      r |= j;
+
+      k = r >> 1;
+      j = i << 3;
+      BUTTERFLY_10 (X[k], X[k + 1],
+		    X[k + n], X[k + n + 1],
+		    Y[j], Y[j + 1],
+		    Y[j + 2], Y[j + 3],
+		    __0, __1);
+      k += n2;
+      j += 4;
+      BUTTERFLY_10 (X[k], X[k + 1],
+		    X[k + n], X[k + n + 1],
+		    Y[j], Y[j + 1],
+		    Y[j + 2], Y[j + 3],
+		    __0, __1);
     }
 }
-", $ieee_type, $ieee_type, $ieee_type, $ieee_type, $tmp_ieee_type;
+", $ieee_type, $ieee_type, $ieee_type, $ieee_type, $tmp_ieee_type, $ieee_type, $ieee_type, $tmp_ieee_type;
 
     # testing:
     # define BUTTERFLY_10(X1re,X1im,X2re,X2im,Y1re,Y1im,Y2re,Y2im,Wre,Wim,T1re,T1im,T2re,T2im) \
@@ -698,7 +748,7 @@ sub gen_stage {
     if ($n_points == 2) {   # input stage needs special handling
 	if ($kind eq 'L') {
 	    printf "\n%s/* perform fft2 and bitreverse input */\n", $indent;
-	    printf "%sbitreverse_fft2%s (%u, X, Y);\n", $indent, $negate_sign ? "synthesis" : "analysis", $fft_size;
+	    printf "%sbitreverse_fft2%s (%u, X, Y);\n", $indent, $func_name, $fft_size;
 	} elsif ($kind eq 'F') {
 	    printf "\n%s/* perform %u times fft2 */\n", $indent, $times;
 	    unroll_stage ($fft_size, $n_points, $times, 1);
@@ -723,10 +773,10 @@ sub gen_stage {
 	    for (my $i = 0; $i < $times; $i++) {
 		if ($i) {
 		    printf($indent."gsl_power2_fft%u%s_skip2 (X + %u, Y + %u);\n",
-			   $n_points, $negate_sign ? "synthesis" : "analysis", $n_points * $i << 1, $n_points * $i << 1);
+			   $n_points, $func_name, $n_points * $i << 1, $n_points * $i << 1);
 		} else {
 		    printf($indent."gsl_power2_fft%u%s_skip2 (X, Y);\n",
-			   $n_points, $negate_sign ? "synthesis" : "analysis");
+			   $n_points, $func_name);
 		}
 	    }
 	} else {
@@ -786,7 +836,7 @@ print " */\n";
     
     printf "static void\n";
     printf("gsl_power2_fft%u%s%s (const %s *X, %s *Y)\n{\n",
-	   $fft_size, $negate_sign ? "synthesis" : "analysis",
+	   $fft_size, $func_name,
 	   $skip2 ? "_skip2" : "",
 	   $ieee_type, $ieee_type);
     printf "%sregister uint butterfly, block, offset;\n", $indent;
diff --git a/bse/gslfft.h b/bse/gslfft.h
index f65873e..7f9436f 100644
--- a/bse/gslfft.h
+++ b/bse/gslfft.h
@@ -42,15 +42,18 @@ extern "C" {
  * equivalent).
  *
  * In general for the gsl_power2_fft*() family of functions, normalization is
- * only performed during backward transform. However, a popular mathematical
- * strategy of defining the FFT and IFFT in a way that the formulas are
- * symmetric is normalizing both, the forward and backward transform with
- * 1/sqrt(N) - where N is the number of complex values (n_values).
+ * only performed during backward transform if the gsl_power2_fftsc_scale()
+ * is used. No normalization is performed if gsl_power2_fftsc() is used.
+ *
+ * However, a popular mathematical strategy of defining the FFT and IFFT in a
+ * way that the formulas are symmetric is normalizing both, the forward and
+ * backward transform with 1/sqrt(N) - where N is the number of complex values
+ * (n_values).
  *
  * Compared to the above definition, in this implementation, the analyzed
  * values produced by gsl_power2_fftac()/gsl_power2_fftar() will be too large
  * by a factor of sqrt(N), which however are cancelled out on the backward
- * transform.
+ * transform (for _scale variants).
  *
  * Note that the transformation is performed out of place, the input
  * array is not modified, and may not overlap with the output array.
@@ -69,11 +72,10 @@ void	gsl_power2_fftac (const uint         n_values,
  * represents the counterpart to gsl_power2_fftac(), that is, a value
  * array which is transformed into the frequency domain with
  * gsl_power2_fftac() can be reconstructed by issuing gsl_power2_fftsc()
- * on the transform.
+ * on the transform. This function does not perform scaling, so calling
+ * gsl_power2_fftac() and gsl_power2_fftsc() will scale the data with a factor
+ * of n_values.  See also gsl_power2_fftsc_scale().
  *
- * This function also scales the time domain coefficients by a
- * factor of 1.0/n_values which is required for perfect reconstruction
- * of time domain data formerly transformed via gsl_power2_fftac().
  * More details on normalization can be found in the documentation of
  * gsl_power2_fftac().
  *
@@ -85,10 +87,34 @@ void	gsl_power2_fftsc (const uint         n_values,
 			  double            *ri_values_out);
 
 /**
- * @param n_values      Number of real sample values
- * @param r_values_in   Real sample values [0..n_values-1]
- * @param ri_values_out Complex frequency values [0..n_values-1]
+ * gsl_power2_fftsc_scale
+ * @n_values:      Number of complex values
+ * @ri_values_in:  Complex frequency values [0..n_values*2-1]
+ * @ri_values_out: Complex sample values [0..n_values*2-1]
+ * This function performs a decimation in time fourier transformation
+ * in backwards direction with normalization. As such, this function
+ * represents the counterpart to gsl_power2_fftac(), that is, a value
+ * array which is transformed into the frequency domain with
+ * gsl_power2_fftac() can be reconstructed by issuing gsl_power2_fftsc()
+ * on the transform.
  *
+ * This function also scales the time domain coefficients by a
+ * factor of 1.0/n_values which is required for perfect reconstruction
+ * of time domain data formerly transformed via gsl_power2_fftac().
+ * More details on normalization can be found in the documentation of
+ * gsl_power2_fftac().
+ *
+ * Note that the transformation is performed out of place, the input
+ * array is not modified, and may not overlap with the output array.
+ */
+void    gsl_power2_fftsc_scale (const unsigned int n_values,
+			        const double      *ri_values_in,
+			        double            *ri_values_out);
+/**
+ * gsl_power2_fftar
+ * @n_values:      Number of real sample values
+ * @r_values_in:   Real sample values [0..n_values-1]
+ * @ri_values_out: Complex frequency values [0..n_values-1]
  * Real valued variant of gsl_power2_fftac(), the input array contains
  * real valued equidistant sampled data [0..n_values-1], and the output
  * array contains the positive frequency half of the complex valued
@@ -124,6 +150,30 @@ void	gsl_power2_fftar (const uint         n_values,
  * A real valued data set transformed into the frequency domain
  * with gsl_power2_fftar() can be reconstructed using this function.
  *
+ * This function does not perform normalization, so data that is transformed
+ * back from gsl_power2_fftar() will be scaled by a factor of n_values. See
+ * also gsl_power2_fftsr_scale().
+ *
+ * More details on normalization can be found in the documentation of
+ * gsl_power2_fftac().
+ *
+ * Note that the transformation is performed out of place, the input
+ * array is not modified, and may not overlap with the output array.
+ */
+void	gsl_power2_fftsr (const unsigned int n_values,
+			  const double      *ri_values_in,
+			  double            *r_values_out);
+
+/**
+ * gsl_power2_fftsr_scale
+ * @n_values:     Number of real sample values
+ * @ri_values_in: Complex frequency values [0..n_values-1]
+ * @r_values_out: Real sample values [0..n_values-1]
+ * Real valued variant of gsl_power2_fftsc(), counterpart to
+ * gsl_power2_fftar(), using the same frequency storage format.
+ * A real valued data set transformed into the frequency domain
+ * with gsl_power2_fftar() can be reconstructed using this function.
+ *
  * This function also scales the time domain coefficients by a
  * factor of 1.0/(n_values/2) which is required for perfect
  * reconstruction of time domain data formerly transformed via
@@ -134,9 +184,9 @@ void	gsl_power2_fftar (const uint         n_values,
  * Note that the transformation is performed out of place, the input
  * array is not modified, and may not overlap with the output array.
  */
-void	gsl_power2_fftsr (const uint         n_values,
-			  const double      *ri_values_in,
-			  double            *r_values_out);
+void	gsl_power2_fftsr_scale (const unsigned int n_values,
+			        const double      *ri_values_in,
+			        double            *r_values_out);
 
 
 /* --- convenience wrappers --- */
@@ -146,7 +196,10 @@ void	gsl_power2_fftar_simple	(const uint         n_values,
 void	gsl_power2_fftsr_simple	(const uint         n_values,
 				 const float	   *complex_values,
 				 float             *real_values);
-     
+void	gsl_power2_fftsr_scale_simple (const unsigned int n_values,
+				       const float	   *complex_values,
+				       float             *real_values);
+
      
      
 #ifdef __cplusplus



[Date Prev][Date Next]   [Thread Prev][Thread Next]   [Thread Index] [Date Index] [Author Index]