Skip to content

Commit

Permalink
Merge 2d102e5 into 9fe0576
Browse files Browse the repository at this point in the history
  • Loading branch information
@nkern
nkern committed May 12, 2019
2 parents 9fe0576 + 2d102e5 commit 8c6f8f5
Show file tree
Hide file tree
Showing 4 changed files with 41 additions and 48 deletions.
40 changes: 22 additions & 18 deletions hera_pspec/noise.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -6,9 +6,10 @@
from collections import OrderedDict as odict


def calc_P_N(scalar, Tsys, t_int, Ncoherent=1, Nincoherent=None, form='Pk', k=None):
def calc_P_N(scalar, Tsys, t_int, Ncoherent=1, Nincoherent=None, form='Pk', k=None, component='real'):
"""
Calculate the noise power spectrum via Eqn. (21) of Cheng et al. 2018
Calculate the noise power spectrum via Eqn. (22) of Cheng et al. 2018 for a specified
component of the power spectrum.
The noise power spectrum is written as
Expand All @@ -17,24 +18,22 @@ def calc_P_N(scalar, Tsys, t_int, Ncoherent=1, Nincoherent=None, form='Pk', k=No
where scalar is a nomalization given by the cosmological model and beam response, i.e. X2Y * Omega_eff
Tsys is the system temp in Kelvin, t_int is the integration time of the underlying data [sec],
Ncoherent is the number of coherent averages before forming power spectra, and Nincoherent is the
number of incoherent averages after squaring.
number of incoherent averages after squaring. If component is 'real' or 'imag' an additional factor
of 1/sqrt(2) is multiplied.
For an ensemble of power spectra with the same Tsys, this estimate should equal their RMS value.
Parameters
----------
scalar : float, Power spectrum normalization factor: X2Y(z) * Omega_P^2 / Omega_PP
Tsys : float, System temperature in Kelvin
t_int : float, integration time of power spectra in seconds
Ncoherent : int, number of coherent averages of visibility data with integration time t_int
Total integration time is t_int * Ncoherent
Nincoherent : int, number of incoherent averages of pspectra (i.e. after squaring).
form : str, power spectra form 'Pk' for P(k) and 'DelSq' for Delta^2(k)
k : float ndarray, cosmological wave-vectors in h Mpc^-1, only needed if form == 'DelSq'
component : str, options=['real', 'imag', 'abs']
If component is real or imag, divide by an extra factor of sqrt(2)
Returns (P_N)
-------
Expand All @@ -43,6 +42,7 @@ def calc_P_N(scalar, Tsys, t_int, Ncoherent=1, Nincoherent=None, form='Pk', k=No
"""
# assert form
assert form in ('Pk', 'DelSq'), "form must be either 'Pk' or 'DelSq' for P(k) or Delta^2(k) respectively"
assert component in ['abs', 'real', 'imag'], "component must be one of 'real', 'imag', 'abs'"

