Commit ffd03728 by Reimar H Leike

### added a function that computes the typical signal variance for a given power...

`added a function that computes the typical signal variance for a given power spectrum. Now the wiener filter easy demo uses this function to set its signal level instead of the even worse practice of letting the user specify a signal to nosie ratio`
parent ca2ceecf
Pipeline #28389 passed with stages
in 2 minutes and 49 seconds
 import numpy as np import numpy as np import nifty4 as ift import nifty4 as ift if __name__ == "__main__": if __name__ == "__main__": np.random.seed(43) np.random.seed(43) # Set up physical constants # Set up physical constants ... @@ -12,21 +13,25 @@ if __name__ == "__main__": ... @@ -12,21 +13,25 @@ if __name__ == "__main__": field_variance = 2. field_variance = 2. # Smoothing length of response (in same unit as L) # Smoothing length of response (in same unit as L) response_sigma = 0.01 response_sigma = 0.01 # typical noise amplitude of the measurement noise_level = 1. # Define resolution (pixels per dimension) # Define resolution (pixels per dimension) N_pixels = 256 N_pixels = 256 # Set up derived constants # Set up derived constants k_0 = 1./correlation_length k_0 = 1./correlation_length # Note that field_variance**2 = a*k_0/4. for this analytic form of power #defining a power spectrum with the right correlation length # spectrum #we later set the field variance to the desired value a = field_variance**2/k_0*4. unscaled_pow_spec = (lambda k: 1. / (1 + k/k_0) ** 4) pow_spec = (lambda k: a / (1 + k/k_0) ** 4) pixel_width = L/N_pixels pixel_width = L/N_pixels # Set up the geometry # Set up the geometry s_space = ift.RGSpace([N_pixels, N_pixels], distances=pixel_width) s_space = ift.RGSpace([N_pixels, N_pixels], distances=pixel_width) h_space = s_space.get_default_codomain() h_space = s_space.get_default_codomain() s_var = ift.get_signal_variance(unscaled_pow_spec, h_space) pow_spec = (lambda k: unscaled_pow_spec(k)/s_var*field_variance**2) HT = ift.HarmonicTransformOperator(h_space, s_space) HT = ift.HarmonicTransformOperator(h_space, s_space) # Create mock data # Create mock data ... @@ -36,11 +41,8 @@ if __name__ == "__main__": ... @@ -36,11 +41,8 @@ if __name__ == "__main__": R = HT*ift.create_harmonic_smoothing_operator((h_space,), 0, R = HT*ift.create_harmonic_smoothing_operator((h_space,), 0, response_sigma) response_sigma) noiseless_data = R(sh) noiseless_data = R(sh) signal_to_noise = 1. N = ift.ScalingOperator(noise_level**2, s_space) noise_amplitude = noiseless_data.val.std()/signal_to_noise N = ift.ScalingOperator(noise_amplitude**2, s_space) n = N.draw_sample() n = N.draw_sample() d = noiseless_data + n d = noiseless_data + n ... ...
 ... @@ -28,7 +28,9 @@ from .logger import logger ... @@ -28,7 +28,9 @@ from .logger import logger __all__ = ['PS_field', __all__ = ['PS_field', 'power_analyze', 'power_analyze', 'create_power_operator', 'create_power_operator', 'create_harmonic_smoothing_operator'] 'create_harmonic_smoothing_operator', 'get_signal_variance'] def PS_field(pspace, func): def PS_field(pspace, func): ... @@ -38,6 +40,33 @@ def PS_field(pspace, func): ... @@ -38,6 +40,33 @@ def PS_field(pspace, func): return Field(pspace, val=data) return Field(pspace, val=data) def get_signal_variance(spec, space): """ Computes how much a field with a given power spectrum will vary in space This is a small helper function that computes how the expected variance of a harmonically transformed sample of this power spectrum. Parameters --------- spec: method a method that takes one k-value and returns the power spectrum at that location space: PowerSpace or any harmonic Domain If this function is given a harmonic domain, it creates the naturally binned PowerSpace to that domain. The field, for which the signal variance is then computed, is assumed to have this PowerSpace as naturally binned PowerSpace """ if space.harmonic: space = PowerSpace(space) if not isinstance(space, PowerSpace): raise ValueError("space must be either a harmonic space or Power space.") field = PS_field(space, spec) dist = PowerDistributor(space.harmonic_partner, space) k_field = dist(field) return k_field.weight(2).sum() def _single_power_analyze(field, idx, binbounds): def _single_power_analyze(field, idx, binbounds): power_domain = PowerSpace(field.domain[idx], binbounds) power_domain = PowerSpace(field.domain[idx], binbounds) pd = PowerDistributor(field.domain, power_domain, idx) pd = PowerDistributor(field.domain, power_domain, idx) ... ...
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!