progress

demos/wiener_filter_via_curvature.py

@@ -17,7 +17,7 @@ if __name__ == "__main__":
     # smoothing length of response
     response_sigma = 0.01*nu.m
     # The signal to noise ratio
-    signal_to_noise = 0.7
+    signal_to_noise = 70.7
 
     # note that field_variance**2 = a*k_0/4. for this analytic form of power
     # spectrum
@@ -49,27 +49,31 @@ if __name__ == "__main__":
     mock_harmonic = ift.power_synthesize(mock_power, real_signal=True)
     print mock_harmonic.val[0]/nu.K/(nu.m**dimensionality)
     mock_signal = fft(mock_harmonic)
-    print mock_signal.val[0]/nu.K
+    print "msig",mock_signal.val[0:10]/nu.K
 
-    exposure = 1./nu.K
-    R = ift.ResponseOperator(signal_space, sigma=(response_sigma,),
-                             exposure=(exposure,))
+    sensitivity = (1./nu.m)**dimensionality/nu.K
+    R = ift.ResponseOperator(signal_space, sigma=(0.*response_sigma,),
+                             sensitivity=(sensitivity,))
     data_domain = R.target[0]
     R_harmonic = R*fft
 
+    noise_amplitude = 1./signal_to_noise*field_sigma*sensitivity*((L/N_pixels)**dimensionality)
+    print noise_amplitude
     N = ift.DiagonalOperator(
-        ift.Field.full(data_domain,
-                       mock_signal.var()/signal_to_noise))
+        ift.Field.full(data_domain, noise_amplitude**2))
     noise = ift.Field.from_random(
         domain=data_domain, random_type='normal',
-        std=mock_signal.std()/np.sqrt(signal_to_noise), mean=0)
-    data = R(mock_signal) #+ noise
-    print data.val[5]
+        std=noise_amplitude, mean=0)
+    data = R(mock_signal)
+    print data.val[5:10]
+    data += noise
+    print data.val[5:10]
      # Wiener filter
 
     j = R_harmonic.adjoint_times(N.inverse_times(data))
+    print "xx",j.val[0]*nu.K*(nu.m**dimensionality)
     ctrl = ift.GradientNormController(
-        verbose=True, tol_abs_gradnorm=1e-4/nu.K)
+        verbose=True, tol_abs_gradnorm=1e-40/(nu.K*(nu.m**dimensionality)))
     inverter = ift.ConjugateGradient(controller=ctrl)
     wiener_curvature = ift.library.WienerFilterCurvature(
         S=S, N=N, R=R_harmonic, inverter=inverter)

nifty/operators/response_operator.py

@@ -32,33 +32,36 @@ class GeometryRemover(LinearOperator):
         return Field(self._domain, val=x.val).weight(1)
 
 
-def ResponseOperator(domain, sigma, exposure):
+def ResponseOperator(domain, sigma, sensitivity):
+    # sensitivity has units 1/field/volume and gives a measure of how much
+    # the instrument will excited when it is exposed to a certain field
+    # volume amplitude
     domain = DomainTuple.make(domain)
-    ncomp = len(exposure)
+    ncomp = len(sensitivity)
     if len(sigma) != ncomp or len(domain) != ncomp:
-        raise ValueError("length mismatch between sigma, exposure "
+        raise ValueError("length mismatch between sigma, sensitivity "
                          "and domain")
     ncomp = len(sigma)
     if ncomp == 0:
         raise ValueError("Empty response operator not allowed")
 
     kernel = None
-    expo = None
+    sensi = None
     for i in range(ncomp):
         if sigma[i] > 0:
             op = FFTSmoothingOperator(domain, sigma[i], space=i)
             kernel = op if kernel is None else op*kernel
-        if np.isscalar(exposure[i]):
-            if exposure[i] != 1.:
-                op = ScalingOperator(exposure[i], domain)
-                expo = op if expo is None else op*expo
-        elif isinstance(exposure[i], Field):
-            op = DiagonalOperator(exposure[i], domain=domain, spaces=i)
-            expo = op if expo is None else op*expo
+        if np.isscalar(sensitivity[i]):
+            if sensitivity[i] != 1.:
+                op = ScalingOperator(sensitivity[i], domain)
+                sensi = op if sensi is None else op*sensi
+        elif isinstance(sensitivity[i], Field):
+            op = DiagonalOperator(sensitivity[i], domain=domain, spaces=i)
+            sensi = op if sensi is None else op*sensi
 
     res = GeometryRemover(domain)
-    if expo is not None:
-        res = res * expo
+    if sensi is not None:
+        res = res * sensi
     if kernel is not None:
         res = res * kernel
     return res