# convert to mK
Tsys *= 1e3
Expand All @@ -57,6 +57,10 @@ def calc_P_N(scalar, Tsys, t_int, Ncoherent=1, Nincoherent=None, form='Pk', k=No
if Nincoherent is not None:
P_N /= np.sqrt(Nincoherent)

# parse component
if component in ['real', 'imag']:
P_N /= np.sqrt(2)

# Convert to Delta Sq
if form == 'DelSq':
assert k is not None, "if form == 'DelSq' then k must be fed"
Expand Down Expand Up @@ -160,9 +164,10 @@ def calc_scalar(self, freqs, pol, num_steps=5000, little_h=True):
self.subband = freqs
self.pol = pol

def calc_P_N(self, Tsys, t_int, Ncoherent=1, Nincoherent=None, form='Pk', k=None):
def calc_P_N(self, Tsys, t_int, Ncoherent=1, Nincoherent=None, form='Pk', k=None, component='real'):
"""
Calculate the noise power spectrum via Eqn. (21) of Cheng et al. 2018
Calculate the noise power spectrum via Eqn. (22) of Cheng et al. 2018 for a specified
component of the power spectrum.
The noise power spectrum is written as
Expand All @@ -171,22 +176,21 @@ def calc_P_N(self, Tsys, t_int, Ncoherent=1, Nincoherent=None, form='Pk', k=None
where scalar is a nomalization given by the cosmological model and beam response, i.e. X2Y * Omega_eff
Tsys is the system temp in Kelvin, t_int is the integration time of the underlying data [sec],
Ncoherent is the number of coherent averages before forming power spectra, and Nincoherent is the
number of incoherent averages after squaring.
number of incoherent averages after squaring. If component is 'real' or 'imag' a factor of 1/sqrt(2)
is multiplied.
Parameters
----------
scalar : float, Power spectrum normalization factor: X2Y(z) * Omega_P^2 / Omega_PP
Tsys : float, System temperature in Kelvin
t_int : float, integration time of power spectra in seconds
Ncoherent : int, number of coherent averages of visibility data with integration time t_int
Total integration time is t_int * Ncoherent
Nincoherent : int, number of incoherent averages of pspectra (i.e. after squaring).
form : str, power spectra form 'Pk' for P(k) and 'DelSq' for Delta^2(k)
k : float ndarray, cosmological wave-vectors in h Mpc^-1, only needed if form == 'DelSq'
component : str, options=['real', 'imag', 'abs']
If component is real or imag, divide by an extra factor of sqrt(2)
Returns (P_N)
-------
Expand All @@ -201,7 +205,7 @@ def calc_P_N(self, Tsys, t_int, Ncoherent=1, Nincoherent=None, form='Pk', k=None

# calculate P_N
P_N = calc_P_N(self.scalar, Tsys, t_int, Ncoherent=Ncoherent, Nincoherent=Nincoherent, form=form,
k=k)
k=k, component=component)

return P_N

Expand Down
23 changes: 8 additions & 15 deletions hera_pspec/tests/test_noise.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -67,7 +67,7 @@ def test_calc_P_N(self):
P_N = self.sense.calc_P_N(Tsys, t_int, Ncoherent=1, Nincoherent=1,
form='Pk')
nt.assert_true(isinstance(P_N, (float, np.float)))
nt.assert_true(np.isclose(P_N, 908472312787.53491))
nt.assert_true(np.isclose(P_N, 642386932892.2921))
# calculate DelSq
Dsq = self.sense.calc_P_N(Tsys, t_int, k=k, Ncoherent=1,
Nincoherent=1, form='DelSq')
Expand Down Expand Up @@ -110,22 +110,15 @@ def test_noise_validation():

# get noise spectra from one of the blpairs
P_N = list(uvp.generate_noise_spectra(0, ('xx','xx'), Tsys,
blpairs=uvp.get_blpairs()[:1],
num_steps=2000).values())[0][0, 0]
blpairs=uvp.get_blpairs()[:1], num_steps=2000,
component='real').values())[0][0, 0]

# get P_std of real spectra for each baseline across time axis
P_stds = np.array([np.std(uvp.get_data((0, bl, ('xx','xx'))).real, axis=1)
for bl in uvp.get_blpairs()])
# get P_rms of real spectra for each baseline across time axis
Pspec = np.array([uvp.get_data((0, bl, ('xx', 'xx'))).real for bl in uvp.get_blpairs()])
P_rms = np.sqrt(np.mean(np.abs(Pspec)**2))

# get average P_std_avg and its standard error
P_std_avg = np.mean(P_stds)

# assert close to P_N: 2%
# This should be updated to be within standard error on P_std_avg
# This should be updated to be within standard error on P_rms
# when the spw_range-variable pspec amplitude bug is resolved
nt.assert_true(np.abs(P_std_avg - P_N) / P_N < 0.02)




nt.assert_true(np.abs(P_rms - P_N) / P_N < 0.02)

10 changes: 5 additions & 5 deletions hera_pspec/tests/test_uvpspec.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -335,24 +335,24 @@ def test_sense(self):
uvp = copy.deepcopy(self.uvp)

# test generate noise spectra
P_N = uvp.generate_noise_spectra(0, 1515, 500, form='Pk', real=True)
P_N = uvp.generate_noise_spectra(0, 1515, 500, form='Pk', component='real')
nt.assert_equal(P_N[101102101102].shape, (10, 30))

# test smaller system temp
P_N2 = uvp.generate_noise_spectra(0, 1515, 400, form='Pk', real=True)
P_N2 = uvp.generate_noise_spectra(0, 1515, 400, form='Pk', component='real')
nt.assert_true((P_N[101102101102] > P_N2[101102101102]).all())

# test complex
P_N2 = uvp.generate_noise_spectra(0, 1515, 500, form='Pk', real=False)
P_N2 = uvp.generate_noise_spectra(0, 1515, 500, form='Pk', component='abs')
nt.assert_true((P_N[101102101102] < P_N2[101102101102]).all())

# test Dsq
Dsq = uvp.generate_noise_spectra(0, 1515, 500, form='DelSq', real=True)
Dsq = uvp.generate_noise_spectra(0, 1515, 500, form='DelSq', component='real')
nt.assert_equal(Dsq[101102101102].shape, (10, 30))
nt.assert_true(Dsq[101102101102][0, 1] < P_N[101102101102][0, 1])

# test a blpair selection
P_N = uvp.generate_noise_spectra(0, 1515, 500, form='Pk', real=True)
P_N = uvp.generate_noise_spectra(0, 1515, 500, form='Pk', component='real')

def test_average_spectra(self):
uvp = copy.deepcopy(self.uvp)
Expand Down
16 changes: 6 additions & 10 deletions hera_pspec/uvpspec.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1570,7 +1570,7 @@ def units(self):

def generate_noise_spectra(self, spw, polpair, Tsys, blpairs=None,
little_h=True, form='Pk', num_steps=2000,
real=True):
component='real'):
"""
Generate the expected 1-sigma noise power spectrum given a selection of
spectral window, system temp., and polarization. This estimate is
Expand All @@ -1583,7 +1583,8 @@ def generate_noise_spectra(self, spw, polpair, Tsys, blpairs=None,
calculated from pspecbeam with noise_scalar = True, integration_time is
in seconds and comes from self.integration_array and Nincoherent is the
number of incoherent averaging samples and comes from
self.nsample_array.
self.nsample_array. If component is 'real' or 'imag', P_N is divided by
an additional factor of sqrt(2).
If the polarizations specified are pseudo Stokes pol (I, Q, U or V)
then an extra factor of 2 is divided.
Expand Down Expand Up @@ -1626,10 +1627,8 @@ def generate_noise_spectra(self, spw, polpair, Tsys, blpairs=None,
Number of frequency bins to use in integrating power spectrum
scalar in pspecbeam. Default: 2000.
real : bool, optional
If True assumes the real component of complex power spectrum is
used, and will divide P_N by an extra sqrt(2). Otherwise, assume
power spectra are complex and keep P_N as is. Default: True.
component : str, options=['real', 'imag', 'abs']
If component is real or imag, divide by an extra factor of sqrt(2)
Returns
-------
Expand Down Expand Up @@ -1689,15 +1688,12 @@ def generate_noise_spectra(self, spw, polpair, Tsys, blpairs=None,

# Get noise power spectrum
pn = noise.calc_P_N(scalar, Tsys, t_int, k=k,
Nincoherent=n_samp, form=form)
Nincoherent=n_samp, form=form, component=component)

# Put into appropriate form
if form == 'Pk':
pn = np.ones(len(dlys), np.float) * pn

if real:
pn /= np.sqrt(2) # if real divide by sqrt(2)

# append to P_blp
P_blp.append(pn)

Expand Down

0 comments on commit 8c6f8f5

Please sign in to comment.