diff --git a/CMakeLists.txt b/CMakeLists.txt index cd3fab4..932b14b 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -21,7 +21,7 @@ elseif(${RPI}) set(DEFAULT_RTAUDIO OFF) set(DEFAULT_PORTAUDIO ON) set(DEFAULT_ULDAQ OFF) - + set(DEFAULT_DOUBLE_PRECISION OFF) else() set(DEFAULT_RTAUDIO OFF) set(DEFAULT_PORTAUDIO ON) diff --git a/cpp_src/dsp/lasp_window.cpp b/cpp_src/dsp/lasp_window.cpp index 9dedb95..fb5d7e4 100644 --- a/cpp_src/dsp/lasp_window.cpp +++ b/cpp_src/dsp/lasp_window.cpp @@ -11,7 +11,7 @@ using std::cerr; using std::endl; // Safe some typing. Linspace form 0 up to (and NOT including N). -#define lin0N arma::linspace(0, N - 1, N) +#define lin0N arma::linspace(0, N - 1, N) vd Window::hann(const us N) { return arma::pow(arma::sin((arma::datum::pi/N) * lin0N), 2); @@ -24,8 +24,8 @@ vd Window::blackman(const us N) { d a0 = 7938. / 18608.; d a1 = 9240. / 18608.; d a2 = 1430. / 18608.; - return a0 - a1 * d_cos((2 * number_pi/N) * lin0N) + - a2 * d_cos((4 * number_pi / N)* lin0N ); + return a0 - a1 * arma::cos((2 * number_pi/N) * lin0N) + + a2 * arma::cos((4 * number_pi / N)* lin0N ); } vd Window::rectangular(const us N) { return arma::ones(N); } diff --git a/cpp_src/lasp_cpp.cpp b/cpp_src/lasp_cpp.cpp index d5c7f76..481d690 100644 --- a/cpp_src/lasp_cpp.cpp +++ b/cpp_src/lasp_cpp.cpp @@ -44,6 +44,12 @@ void init_siggen(py::module &m); PYBIND11_MODULE(lasp_cpp, m) { +#if LASP_DOUBLE_PRECISION == 1 + m.attr("LASP_DOUBLE_PRECISION") = true; +#else + m.attr("LASP_DOUBLE_PRECISION") = false; +#endif + init_dsp(m); init_deviceinfo(m); init_daqconfiguration(m); @@ -51,6 +57,5 @@ PYBIND11_MODULE(lasp_cpp, m) { init_streammgr(m); init_datahandler(m); init_siggen(m); - } -/** @} */ +/** @} */ \ No newline at end of file diff --git a/python_src/lasp/filter/fir_design.py b/python_src/lasp/filter/fir_design.py index 5dd26bc..ad742fb 100644 --- a/python_src/lasp/filter/fir_design.py +++ b/python_src/lasp/filter/fir_design.py @@ -12,6 +12,7 @@ __all__ = ['freqResponse', 'bandpass_fir_design', 'lowpass_fir_design', import numpy as np from scipy.signal import freqz, hann, firwin2 +from ..lasp_config import empty def freqResponse(fs, freq, coefs_b, coefs_a=1.): @@ -44,7 +45,7 @@ def bandpass_fir_design(L, fs, fl, fu, window=hann): Omg2 = 2*np.pi*fu/fs Omg1 = 2*np.pi*fl/fs - fir = np.empty(L, dtype=float) + fir = empty(L, dtype=float) # First Create ideal band-pass filter fir[L//2] = (Omg2-Omg1)/np.pi @@ -64,7 +65,7 @@ def lowpass_fir_design(L, fs, fc, window=hann): " than upper cut-off" Omgc = 2*np.pi*fc/fs - fir = np.empty(L, dtype=float) + fir = empty(L, dtype=float) # First Create ideal band-pass filter fir[L//2] = Omgc/np.pi diff --git a/python_src/lasp/lasp_config.py b/python_src/lasp/lasp_config.py index 46e7af8..24074e4 100644 --- a/python_src/lasp/lasp_config.py +++ b/python_src/lasp/lasp_config.py @@ -6,8 +6,14 @@ Author: J.A. de Jong - ASCEE Description: LASP configuration """ import numpy as np +from .lasp_cpp import LASP_DOUBLE_PRECISION -LASP_NUMPY_FLOAT_TYPE = np.float64 +if LASP_DOUBLE_PRECISION: + LASP_NUMPY_FLOAT_TYPE = np.float64 + LASP_NUMPY_COMPLEX_TYPE = np.float128 +else: + LASP_NUMPY_FLOAT_TYPE = np.float32 + LASP_NUMPY_COMPLEX_TYPE = np.float64 def zeros(shape): diff --git a/python_src/lasp/lasp_measurement.py b/python_src/lasp/lasp_measurement.py index 9fb3614..5259575 100644 --- a/python_src/lasp/lasp_measurement.py +++ b/python_src/lasp/lasp_measurement.py @@ -64,6 +64,7 @@ from .lasp_version import LASP_VERSION_MAJOR, LASP_VERSION_MINOR from .lasp_cpp import Window, DaqChannel, AvPowerSpectra from typing import List from functools import lru_cache +from .lasp_config import ones # Measurement file extension MEXT = 'h5' @@ -223,7 +224,7 @@ class IterData(IterRawData): def __init__(self, fa, channels, sensitivity, **kwargs): super().__init__(fa, channels, **kwargs) - self.sens = np.asarray(sensitivity)[self.channels] + self.sens = np.asarray(sensitivity, dtype=LASP_NUMPY_FLOAT_TYPE)[self.channels] assert self.sens.ndim == 1 def __next__(self): @@ -329,10 +330,10 @@ class Measurement: try: sens = f.attrs["sensitivity"] self._sens = ( - sens * np.ones(self.nchannels) if isinstance(sens, float) else sens + sens * ones(self.nchannels) if isinstance(sens, float) else sens ) except KeyError: - self._sens = np.ones(self.nchannels) + self._sens = ones(self.nchannels) # The time is cached AND ALWAYS ASSUMED TO BE AN IMMUTABLE OBJECT. # It is also cached. Changing the measurement timestamp should not @@ -885,7 +886,7 @@ class Measurement: """ if isinstance(sens, float): # Put all sensitivities equal - sens = sens * np.ones(self.nchannels) + sens = sens * ones(self.nchannels) elif isinstance(sens, list): sens = np.asarray(sens) @@ -1199,7 +1200,7 @@ class Measurement: nchannels = 1 nframes = len(data) data = data[:, np.newaxis] - sensitivity = np.ones(nchannels) + sensitivity = ones(nchannels) with h5.File(newfn, "w") as hf: hf.attrs["samplerate"] = samplerate diff --git a/python_src/lasp/lasp_octavefilter.py b/python_src/lasp/lasp_octavefilter.py index 2a183ac..f8aeb43 100644 --- a/python_src/lasp/lasp_octavefilter.py +++ b/python_src/lasp/lasp_octavefilter.py @@ -77,6 +77,7 @@ class FirFilterBank: self.fs = fs self.xs = list(range(xmin, xmax + 1)) + raise RuntimeError('Not working code anymore') maxdecimation = self.designer.firDecimation(self.xs[0]) self.decimators = [] @@ -245,7 +246,7 @@ class SosFilterBank: for i, x in enumerate(self.xs): channel = self.designer.createSOSFilter(x) if sos is None: - sos = np.empty((channel.size, len(self.xs))) + sos = empty((channel.size, len(self.xs))) sos[:, i] = channel.flatten() self._fb = BiquadBank(sos) diff --git a/python_src/lasp/lasp_reverb.py b/python_src/lasp/lasp_reverb.py index 9ac5831..88ed057 100644 --- a/python_src/lasp/lasp_reverb.py +++ b/python_src/lasp/lasp_reverb.py @@ -6,6 +6,7 @@ Description: Reverberation time estimation tool using least squares """ from .lasp_common import getTime +from .lasp_config import ones import numpy as np @@ -56,7 +57,7 @@ class ReverbTime: x = self._t[istart:istop][:, np.newaxis] # Solve the least-squares problem, by creating a matrix of - A = np.hstack([x, np.ones(x.shape)]) + A = np.hstack([x, ones(x.shape)]) # print(A.shape) # print(points.shape) diff --git a/python_src/lasp/lasp_slm.py b/python_src/lasp/lasp_slm.py index 63f974b..052ea07 100644 --- a/python_src/lasp/lasp_slm.py +++ b/python_src/lasp/lasp_slm.py @@ -5,6 +5,7 @@ Sound level meter implementation @author: J.A. de Jong - ASCEE """ from .lasp_cpp import cppSLM +from .lasp_config import empty import numpy as np from .lasp_common import (TimeWeighting, FreqWeighting, P_REF) from .filter import SPLFilterDesigner @@ -101,7 +102,7 @@ class SLM: assert fbdesigner.fs == fs sos_firstx = fbdesigner.createSOSFilter(self.xs[0]).flatten() self.nom_txt.append(fbdesigner.nominal_txt(self.xs[0])) - sos = np.empty((sos_firstx.size, nfilters), dtype=float, order='C') + sos = empty((sos_firstx.size, nfilters), dtype=float, order='C') sos[:, 0] = sos_firstx for i, x in enumerate(self.xs[1:]): diff --git a/python_src/lasp/lasp_weighcal.py b/python_src/lasp/lasp_weighcal.py index dcff2bc..38e7246 100644 --- a/python_src/lasp/lasp_weighcal.py +++ b/python_src/lasp/lasp_weighcal.py @@ -47,7 +47,7 @@ class WeighCal: P = 2048 # Filter length (number of taps) - self._firs = np.empty((P, self.nchannels)) + self._firs = empty((P, self.nchannels)) self._fbs = [] for chan in range(self.nchannels): fir = arbitrary_fir_design(fs, P, freq_design, diff --git a/python_src/lasp/tools/tools.py b/python_src/lasp/tools/tools.py index 8caee6b..64e4cd9 100644 --- a/python_src/lasp/tools/tools.py +++ b/python_src/lasp/tools/tools.py @@ -21,6 +21,7 @@ import copy import numpy as np from numpy import log2, pi, sin from ..lasp_cpp import freqSmooth +from ..lasp_config import zeros @unique diff --git a/test/test_fft.py b/test/test_fft.py index 09e0b4e..f4f4176 100644 --- a/test/test_fft.py +++ b/test/test_fft.py @@ -31,7 +31,7 @@ def test_backward_fft(): nfft = 2048 freq = getFreq(nfft, nfft) - # Sig = np.zeros(nfft//2+1, dtype=complex) + # Sig = zeros(nfft//2+1, dtype=complex) Sigr = np.random.randn(nfft//2+1) Sigi = np.random.randn(nfft//2+1)