diff --git a/test/noise_realisation_tests.py b/test/noise_realisation_tests.py
index 8f96cb30c69383527ae77cdee4c765ed5d95d190..d69d4921bcb9b0cb301ceea5edcc8c5afee49254 100644
--- a/test/noise_realisation_tests.py
+++ b/test/noise_realisation_tests.py
@@ -8,7 +8,7 @@ class TestNoiseRealisation(unittest.TestCase):
def test_averaged_noise(self):
time_duration = 1.
sampling_frequency = 4096.
- factor = np.sqrt(2./time_duration)
+ factor = np.sqrt(2 / time_duration)
n_avg = 1000
psd_avg = 0
interferometer = tupak.gw.detector.get_empty_interferometer('H1')
@@ -16,11 +16,12 @@ class TestNoiseRealisation(unittest.TestCase):
interferometer.set_data(sampling_frequency, time_duration, from_power_spectral_density=True)
psd_avg += abs(interferometer.data)**2
- psd_avg = psd_avg/n_avg
- asd_avg = np.sqrt(abs(psd_avg))
+ psd_avg = psd_avg / n_avg
+ asd_avg = np.sqrt(abs(psd_avg)) * interferometer.frequency_mask
- a = interferometer.amplitude_spectral_density_array/factor
+ a = np.nan_to_num(interferometer.amplitude_spectral_density_array / factor * interferometer.frequency_mask)
b = asd_avg
+ print(a, b)
self.assertTrue(np.allclose(a, b, rtol=1e-1))
def test_noise_normalisation(self):
diff --git a/tupak/gw/detector.py b/tupak/gw/detector.py
index 34a868a30f547236b54c126bc832e17e8a05b289..bf7d8bd380fbf3866b534afce04d2279da36f5ec 100644
--- a/tupak/gw/detector.py
+++ b/tupak/gw/detector.py
@@ -574,7 +574,7 @@ class Interferometer(object):
header='f h(f)')
-class PowerSpectralDensity:
+class PowerSpectralDensity(object):
def __init__(self, asd_file=None, psd_file='aLIGO_ZERO_DET_high_P_psd.txt', frame_file=None, asd_array=None,
psd_array=None, frequencies=None, epoch=0,
@@ -622,10 +622,6 @@ class PowerSpectralDensity:
Name of the ASD file
frequencies: array_like
Array containing the frequencies of the ASD/PSD values
- frequency_noise_realization: list
- TODO: This isn't doing anything right now
- interpolated_frequency: list
- TODO: This isn't doing anything right now
power_spectral_density: array_like
Array representation of the PSD
power_spectral_density_file: str
@@ -633,12 +629,10 @@ class PowerSpectralDensity:
power_spectral_density_interpolated: scipy.interpolated.interp1d
Interpolated function of the PSD
"""
+ self.__power_spectral_density = None
+ self.__amplitude_spectral_density = None
self.frequencies = []
- self.power_spectral_density = []
- self.amplitude_spectral_density = []
- self.frequency_noise_realization = []
- self.interpolated_frequency = []
self.power_spectral_density_interpolated = None
if asd_file is not None:
@@ -666,13 +660,12 @@ class PowerSpectralDensity:
psd = strain.psd(fftlength=fft_length, window=window)
self.frequencies = psd.frequencies
self.power_spectral_density = psd.value
- self.amplitude_spectral_density = self.power_spectral_density ** 0.5
- self.interpolate_power_spectral_density()
elif frequencies is not None:
+ self.frequencies = frequencies
if asd_array is not None:
- self._set_from_amplitude_spectral_density(frequencies, asd_array)
+ self.amplitude_spectral_density = asd_array
elif psd_array is not None:
- self._set_from_amplitude_spectral_density(frequencies, psd_array)
+ self.power_spectral_density = psd_array
else:
if psd_file is None:
logging.info("No power spectral density provided, using aLIGO, zero detuning, high power.")
@@ -685,6 +678,28 @@ class PowerSpectralDensity:
min(self.power_spectral_density)))
logging.warning("You may have intended to provide this as an amplitude spectral density.")
+ @property
+ def power_spectral_density(self):
+ return self.__power_spectral_density
+
+ @power_spectral_density.setter
+ def power_spectral_density(self, power_spectral_density):
+ self._check_frequency_array_matches_density_array(power_spectral_density)
+ self.__power_spectral_density = power_spectral_density
+ self._interpolate_power_spectral_density()
+ self.__amplitude_spectral_density = power_spectral_density**0.5
+
+ @property
+ def amplitude_spectral_density(self):
+ return self.__amplitude_spectral_density
+
+ @amplitude_spectral_density.setter
+ def amplitude_spectral_density(self, amplitude_spectral_density):
+ self._check_frequency_array_matches_density_array(amplitude_spectral_density)
+ self.__amplitude_spectral_density = amplitude_spectral_density
+ self.__power_spectral_density = amplitude_spectral_density**2
+ self._interpolate_power_spectral_density()
+
def import_amplitude_spectral_density(self):
"""
Automagically load one of the amplitude spectral density curves contained in the noise_curves directory.
@@ -695,8 +710,7 @@ class PowerSpectralDensity:
if '/' not in self.amplitude_spectral_density_file:
self.amplitude_spectral_density_file = os.path.join(os.path.dirname(__file__), 'noise_curves',
self.amplitude_spectral_density_file)
- spectral_density = np.genfromtxt(self.amplitude_spectral_density_file)
- self._set_from_amplitude_spectral_density(spectral_density[:, 0], spectral_density[:, 1])
+ self.frequencies, self.amplitude_spectral_density = np.genfromtxt(self.amplitude_spectral_density_file).T
def import_power_spectral_density(self):
"""
@@ -708,26 +722,20 @@ class PowerSpectralDensity:
if '/' not in self.power_spectral_density_file:
self.power_spectral_density_file = os.path.join(os.path.dirname(__file__), 'noise_curves',
self.power_spectral_density_file)
- spectral_density = np.genfromtxt(self.power_spectral_density_file)
- self._set_from_power_spectral_density(spectral_density[:, 0], spectral_density[:, 1])
-
- def _set_from_amplitude_spectral_density(self, frequencies, amplitude_spectral_density):
- self.frequencies = frequencies
- self.amplitude_spectral_density = amplitude_spectral_density
- self.power_spectral_density = self.amplitude_spectral_density ** 2
- self.interpolate_power_spectral_density()
+ self.frequencies, self.power_spectral_density = np.genfromtxt(self.power_spectral_density_file).T
- def _set_from_power_spectral_density(self, frequencies, power_spectral_density):
- self.frequencies = frequencies
- self.power_spectral_density = power_spectral_density
- self.amplitude_spectral_density = self.power_spectral_density ** 0.5
- self.interpolate_power_spectral_density()
+ def _check_frequency_array_matches_density_array(self, density_array):
+ """Check if the provided frequency and spectral density arrays match."""
+ try:
+ self.frequencies - density_array
+ except ValueError as e:
+ raise(e, 'Provided spectral density does not match frequency array. Not updating.')
- def interpolate_power_spectral_density(self):
+ def _interpolate_power_spectral_density(self):
"""Interpolate the loaded PSD so it can be resampled for arbitrary frequency arrays."""
self.power_spectral_density_interpolated = interp1d(self.frequencies, self.power_spectral_density,
bounds_error=False,
- fill_value=max(self.power_spectral_density))
+ fill_value=np.inf)
def get_noise_realisation(self, sampling_frequency, duration):
"""
@@ -749,6 +757,8 @@ class PowerSpectralDensity:
white_noise, frequencies = utils.create_white_noise(sampling_frequency, duration)
interpolated_power_spectral_density = self.power_spectral_density_interpolated(frequencies)
frequency_domain_strain = interpolated_power_spectral_density ** 0.5 * white_noise
+ out_of_bounds = (frequencies < min(self.frequencies)) | (frequencies > max(self.frequencies))
+ frequency_domain_strain[out_of_bounds] = 0 * (1 + 1j)
return frequency_domain_strain, frequencies