Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Dont repeat dpss window calc #135

Open
wants to merge 16 commits into
base: main
Choose a base branch
from
18 changes: 10 additions & 8 deletions uvtools/dspec.py
Original file line number Diff line number Diff line change
Expand Up @@ -2182,32 +2182,34 @@ def dpss_operator(x, filter_centers, filter_half_widths, cache=None, eigenval_cu
if xc is None:
xc = x[nf//2]
#determine cutoffs
dpss_vectors = []
for fw in filter_half_widths:
dpss_vectors.append(windows.dpss(nf, nf * df * fw, nf))
if nterms is None:
nterms = []
for fn,fw in enumerate(filter_half_widths):
dpss_vectors = windows.dpss(nf, nf * df * fw, nf)
if not eigenval_cutoff is None:
smat = np.sinc(2 * fw * (xg-yg)) * 2 * df * fw
eigvals = np.sum((smat @ dpss_vectors.T) * dpss_vectors.T, axis=0)
eigvals = np.sum((smat @ dpss_vectors[fn].T) * dpss_vectors[fn].T, axis=0)
nterms.append(np.max(np.where(eigvals>=eigenval_cutoff[fn])))
if not edge_suppression is None:
z0=fw * df
edge_tone=np.exp(-2j*np.pi*np.arange(nf)*z0)
fit_components = dpss_vectors * (dpss_vectors @ edge_tone)
fit_components = dpss_vectors[fn] * (dpss_vectors[fn] @ edge_tone)
#this is a vector of RMS residuals of a tone at the edge of the delay window being fitted between 0 to nf DPSS components.
rms_residuals = np.asarray([ np.sqrt(np.mean(np.abs(edge_tone - np.sum(fit_components[:k],axis=0))**2.)) for k in range(nf)])
nterms.append(np.max(np.where(rms_residuals>=edge_suppression[fn])))
if not avg_suppression is None:
sinc_vector=np.sinc(2 * fw * df * (np.arange(nf)-nf/2.))
sinc_vector = sinc_vector / np.sqrt(np.mean(sinc_vector**2.))
fit_components = dpss_vectors * (dpss_vectors @ sinc_vector)
fit_components = dpss_vectors[fn] * (dpss_vectors[fn] @ sinc_vector)
#this is a vector of RMS residuals of vector with equal contributions from all tones within -fw and fw.
rms_residuals = np.asarray([ np.sqrt(np.mean(np.abs(sinc_vector - np.sum(fit_components[:k],axis=0))**2.)) for k in range(nf)])
nterms.append(np.max(np.where(rms_residuals>=avg_suppression[fn])))
#next, construct A matrix.
amat = []
for fc, fw, nt in zip(filter_centers,filter_half_widths, nterms):
amat.append(np.exp(2j * np.pi * (yg[:,:nt]-xc) * fc ) * windows.dpss(nf, nf * df * fw, nt).T )
for fn, (fc, fw, nt) in enumerate(zip(filter_centers,filter_half_widths, nterms)):
amat.append(np.exp(2j * np.pi * (yg[:,:nt]-xc) * fc ) * dpss_vectors[fn][:nt].T )
if len(amat) > 1:
amat = np.hstack(amat)
else:
Expand Down Expand Up @@ -2253,9 +2255,9 @@ def dft_operator(x, filter_centers, filter_half_widths,
if cache is None:
cache = {}
#if no fundamental fourier period is provided, set fundamental period equal to measurement
#bandwidth.
#bandwidth * 2.
if fundamental_period is None:
fundamental_period = np.median(np.diff(x)) * len(x)
fundamental_period = np.median(np.diff(x)) * len(x) * 2
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I don't understand this change well enough to comment on it, so it's a bit concerning that it's not covered by unit tests.

if xc is None:
xc = x[int(np.round(len(x)/2))]
if isinstance(filter_centers, float):
Expand Down
8 changes: 6 additions & 2 deletions uvtools/tests/test_dspec.py
Original file line number Diff line number Diff line change
Expand Up @@ -213,7 +213,7 @@ def test_dft_operator():
freqs = np.arange(-NF/2, NF/2)*DF + 150e6
#test dft_operator by checking whether
#it gives us expected values.
fop = dspec.dft_operator(freqs, 0., 1e-6)
fop = dspec.dft_operator(freqs, 0., 1e-6, fundamental_period=len(freqs) * np.mean(np.diff(freqs)))
fg, dg = np.meshgrid(freqs-150e6, np.arange(-10, 10) * (1./DF/NF) , indexing='ij')
y = np.exp(2j * np.pi * fg * dg )
np.testing.assert_allclose(fop, y)
Expand Down Expand Up @@ -253,7 +253,11 @@ def test_dpss_operator():

dpss_mat = windows.dpss(NF, NF * DF * 100e-9, ncolmax).T
for m in range(ncolmax):
np.testing.assert_allclose(amat4[:,m], dpss_mat[:,m])
# deal with -1 degeneracy which can come up since vectors are obtained from spectral decomposition.
# and there is degeneracy between +/-1 eigenval and eigenvector.
# This simply effects the sign of the vector to be fitted
# and since fitting will cancel by obtaining a -1 on the coefficient, this sign is meaningless
assert np.allclose(amat4[:, m], dpss_mat[:, m]) or np.allclose(amat4[:, m], -dpss_mat[:, m])


def test_fit_solution_matrix():
Expand Down