Corrected power_synthesize issues.

nifty/field.py

@@ -348,10 +348,15 @@ class Field(Loggable, Versionable, object):
         if real_signal:
             result_val_list = [harmonic_domain.hermitian_decomposition(
                                     x.val,
-                                    axes=x.domain_axes[power_space_index])[0]
+                                    axes=x.domain_axes[power_space_index],
+                                    preserve_gaussian_variance=True)[0]
                                for x in result_list]
         else:
-            result_val_list = [x.val for x in result_list]
+#            # if the synthesized field is complex in signal space,
+#            # one must correct the variance here, since one draws
+#            # sqrt(twice) the power via real- and imaginary-part
+#            result_val_list = [x.val*np.sqrt(0.5) for x in result_list]
+             result_val_list = [x.val for x in result_list]
 
         # weight the random fields with the power spectrum
         # therefore get the pindex from the power space

nifty/spaces/lm_space/lm_space.py

@@ -111,7 +111,8 @@ class LMSpace(Space):
         super(LMSpace, self).__init__(dtype)
         self._lmax = self._parse_lmax(lmax)
 
-    def hermitian_decomposition(self, x, axes=None):
+    def hermitian_decomposition(self, x, axes=None,
+                                preserve_gaussian_variance=False):
         hermitian_part = x.copy_empty()
         anti_hermitian_part = x.copy_empty()
         hermitian_part[:] = x.real

nifty/spaces/rg_space/rg_space.py

@@ -150,7 +150,8 @@ class RGSpace(Space):
         self._distances = self._parse_distances(distances)
         self._zerocenter = self._parse_zerocenter(zerocenter)
 
-    def hermitian_decomposition(self, x, axes=None):
+    def hermitian_decomposition(self, x, axes=None,
+                                preserve_gaussian_variance=False):
         # compute the hermitian part
         flipped_x = self._hermitianize_inverter(x, axes=axes)
         flipped_x = flipped_x.conjugate()
@@ -160,8 +161,39 @@ class RGSpace(Space):
 
         # use subtraction since it is faster than flipping another time
         anti_hermitian_part = (x-hermitian_part)/1j
+
+        if preserve_gaussian_variance:
+            hermitian_part, anti_hermitian_part = \
+                self._hermitianize_correct_variance(hermitian_part,
+                                                    anti_hermitian_part,
+                                                    axes=axes)
+
         return (hermitian_part, anti_hermitian_part)
 
+    def _hermitianize_correct_variance(self, hermitian_part,
+                                       anti_hermitian_part, axes):
+        # Correct the variance by multiplying sqrt(2)
+        hermitian_part = hermitian_part * np.sqrt(2)
+        anti_hermitian_part = anti_hermitian_part * np.sqrt(2)
+
+        # The fixed points of the point inversion must not be avaraged.
+        # Hence one must divide out the sqrt(2) again
+        # -> Get the middle index of the array
+        mid_index = np.array(hermitian_part.shape, dtype=np.int) // 2
+        dimensions = mid_index.size
+        # Use ndindex to iterate over all combinations of zeros and the
+        # mid_index in order to correct all fixed points.
+        if axes is None:
+            axes = xrange(dimensions)
+
+        ndlist = [2 if i in axes else 1 for i in xrange(dimensions)]
+        ndlist = tuple(ndlist)
+        for i in np.ndindex(ndlist):
+            temp_index = tuple(i * mid_index)
+            hermitian_part[temp_index] /= np.sqrt(2)
+            anti_hermitian_part[temp_index] /= np.sqrt(2)
+        return hermitian_part, anti_hermitian_part
+
     def _hermitianize_inverter(self, x, axes):
         # calculate the number of dimensions the input array has
         dimensions = len(x.shape)

nifty/spaces/space/space.py

@@ -229,7 +229,8 @@ class Space(DomainObject):
         raise NotImplementedError(
             "There is no generic co-smoothing kernel for Space base class.")
 
-    def hermitian_decomposition(self, x, axes=None):
+    def hermitian_decomposition(self, x, axes=None,
+                                preserve_gaussian_variance=False):
         raise NotImplementedError
 
     def __repr__(self):