From 310c70985476cca0ca133e6e965e3ecd81506487 Mon Sep 17 00:00:00 2001
From: Kipp Cannon <kcannon@cita.utoronto.ca>
Date: Thu, 31 May 2018 04:43:42 +0900
Subject: [PATCH] PSDFirKernel: clean up
- remove redundant operations from linear_phase_fir_kernel_to_minimum_phase_whitening_fir_kernel(). lots of rolling arrays one way, then back again, multiplying by normalization constants and taking them back out again.
- NOTE: there was probably a bug in this code, roll-left and roll-right operations were not the inverses of each other, leading to a net shift of the frequency series samples. re-testing of the zero-phase whitening filter will be required after this change.
- NOTE: this also fixes the returned phase vector's length
---
gstlal/python/reference_psd.py | 32 ++++++++++++++++++--------------
1 file changed, 18 insertions(+), 14 deletions(-)
diff --git a/gstlal/python/reference_psd.py b/gstlal/python/reference_psd.py
index 6557df44b7..fd36a8a5c1 100644
--- a/gstlal/python/reference_psd.py
+++ b/gstlal/python/reference_psd.py
@@ -29,7 +29,6 @@
import math
import numpy
import os
-import scipy
try:
from pyfftw.interfaces import scipy_fftpack as fftpack
except ImportError:
@@ -529,7 +528,6 @@ class PSDFirKernel(object):
sampleUnits = lal.Unit("strain")
)
-
# FIXME check for change in length
if self.revplan is None:
self.revplan = lal.CreateReverseCOMPLEX16FFTPlan(len(data), 1)
@@ -646,39 +644,35 @@ class PSDFirKernel(object):
)
lal.COMPLEX16TimeFreqFFT(absX, kernel_tseries, self.fwdplan)
- absX.data.data[:] = numpy.roll(abs(absX.data.data), -working_length // 2 + 1) * sample_rate
+ absX.data.data[:] = abs(absX.data.data)
#
# compute the cepstrum of the kernel (i.e., the iFFT of the
# log of the abs of the FFT of the kernel)
#
- logabsX.data.data[:] = numpy.roll(numpy.log(absX.data.data), +working_length // 2 + 1)
+ logabsX.data.data[:] = numpy.log(absX.data.data)
lal.COMPLEX16FreqTimeFFT(cepstrum, logabsX, self.revplan)
- cepstrum.data.data /= sample_rate
#
- # compute sgn
+ # multiply cepstrum by sgn
#
- sgn = scipy.ones(working_length)
- sgn[0] = 0.
- sgn[(len(sgn) + 1) // 2] = 0.
- sgn[(len(sgn) + 1) // 2:] *= -1
- cepstrum.data.data *= sgn
+ cepstrum.data.data[0] = 0.
+ cepstrum.data.data[working_length // 2] = 0.
+ cepstrum.data.data[working_length // 2 + 1:] = -cepstrum.data.data[working_length // 2 + 1:]
#
# compute theta
#
lal.COMPLEX16TimeFreqFFT(theta, cepstrum, self.fwdplan)
- theta.data.data[:] = numpy.roll(-1.j * theta.data.data, -working_length // 2) * sample_rate
#
# compute minimum phase kernel
#
- absX.data.data[:] = numpy.roll(absX.data.data * scipy.exp(1.j * theta.data.data), +working_length // 2 + 1) / sample_rate
+ absX.data.data *= numpy.exp(theta.data.data)
lal.COMPLEX16FreqTimeFFT(min_phase_kernel, absX, self.revplan)
kernel = min_phase_kernel.data.data.real
@@ -690,11 +684,21 @@ class PSDFirKernel(object):
kernel = kernel[-1::-1]
+ #
+ # compute the gain and phase of the zero-phase
+ # approximation relative to the original linear-phase
+ # filter
+ #
+
+ theta = theta.data.data[working_length // 2:]
+ #gain = numpy.exp(theta.real)
+ phase = -theta.imag
+
#
# done
#
- return kernel, -theta.data.data
+ return kernel, phase
def interpolate_psd(psd, deltaF):
--
GitLab