diff --git a/demos/paper_demos/cartesian_wiener_filter.py b/demos/paper_demos/cartesian_wiener_filter.py
index 030ffb320d3a9703f7071ce96179c5c0166ed261..655d7cf015caab288107dbb4427496064be342b9 100644
--- a/demos/paper_demos/cartesian_wiener_filter.py
+++ b/demos/paper_demos/cartesian_wiener_filter.py
@@ -52,7 +52,7 @@ if __name__ == "__main__":
             np.outer(ift.dobj.to_global_data(mock_power_1.val),
                      ift.dobj.to_global_data(mock_power_2.val))))
 
-    diagonal = ift.power_synthesize_special(mock_power, spaces=(0, 1))**2
+    diagonal = ift.power_synthesize_nonrandom(mock_power, spaces=(0, 1))**2
 
     S = ift.DiagonalOperator(diagonal)
 
diff --git a/demos/wiener_filter_easy.py b/demos/wiener_filter_easy.py
index 7fe7fb08c4780a7aa286a58323a0be4cb9bad26b..ab9d181965636f0f8bb96435ff59cfb329fe2473 100644
--- a/demos/wiener_filter_easy.py
+++ b/demos/wiener_filter_easy.py
@@ -13,7 +13,7 @@ class PropagatorOperator(ift.InversionEnabler, ift.EndomorphicOperator):
         self.R = R
         self.N = N
         self.Sh = Sh
-        self.fft = ift.FFTOperator(R.domain, target=Sh.domain)
+        self.fft = ift.FFTOperator(R.domain, target=Sh.domain[0])
 
     def _inverse_times(self, x):
         return self.R.adjoint_times(self.N.inverse_times(self.R(x))) \
diff --git a/nifty/operators/fft_operator.py b/nifty/operators/fft_operator.py
index 2fac5f1b249c9a624659afe86e58fd7c40a88569..0e50df1269655ab26093471813d3f472f05487ba 100644
--- a/nifty/operators/fft_operator.py
+++ b/nifty/operators/fft_operator.py
@@ -72,7 +72,6 @@ class FFTOperator(LinearOperator):
         if "domain" or "target" are not of the proper type.
     """
 
-    # MR FIXME: target should only be a single DomainObject, not the full tuple
     def __init__(self, domain, target=None, space=None):
         super(FFTOperator, self).__init__()
 
@@ -89,13 +88,13 @@ class FFTOperator(LinearOperator):
 
         adom = self.domain[self._space]
         if target is None:
-            target = [dom for dom in self.domain]
-            target[self._space] = adom.get_default_codomain()
+            target = adom.get_default_codomain()
 
-        self._target = DomainTuple.make(target)
-        atgt = self._target[space]
-        adom.check_codomain(atgt)
-        atgt.check_codomain(adom)
+        self._target = [dom for dom in self.domain]
+        self._target[self._space] = target
+        self._target = DomainTuple.make(self._target)
+        adom.check_codomain(target)
+        target.check_codomain(adom)
 
         if self._target[space].harmonic:
             pdom, hdom = (self._domain, self._target)
diff --git a/nifty/sugar.py b/nifty/sugar.py
index 6e99c073d516bbfea9cfa76b912635c672ba5aea..98613ec77b99c79f1f728843552591d46606de89 100644
--- a/nifty/sugar.py
+++ b/nifty/sugar.py
@@ -24,7 +24,7 @@ from . import Space, PowerSpace, Field, ComposedOperator, DiagonalOperator,\
 __all__ = ['PS_field',
            'power_analyze',
            'power_synthesize',
-           'power_synthesize_special',
+           'power_synthesize_nonrandom',
            'create_power_field',
            'create_power_operator',
            'create_composed_fft_operator']
@@ -58,7 +58,7 @@ def power_analyze(field, spaces=None, binbounds=None,
     field : Field
         The field to be analyzed
     spaces : int *optional*
-        The subspace for which the powerspectrum shall be computed.
+        The set of subspaces for which the powerspectrum shall be computed.
         (default : None).
     binbounds : array-like *optional*
         Inner bounds of the bins (default : None).
@@ -116,7 +116,7 @@ def power_analyze(field, spaces=None, binbounds=None,
     return parts[0] + 1j*parts[1] if keep_phase_information else parts[0]
 
 
-def _compute_spec(field, spaces):
+def power_synthesize_nonrandom(field, spaces=None):
     spaces = range(len(field.domain)) if spaces is None \
              else utilities.cast_iseq_to_tuple(spaces)
 
@@ -168,7 +168,7 @@ def power_synthesize(field, spaces=None, real_power=True, real_signal=True):
     ValueError : If domain specified by `spaces` is not a PowerSpace.
     """
 
-    spec = _compute_spec(field, spaces)
+    spec = power_synthesize_nonrandom(field, spaces)
     self_real = not np.issubdtype(spec.dtype, np.complexfloating)
     if (not real_power) and self_real:
         raise ValueError("can't draw complex realizations from real spectrum")
@@ -181,11 +181,6 @@ def power_synthesize(field, spaces=None, real_power=True, real_signal=True):
                                 else np.complex)
               for x in range(1 if real_power else 2)]
 
-    # MR FIXME: dummy call - will be removed soon
-    if real_signal:
-        field.from_random('normal', mean=0., std=1.,
-                          domain=spec.domain, dtype=np.float)
-
     result[0] *= spec if self_real else spec.real
     if not real_power:
         result[1] *= spec.imag
@@ -193,16 +188,6 @@ def power_synthesize(field, spaces=None, real_power=True, real_signal=True):
     return result[0] if real_power else result[0] + 1j*result[1]
 
 
-def power_synthesize_special(field, spaces=None):
-    spec = _compute_spec(field, spaces)
-
-    # MR FIXME: dummy call - will be removed soon
-    field.from_random('normal', mean=0., std=1.,
-                      domain=spec.domain, dtype=np.complex)
-
-    return spec
-
-
 def create_power_field(domain, power_spectrum, dtype=None):
     if not callable(power_spectrum):  # we have a Field living on a PowerSpace
         if not isinstance(power_spectrum, Field):
@@ -241,7 +226,6 @@ def create_power_operator(domain, power_spectrum, space=None, dtype=None):
     Returns
     -------
     DiagonalOperator : An operator that implements the given power spectrum.
-
     """
     domain = DomainTuple.make(domain)
     if space is None:
@@ -261,7 +245,7 @@ def create_composed_fft_operator(domain, codomain=None, all_to='other'):
 
     if codomain is None:
         codomain = [None]*len(domain)
-    interdomain = list(domain)
+    tdom = list(domain)
     for i, space in enumerate(domain):
         if not isinstance(space, Space):
             continue
@@ -269,11 +253,9 @@ def create_composed_fft_operator(domain, codomain=None, all_to='other'):
                 (all_to == 'position' and space.harmonic) or
                 (all_to == 'harmonic' and not space.harmonic)):
             if codomain[i] is None:
-                interdomain[i] = domain[i].get_default_codomain()
+                tgt = domain[i].get_default_codomain()
             else:
-                interdomain[i] = codomain[i]
-            fft_op_list += [FFTOperator(domain=domain, target=interdomain,
-                                        space=i)]
-        # MR FIXME: this looks slightly fishy...
-        domain = interdomain
+                tgt = codomain[i]
+            fft_op_list += [FFTOperator(domain=tdom, target=tgt, space=i)]
+            tdom[i] = tgt
     return ComposedOperator(fft_op_list)