From a0d7abc69e52b2677b9b09e4a2171da6effe95e5 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Aur=C3=A9lien=20Ponte?= Date: Tue, 25 Jun 2019 14:59:05 +0200 Subject: [PATCH 1/4] init implementation of complex mtspec --- mtspec/multitaper.py | 35 ++++++++++++++++++++++++--------- mtspec/tests/test_multitaper.py | 22 +++++++++++++++++++++ 2 files changed, 48 insertions(+), 9 deletions(-) diff --git a/mtspec/multitaper.py b/mtspec/multitaper.py index 8ebbcc9..741d5ff 100644 --- a/mtspec/multitaper.py +++ b/mtspec/multitaper.py @@ -93,9 +93,13 @@ def mtspec(data, delta, time_bandwidth, nfft=None, number_of_tapers=None, # Depending if nfft is specified or not initialte MtspecTytpe # for mtspec_pad_ or mtspec_d_ + complex = any(np.iscomplex(data)) if nfft is None or nfft == npts: nfft = npts - mt = _MtspecType("float64") # mtspec_d_ + if complex: + mt = _MtspecType("complex64") # mtspec_c_ + else: + mt = _MtspecType("float64") # mtspec_d_ else: mt = _MtspecType("float32") # mtspec_pad_ quadratic = False @@ -103,9 +107,15 @@ def mtspec(data, delta, time_bandwidth, nfft=None, number_of_tapers=None, if number_of_tapers is None: number_of_tapers = int(2 * time_bandwidth) - 1 # Transform the data to work with the library. - data = np.require(data, dtype=mt.float, requirements=[mt.order]) + if complex: + data = np.require(data, dtype=mt.complex, requirements=[mt.order]) + else: + data = np.require(data, dtype=mt.real, requirements=[mt.order]) # Get some information necessary for the call to the Fortran library. - number_of_frequency_bins = int(nfft / 2) + 1 + if complex: + number_of_frequency_bins = nfft + else: + number_of_frequency_bins = int(nfft / 2) + 1 # Create output arrays. spectrum = mt.empty(number_of_frequency_bins) frequency_bins = mt.empty(number_of_frequency_bins) @@ -155,7 +165,7 @@ def mtspec(data, delta, time_bandwidth, nfft=None, number_of_tapers=None, fcrit = None # Call the library. Fortran passes pointers! args = [C.byref(C.c_int(npts)), C.byref(C.c_int(nfft)), - C.byref(mt.c_float(delta)), mt.p(data), + C.byref(mt.c_float(delta)), mt.p(data, complex), C.byref(mt.c_float(time_bandwidth)), C.byref(C.c_int(number_of_tapers)), C.byref(C.c_int(number_of_frequency_bins)), mt.p(frequency_bins), @@ -755,7 +765,8 @@ class _MtspecType(object): the fortran code and provides some helper functions. """ struct = {"float32": (C.c_float, mtspeclib.mtspec_pad_), - "float64": (C.c_double, mtspeclib.mtspec_d_)} + "float64": (C.c_double, mtspeclib.mtspec_d_), + "complex64": (C.c_double, mtspeclib.mtspec_c_)} def __init__(self, dtype): """ @@ -764,10 +775,13 @@ def __init__(self, dtype): :param dtype: 'float32' or 'float64' """ if dtype not in self.struct.keys(): - raise ValueError("dtype must be either 'float32' or 'float64'") + raise ValueError("dtype must be either 'float32' or 'float64'" + + " or 'complex32'") self.float = dtype + self.real = 'float%d' % (float(dtype[-2:])) self.complex = 'complex%d' % (2 * float(dtype[-2:])) - self.c_float = self.struct[dtype][0] + #self.c_float = self.struct[dtype][0] + self.c_float = C.c_double self.pointer = C.POINTER(self.c_float) self.order = "F" self.mtspec = self.struct[dtype][1] @@ -781,9 +795,9 @@ def empty(self, shape, complex=False): """ if complex: return np.empty(shape, dtype=self.complex, order=self.order) - return np.empty(shape, dtype=self.float, order=self.order) + return np.empty(shape, dtype=self.real, order=self.order) - def p(self, ndarray): + def p(self, ndarray, complex=False): """ A wrapper around ctypes.data_as which automatically sets the correct type. Returns none if ndarray is None. @@ -793,4 +807,7 @@ def p(self, ndarray): # short variable name for passing as argument in function calls if ndarray is None: return None + #if complex: return ndarray.ctypes.data_as(self.pointer) + #else: + # return ndarray.ctypes.data_as(self.pointer) diff --git a/mtspec/tests/test_multitaper.py b/mtspec/tests/test_multitaper.py index f396cc9..ed0c70d 100644 --- a/mtspec/tests/test_multitaper.py +++ b/mtspec/tests/test_multitaper.py @@ -71,6 +71,28 @@ def test_multitaper_spectrum(self): np.testing.assert_almost_equal(freq, freq2) np.testing.assert_almost_equal(spec / spec, spec2 / spec, 5) + def test_multitaper_spectrum_complex_input(self): + """ + Test for mtspec. The result is compared to the output of + test_recreatePaperFigures.py in the same directory. This is assumed to + be correct because they are identical to the figures in the paper on + the machine that created these. + """ + data = _load_mtdata('PASC.dat.gz') + # Calculate the spectra. + spec, freq = mtspec(data+1j*data, 1.0, 4.5, number_of_tapers=5) + # No NaNs are supposed to be in the output. + self.assertEqual(np.isnan(spec).any(), False) + self.assertEqual(np.isnan(spec).any(), False) + # Load the good data. + datafile = os.path.join(os.path.dirname(__file__), 'data', + 'multitaper.npz') + spec2 = np.load(datafile)['spec'] + freq2 = np.arange(43201) * 1.15740741e-05 + # Compare, normalize for subdigit comparision + np.testing.assert_almost_equal(freq, freq2) + np.testing.assert_almost_equal(spec / spec, 2*spec2 / spec, 5) + def test_multitaper_spectrum_optional_output(self): """ Test for mtspec. The result is compared to the output of From c78ccd70f8d927b27eb5d51db881c3b14bc8026a Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Aur=C3=A9lien=20Ponte?= Date: Wed, 26 Jun 2019 15:33:58 +0200 Subject: [PATCH 2/4] iteration: mtspec_c in/out variables float32 --- mtspec/multitaper.py | 25 ++++++++++++------------- 1 file changed, 12 insertions(+), 13 deletions(-) diff --git a/mtspec/multitaper.py b/mtspec/multitaper.py index 741d5ff..6aed111 100644 --- a/mtspec/multitaper.py +++ b/mtspec/multitaper.py @@ -96,13 +96,13 @@ def mtspec(data, delta, time_bandwidth, nfft=None, number_of_tapers=None, complex = any(np.iscomplex(data)) if nfft is None or nfft == npts: nfft = npts - if complex: - mt = _MtspecType("complex64") # mtspec_c_ - else: - mt = _MtspecType("float64") # mtspec_d_ + mt = _MtspecType("float64") # mtspec_d_ else: mt = _MtspecType("float32") # mtspec_pad_ quadratic = False + if complex: + mt = _MtspecType("float32") # mtspec_d_ + mt.mtspec = mtspeclib.mtspec_c_ # Use the optimal number of tapers in case no number is specified. if number_of_tapers is None: number_of_tapers = int(2 * time_bandwidth) - 1 @@ -765,8 +765,7 @@ class _MtspecType(object): the fortran code and provides some helper functions. """ struct = {"float32": (C.c_float, mtspeclib.mtspec_pad_), - "float64": (C.c_double, mtspeclib.mtspec_d_), - "complex64": (C.c_double, mtspeclib.mtspec_c_)} + "float64": (C.c_double, mtspeclib.mtspec_d_)} def __init__(self, dtype): """ @@ -779,9 +778,8 @@ def __init__(self, dtype): + " or 'complex32'") self.float = dtype self.real = 'float%d' % (float(dtype[-2:])) - self.complex = 'complex%d' % (2 * float(dtype[-2:])) - #self.c_float = self.struct[dtype][0] - self.c_float = C.c_double + self.complex = 'complex%d' % (2*float(dtype[-2:])) + self.c_float = self.struct[dtype][0] self.pointer = C.POINTER(self.c_float) self.order = "F" self.mtspec = self.struct[dtype][1] @@ -807,7 +805,8 @@ def p(self, ndarray, complex=False): # short variable name for passing as argument in function calls if ndarray is None: return None - #if complex: - return ndarray.ctypes.data_as(self.pointer) - #else: - # return ndarray.ctypes.data_as(self.pointer) + if complex: + pointer_complex = C.POINTER(2*self.c_float) + return ndarray.ctypes.data_as(pointer_complex) + else: + return ndarray.ctypes.data_as(self.pointer) From 4715559e0578a7c9f0cc0cd1ce401fafa06cc656 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Aur=C3=A9lien=20Ponte?= Date: Wed, 26 Jun 2019 15:37:08 +0200 Subject: [PATCH 3/4] minor clean up --- mtspec/multitaper.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/mtspec/multitaper.py b/mtspec/multitaper.py index 6aed111..919af2f 100644 --- a/mtspec/multitaper.py +++ b/mtspec/multitaper.py @@ -777,7 +777,7 @@ def __init__(self, dtype): raise ValueError("dtype must be either 'float32' or 'float64'" + " or 'complex32'") self.float = dtype - self.real = 'float%d' % (float(dtype[-2:])) + #self.real = 'float%d' % (float(dtype[-2:])) self.complex = 'complex%d' % (2*float(dtype[-2:])) self.c_float = self.struct[dtype][0] self.pointer = C.POINTER(self.c_float) @@ -793,7 +793,7 @@ def empty(self, shape, complex=False): """ if complex: return np.empty(shape, dtype=self.complex, order=self.order) - return np.empty(shape, dtype=self.real, order=self.order) + return np.empty(shape, dtype=self.float, order=self.order) def p(self, ndarray, complex=False): """ From 0e64430d6f4405a0f12dc608d7255565637dbd33 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Aur=C3=A9lien=20Ponte?= Date: Mon, 15 Jul 2019 12:10:05 +0200 Subject: [PATCH 4/4] WIP: close but frequencies are not right --- mtspec/multitaper.py | 7 +++++-- 1 file changed, 5 insertions(+), 2 deletions(-) diff --git a/mtspec/multitaper.py b/mtspec/multitaper.py index 919af2f..ff10050 100644 --- a/mtspec/multitaper.py +++ b/mtspec/multitaper.py @@ -101,7 +101,7 @@ def mtspec(data, delta, time_bandwidth, nfft=None, number_of_tapers=None, mt = _MtspecType("float32") # mtspec_pad_ quadratic = False if complex: - mt = _MtspecType("float32") # mtspec_d_ + mt = _MtspecType("float32") # mtspec_c_ mt.mtspec = mtspeclib.mtspec_c_ # Use the optimal number of tapers in case no number is specified. if number_of_tapers is None: @@ -110,7 +110,7 @@ def mtspec(data, delta, time_bandwidth, nfft=None, number_of_tapers=None, if complex: data = np.require(data, dtype=mt.complex, requirements=[mt.order]) else: - data = np.require(data, dtype=mt.real, requirements=[mt.order]) + data = np.require(data, dtype=mt.float, requirements=[mt.order]) # Get some information necessary for the call to the Fortran library. if complex: number_of_frequency_bins = nfft @@ -175,6 +175,7 @@ def mtspec(data, delta, time_bandwidth, nfft=None, number_of_tapers=None, mt.p(eigenspectra), rshape, mt.p(f_statistics), fcrit, None] # diffrent arguments, depending on mtspec_pad_ or mtspec_d_, adapt if npts == nfft: + #print('npts == nfft') args.pop(1) # finally call the shared library function @@ -779,6 +780,8 @@ def __init__(self, dtype): self.float = dtype #self.real = 'float%d' % (float(dtype[-2:])) self.complex = 'complex%d' % (2*float(dtype[-2:])) + #self.complex = 'complex%d' % (float(dtype[-2:])) + # aboves leads to: TypeError: data type "complex32" not understood self.c_float = self.struct[dtype][0] self.pointer = C.POINTER(self.c_float) self.order = "F"