move to a more standard python packaging structure
Build Debian Packages / build (bookworm, debian12) (push) Successful in 13m3s
Build Debian Packages / build (trixie, debian13) (push) Successful in 13m46s
Build Debian Packages / build (bullseye, debian11) (push) Has been cancelled

This commit is contained in:
2026-04-07 17:09:10 +02:00
parent 1322fd3835
commit 6abb338c4a
43 changed files with 245 additions and 83 deletions
+558
View File
@@ -0,0 +1,558 @@
# -*- coding: iso-8859-1 -*-
from .Resultable import Resultable
from .Drawable import Drawable
from .Signalpath import Signalpath
from .DamarisFFT import DamarisFFT
import threading
import numpy
import sys
import datetime
import tables
#############################################################################
# #
# Name: Class ADC_Result #
# #
# Purpose: Specialised class of Resultable and Drawable #
# Contains recorded ADC Data #
# #
#############################################################################
class ADC_Result(Resultable, Drawable, DamarisFFT, Signalpath):
"""
Represents the result of an ADC, encapsulating data and metadata
for processing, visualization, and export.
This class combines data storage and manipulation functionality with interfaces for
drawing and result processing. It manages time-series data across multiple channels,
supports dynamic resizing of datasets, and provides mechanisms to export data in
various formats. Its key roles include ADC result storage, metadata management, and
integration with external systems.
Attributes:
xlabel: Label for the x-axis, used in visualization.
ylabel: Label for the y-axis, used in visualization.
lock: A threading lock for synchronizing access to the data.
nChannels: Number of data channels in the result set.
sampling_rate: Sampling frequency of the ADC data.
job_id: identifier for the job that generated this result.
"""
def __init__(self, x = None, y:list = None, index = None, sampl_freq = None, desc = None, job_id = None, job_date = None):
Resultable.__init__(self)
Drawable.__init__(self)
# Title of this accumulation: set Values: Job-ID and Description (plotted in GUI -> look Drawable)
# Is set in ResultReader.py (or in copy-construktor)
self.__title_pattern = "ADC-Result: job_id = %s, desc = %s"
# Axis-Labels (inherited from Drawable)
self.xlabel = "Time (s)"
self.ylabel = "Samples [Digits]"
self.lock=threading.RLock()
self.nChannels = 0
# using no argument for initialization
if (x is None) and (y is None) and (index is None) and (sampl_freq is None) and (desc is None) and (job_id is None) and (job_date is None):
self.cont_data = False
self.sampling_rate = 0
self.index = []
self.x = []
self.y = []
# using all arguments for initialization
elif (x is not None) and (y is not None) and (index is not None) and (sampl_freq is not None) and (desc is not None) and (job_id is not None) and (job_date is not None):
# TODO: insure integer calculations for ADC_Result operations.
#for ch in y:
# if not numpy.issubdtype(ch.dtype, numpy.integer):
# raise TypeError("TypeError: ADC_Result y data must be a list with integer type channels")
self.x = x
self.y = y
self.index = index
self.sampling_rate = sampl_freq
self.cont_data = True
self.description = desc
self.job_id = job_id
self.job_date = job_date
title="ADC-Result: job-id=%d"%int(self.job_id)
if len(self.description)>0:
for k,v in self.description.items():
# string keys can be made invisible by adding two underscores in front of them
if not (type(k) in (str,) and k[0] == '_' and k[1] == '_'):
title+=", %s=%s"%(k,v)
self.set_title(title)
else:
raise ValueError("Wrong usage of __init__!")
def create_data_space(self, channels, samples):
"Initialises the internal data-structures"
if self.contains_data():
print("Warning ADC-Result: Tried to run \"create_data_space()\" more than once.")
return
raise ValueError("ValueError: You cant create an ADC-Result with less than 1 sample!")
for i in range(channels):
self.y.append(numpy.zeros((samples,), dtype="int16"))
self.x = numpy.zeros((samples,), dtype="float32")
self.index.append((0, samples-1))
self.cont_data = True
def contains_data(self):
"""Returns true if ADC_Result contains data. (-> create_data_space() was called)"""
return self.cont_data
def add_sample_space(self, samples):
"Adds space for n samples, where n can also be negative (deletes space). New space is filled up with \"0\""
self.lock.acquire()
if not self.cont_data:
print("Warning ADC-Result: Tried to resize empty array!")
return
length = len(self.y[0])
self.x = numpy.resize(self.x, (length+samples))
for i in range(self.get_number_of_channels()):
self.y[i] = numpy.resize(self.y[i], (length+samples))
self.index.append((length, len(self.y[0])-1))
self.lock.release()
def get_result_by_index(self, index):
self.lock.acquire()
try:
start = self.index[index][0]
end = self.index[index][1]
except:
self.lock.release()
raise
tmp_x = self.x[start:end+1].copy()
tmp_y = []
for i in range(self.get_number_of_channels()):
tmp_y.append(self.y[i][start:end+1].copy())
r = ADC_Result(x = tmp_x, y = tmp_y, index = [(0,len(tmp_y[0])-1)], sampl_freq = self.sampling_rate, desc = self.description, job_id = self.job_id, job_date = self.job_date)
self.lock.release()
return r
def get_sampling_rate(self):
"""Returns the samplingfrequency"""
return self.sampling_rate + 0
def set_sampling_rate(self, hz):
"""Sets the samplingfrequency in hz"""
self.sampling_rate = float(hz)
def get_nChannels(self):
"""Gets the number of channels"""
return self.nChannels + 0
def set_nChannels(self, channels):
"""Sets the number of channels"""
self.nChannels = int(channels)
def get_index_bounds(self, index):
"Returns a tuple with (start, end) of the wanted result"
return self.index[index]
def uses_statistics(self):
return False
def write_to_csv(self, destination=sys.stdout, delimiter=" "):
"""
writes the data to a file or to sys.stdout
destination can be a file or a filename
suitable for further processing
"""
# write sorted
the_destination=destination
if type(destination) in (str,):
the_destination=open(destination, "w")
the_destination.write("# adc_result\n")
the_destination.write("# t y0 y1 ...\n")
self.lock.acquire()
try:
xdata=self.get_xdata()
ch_no=self.get_number_of_channels()
ydata=list(map(self.get_ydata, range(ch_no)))
#yerr=map(self.get_yerr, xrange(ch_no))
for i in range(len(xdata)):
the_destination.write("%e"%xdata[i])
for j in range(ch_no):
the_destination.write("%s%e"%(delimiter, ydata[j][i]))
the_destination.write("\n")
the_destination=None
xdata=ydata=None
finally:
self.lock.release()
def write_to_simpson(self, destination=sys.stdout, delimiter=" "):
"""
writes the data to a text file or sys.stdout in Simpson format,
for further processing with the NMRnotebook software;
destination can be a file or a filename
"""
# write sorted
the_destination=destination
if type(destination) in (str,):
the_destination=open(destination, "w")
self.lock.acquire()
try:
xdata=self.get_xdata()
the_destination.write("SIMP\n")
the_destination.write("%s%i%s"%("NP=", len(xdata), "\n"))
the_destination.write("%s%i%s"%("SW=", self.get_sampling_rate(), "\n"))
the_destination.write("TYPE=FID\n")
the_destination.write("DATA\n")
ch_no=self.get_number_of_channels()
ydata=list(map(self.get_ydata, range(ch_no)))
for i in range(len(xdata)):
for j in range(ch_no):
the_destination.write("%g%s"%(ydata[j][i], delimiter))
the_destination.write("\n")
the_destination.write("END\n")
the_destination=None
xdata=ydata=None
finally:
self.lock.release()
def write_to_hdf(self, hdffile, where, name, title, complib=None, complevel=None):
accu_group=hdffile.create_group(where=where,name=name,title=title)
accu_group._v_attrs.damaris_type="ADC_Result"
if self.contains_data():
self.lock.acquire()
try:
# save time stamps
if "job_date" in dir(self) and self.job_date is not None:
accu_group._v_attrs.time="%04d%02d%02d %02d:%02d:%02d.%03d"%(self.job_date.year,
self.job_date.month,
self.job_date.day,
self.job_date.hour,
self.job_date.minute,
self.job_date.second,
self.job_date.microsecond/1000)
if self.description is not None:
for (key,value) in self.description.items():
accu_group._v_attrs.__setattr__("description_"+key,str(value))
accu_group._v_attrs.__setattr__("sampling_rate",self.sampling_rate)
# save interval information
filter=None
if complib is not None:
if complevel is None:
complevel=9
filter=tables.Filters(complevel=complevel,complib=complib,shuffle=1)
index_table=hdffile.create_table(where=accu_group,
name="indices",
description={"start": tables.UInt64Col(),
"length": tables.UInt64Col(),
"start_time": tables.Float32Col(),
"dwelltime": tables.Float32Col()},
title="indices of adc data intervals",
filters=filter,
expectedrows=len(self.index))
index_table.flavor="numpy"
# save channel data
new_row=index_table.row
for i in range(len(self.index)):
new_row["start"]=self.index[i][0]
new_row["dwelltime"]=1.0/self.sampling_rate
new_row["start_time"]=1.0/self.sampling_rate*self.index[i][0]
new_row["length"]=self.index[i][1]-self.index[i][0]+1
new_row.append()
index_table.flush()
new_row=None
index_table=None
# prepare saving data
channel_no=len(self.y)
timedata=numpy.empty((len(self.y[0]),channel_no),
dtype = "int16")
for ch in range(channel_no):
timedata[:,ch]=self.get_ydata(ch)
# save data
time_slice_data=None
if filter is not None:
chunkshape = numpy.shape(timedata)
if len(chunkshape) <= 1:
chunkshape = (min(chunkshape[0],1024*8),)
else:
chunkshape = (min(chunkshape[0],1024*8), chunkshape[1])
if tables.__version__[0]=="1":
time_slice_data=hdffile.create_carray(accu_group,
name="adc_data",
shape=timedata.shape,
atom=tables.Int16Atom(shape=chunkshape,
flavor="numpy"),
filters=filter,
title="adc data")
else:
time_slice_data=hdffile.create_carray(accu_group,
name="adc_data",
shape=timedata.shape,
chunkshape=chunkshape,
atom=tables.Int16Atom(),
filters=filter,
title="adc data")
time_slice_data[:]=timedata
else:
time_slice_data=hdffile.create_array(accu_group,
name="adc_data",
obj=timedata,
title="adc data")
finally:
timedata=None
time_slice_data=None
accu_group=None
self.lock.release()
# Ueberladen von Operatoren und Built-Ins -------------------------------------------------------
def __len__(self):
"Redefining len(ADC_Result obj), returns the number of samples in one channel and 0 without data"
if len(self.y)>0:
return len(self.y[0])
return 0
def __repr__(self):
"""
writes job meta data and data to string returned
"""
tmp_string = "Job ID: " + str(self.job_id) + "\n"
tmp_string += "Job Date: " + str(self.job_date) + "\n"
tmp_string += "Description: " + str(self.description) + "\n"
if len(self.y)>0:
tmp_string += "Indexes: " + str(self.index) + "\n"
tmp_string += "Samples per Channel: " + str(len(self.y[0])) + "\n"
tmp_string += "Samplingfrequency: " + str(self.sampling_rate) + "\n"
tmp_string += "X: " + repr(self.x) + "\n"
for i in range(self.get_number_of_channels()):
tmp_string += ("Y(%d): " % i) + repr(self.y[i]) + "\n"
return tmp_string
def __add__(self, other):
"Redefining self + other (scalar)"
if isinstance(other, int) or isinstance(other, float):
self.lock.acquire()
tmp_y = []
for i in range(self.get_number_of_channels()):
tmp_y.append(numpy.array(self.y[i], dtype="float32") + other)
r = ADC_Result(x = self.x[:], y = tmp_y, index = self.index[:], sampl_freq = self.sampling_rate, desc = self.description, job_id = self.job_id, job_date = self.job_date)
self.lock.release()
return r
else:
raise ValueError(f"ValueError: Cannot add \"{other.__class__}\" to ADC-Result!")
def __radd__(self, other):
"Redefining other (scalar) + self"
return self.__add__(other)
def __sub__(self, other):
"Redefining self - other (scalar)"
if isinstance(other, int) or isinstance(other, float):
self.lock.acquire()
tmp_y = []
for i in range(self.get_number_of_channels()):
tmp_y.append(numpy.array(self.y[i], dtype="float32") - other)
r = ADC_Result(x = self.x[:], y = tmp_y, index = self.index[:], sampl_freq = self.sampling_rate, desc = self.description, job_id = self.job_id, job_date = self.job_date)
self.lock.release()
return r
else:
raise ValueError(f"ValueError: Cannot subtract \"{other.__class__}\" to ADC-Result!")
def __rsub__(self, other):
"Redefining other (scalar) - self"
if isinstance(other, int) or isinstance(other, float):
self.lock.acquire()
tmp_y = []
for i in range(self.get_number_of_channels()):
tmp_y.append(other - numpy.array(self.y[i], dtype="float32"))
r = ADC_Result(x = self.x[:], y = tmp_y, index = self.index[:], sampl_freq = self.sampling_rate, desc = self.description, job_id = self.job_id, job_date = self.job_date)
self.lock.release()
return r
else:
raise ValueError(f"ValueError: Cannot subtract \"{other.__class__}\" to ADC-Result!")
def __mul__(self, other):
"Redefining self * other (scalar)"
if isinstance(other, int) or isinstance(other, float):
self.lock.acquire()
tmp_y = []
for i in range(self.get_number_of_channels()):
tmp_y.append(numpy.array(self.y[i], dtype="float32") * other)
r = ADC_Result(x = self.x[:], y = tmp_y, index = self.index[:], sampl_freq = self.sampling_rate, desc = self.description, job_id = self.job_id, job_date = self.job_date)
self.lock.release()
return r
else:
raise ValueError(f"ValueError: Cannot multiply \"{other.__class__}\" to ADC-Result!")
def __rmul__(self, other):
"Redefining other (scalar) * self"
return self.__mul__(other)
def __pow__(self, other):
"Redefining self ** other (scalar)"
if isinstance(other, int) or isinstance(other, float):
self.lock.acquire()
tmp_y = []
for i in range(self.get_number_of_channels()):
tmp_y.append(numpy.array(self.y[i], dtype="float32") ** other)
r = ADC_Result(x = self.x[:], y = tmp_y, index = self.index[:], sampl_freq = self.sampling_rate, desc = self.description, job_id = self.job_id, job_date = self.job_date)
self.lock.release()
return r
else:
raise ValueError(f"ValueError: Cannot power raise \"{other.__class__}\" to ADC-Result!")
def __truediv__(self, other):
"Redefining other (scalar) / self"
if isinstance(other, int) or isinstance(other, float):
self.lock.acquire()
tmp_y = []
for i in range(self.get_number_of_channels()):
tmp_y.append(other / numpy.array(self.y[i], dtype="float32"))
r = ADC_Result(x = self.x[:], y = tmp_y, index = self.index[:], sampl_freq = self.sampling_rate, desc = self.description, job_id = self.job_id, job_date = self.job_date)
self.lock.release()
return r
else:
raise ValueError(f"ValueError: Cannot divide \"{other.__class__}\" to ADC-Result!")
def __rtruediv__(self, other):
"Redefining other (scalar) / self"
return self.__truediv__(other)
def __floordiv__(self, other):
"Redefining other (scalar) / self"
if isinstance(other, float):
raise ValueError("ValueError: Cannot use floor division (//) on floats! Use \"//\" instead of \"/\"! ")
if isinstance(other, int):
self.lock.acquire()
tmp_y = []
for i in range(self.get_number_of_channels()):
tmp_y.append(other / numpy.array(self.y[i], dtype="float32"))
r = ADC_Result(x = self.x[:], y = tmp_y, index = self.index[:], sampl_freq = self.sampling_rate, desc = self.description, job_id = self.job_id, job_date = self.job_date)
self.lock.release()
return r
else:
raise ValueError(f"ValueError: Cannot divide \"{other.__class__}\" to ADC-Result!")
def __rfloordiv__(self, other):
"Redefining other (scalar) / self"
return self.__floordiv__(other)
def __neg__(self):
"Redefining -self"
self.lock.acquire()
tmp_y = []
for i in range(self.get_number_of_channels()):
tmp_y.append(numpy.array(-self.y[i]))
r = ADC_Result(x = self.x[:], y = tmp_y, index = self.index[:], sampl_freq = self.sampling_rate, desc = self.description, job_id = self.job_id, job_date = self.job_date)
self.lock.release()
return r
def read_from_hdf(hdf_node):
"""
read accumulation data from HDF node and return it.
"""
# formal checks first
if not isinstance(hdf_node, tables.Group):
return None
if hdf_node._v_attrs.damaris_type!="ADC_Result":
return None
if not (hdf_node.__contains__("indices") and hdf_node.__contains__("adc_data")):
return None
# job id and x,y titles are missing
adc=ADC_Result()
# populate description dictionary
adc.description={}
for attrname in hdf_node._v_attrs._v_attrnamesuser:
if attrname.startswith("description_"):
adc.description[attrname[12:]]=hdf_node._v_attrs.__getattr__(attrname)
if "time" in dir(hdf_node._v_attrs):
timestring=hdf_node._v_attrs.__getattr__("time")
adc.job_date=datetime.datetime(int(timestring[:4]), # year
int(timestring[4:6]), # month
int(timestring[6:8]), # day
int(timestring[9:11]), # hour
int(timestring[12:14]), # minute
int(timestring[15:17]), # second
int(timestring[18:21])*1000 # microsecond
)
# start with indices
for r in hdf_node.indices.iterrows():
adc.index.append((r["start"],r["start"]+r["length"]-1))
adc.sampling_rate=1.0/r["dwelltime"]
# now really belief there are no data
if len(adc.index)==0:
adc.cont_data=False
return adc
adc.cont_data=True
# now do the real data
adc_data=hdf_node.adc_data.read()
adc.x=numpy.arange(adc_data.shape[0], dtype="float32")/adc.sampling_rate
for ch in range(adc_data.shape[1]):
adc.y.append(adc_data[:,ch])
return adc
File diff suppressed because it is too large Load Diff
+57
View File
@@ -0,0 +1,57 @@
# -*- coding: iso-8859-1 -*-
from .Resultable import Resultable
#############################################################################
# #
# Name: Class Error_Result #
# #
# Purpose: Specialised class of Resultable #
# Contains occured error-messages from the core #
# #
#############################################################################
class Config_Result(Resultable):
def __init__(self, config = None, desc = None, job_id = None, job_date = None):
Resultable.__init__(self)
if config is None:
self.config = { }
if desc is None:
self.description = { }
self.job_id = job_id
self.job_date = job_date
def get_config_dictionary(self):
return self.config
def set_config_dictionary(self, config):
self.config = config
def get_config(self, key):
if key in self.config:
return self.config[key]
else:
return None
def set_config(self, key, value):
if key in self.config:
print("Warning Config_Result: Key \"%s\" will be overwritten with \"%s\"" % (key, value))
self.config[key] = value
# Ueberladen von Operatoren und Built-Ins -------------------------------------------------------
def __repr__(self):
return str(self.config)
def __str__(self):
return str(self.config)
+215
View File
@@ -0,0 +1,215 @@
import numpy
from . import autophase
class DamarisFFT:
"""
Class for Fourier transforming data.
Provides several helper and apodization functions
"""
def clip( self, start=None, stop=None ):
"""
Method for clipping data, returns only the data between start and stop
start and stop can be either time or frequency.
The unit is automatically determined (Hz or s).
:param float start: beginning of clipping in s
:param float stop: end of clipping in s
"""
# check if start/stop order is properly
if start > stop:
start, stop = stop, start
# do nothing if one uses clip as a "placeholder"
if start is None and stop is None:
return self
if start is None:
start = self.x[ 0 ]
if stop is None:
stop = self.x[ -1 ]
# check if data is fft which changes the start/stop units
if self.xlabel == "Frequency / Hz":
start = self.x.size * (0.5 + start / self.sampling_rate)
stop = self.x.size * (0.5 + stop / self.sampling_rate)
else:
# get the corresponding indices
start *= self.sampling_rate
stop *= self.sampling_rate
# check if boundaries make sense, raise exception otherwise
if numpy.abs( int( start ) - int( stop ) ) <= 0:
raise ValueError( "start stop too close: There are no samples in the given boundaries!" )
# clip the data for each channel
for ch in range( len( self.y ) ):
self.y[ ch ] = self.y[ ch ][ int( start ):int( stop ) ]
self.x = self.x[ int( start ):int( stop ) ]
return self
def baseline( self, last_part=0.1 ):
"""
Correct the baseline of your data by subtracting the mean of the
last_part fraction of your data.
:param float last_part: last section of your timesignal used to calculate baseline
last_part defaults to 0.1, i.e. last 10% of your data
"""
# TODO baseline correction for spectra after:
# Heuer, A; Haeberlen, U.: J. Mag. Res.(1989) 85, Is 1, 79-94
n = int( self.x.size * last_part )
for ch in range( len( self.y ) ):
self.y[ ch ] -= self.y[ ch ][ -n: ].mean( )
return self
def exp_window( self, line_broadening=10 ):
"""
Exponential window function
:param float line_broadening: default 10, line broadening factor in Hz
.. math::
\\exp\\left(-\\pi\\cdot \\textsf{line_broadening} \\cdot t\\right)
"""
apod = numpy.exp( -self.x * numpy.pi * line_broadening )
for i in range( 2 ):
self.y[ i ] = self.y[ i ] * apod
return self
def gauss_window( self, line_broadening=10 ):
"""
Gaussian window function
:param float line_broadening: default 10, line broadening factor in Hz
.. math:: \\exp\\left(- (\\textsf{line_broadening} \\cdot t)^2\\right)
"""
apod = numpy.exp( -(self.x * line_broadening) ** 2 )
for i in range( 2 ):
self.y[ i ] = self.y[ i ] * apod
return self
def dexp_window( self, line_broadening=-10, gaussian_multiplicator=0.3 ):
apod = numpy.exp( -(self.x * line_broadening - gaussian_multiplicator * self.x.max( )) ** 2 )
for i in range( 2 ):
self.y[ i ] = self.y[ i ] * apod
return self
def traf_window( self, line_broadening=10 ):
apod = (numpy.exp( -self.x * line_broadening )) ** 2 / ( (numpy.exp( -self.x * line_broadening )) ** 3
+ (
numpy.exp( -self.x.max( ) * line_broadening )) ** 3 )
for i in range( 2 ):
self.y[ i ] = self.y[ i ] * apod
return self
def hanning_window( self ):
"""
Symmetric centered window (hanning)
"""
apod = numpy.hanning( self.x.size )
for i in range( 2 ):
self.y[ i ] = self.y[ i ] * apod
return self
def hamming_window( self ):
"""
Symmetric centered window (hamming)
"""
apod = numpy.hamming( self.x.size )
for i in range( 2 ):
self.y[ i ] = self.y[ i ] * apod
return self
def blackman_window( self ):
"""
Symmetric centered window (blackmann)
"""
apod = numpy.blackman( self.x.size )
for i in range( 2 ):
self.y[ i ] = self.y[ i ] * apod
return self
def bartlett_window( self ):
"""
Symmetric centered window (bartlett)
"""
apod = numpy.bartlett( self.x.size )
for i in range( 2 ):
self.y[ i ] = self.y[ i ] * apod
return self
def kaiser_window( self, beta=4, use_scipy=None ):
"""
Symmetric centered window (kaiser)
"""
apod = numpy.kaiser( self.x.size, beta )
for i in range( 2 ):
self.y[ i ] = self.y[ i ] * apod
return self
def autophase( self ):
"""
Automatically phases the data to maximize real part.
Works nice with a SNR above 20 dB, i.e.
10 V signal to 0.1 V noise amplitude.
"""
autophase.get_phase( self )
return self
def fft( self, samples=None ):
"""
Calculate the Fourier transform of the data inplace.
For zero filling set **samples** to a value higher than your data length,
smaller values will truncate your data.
:param int samples: default=None, if given, number of samples returned
"""
fft_of_signal = numpy.fft.fft( self.y[ 0 ] + 1j * self.y[ 1 ], n=samples )
fft_of_signal = numpy.fft.fftshift( fft_of_signal )
dwell = 1.0 / self.sampling_rate
n = fft_of_signal.size
fft_frequencies = numpy.fft.fftfreq( n, dwell )
self.x = numpy.fft.fftshift( fft_frequencies )
self.y[ 0 ] = fft_of_signal.real
self.y[ 1 ] = fft_of_signal.imag
self.set_xlabel( "Frequency / Hz" )
return self
def magnitude( self ):
"""
Return absolute signal, i.e.:
.. math::
y[0] &= \\sqrt{y[0]^2 + y[1]^2} \\\\
y[1] &= 0
"""
# this should calculate the absolute value, and set the imag channel to zero
self.y[ 0 ] = numpy.sqrt( self.y[ 0 ] ** 2 + self.y[ 1 ] ** 2 )
self.y[ 1 ] *= 0 # self.y[0].copy()
return self
def ppm(self, f_ref):
"""
Return result scaled to PPM compared to f_ref
:param f_ref: larmor frequency in MHz
:return:
"""
self.x /= f_ref
self.set_xlabel( "PPM" )
return self
+188
View File
@@ -0,0 +1,188 @@
# data pool collects data from data handling script
# provides data to experiment script and display
import sys
import tables
import collections
import threading
import traceback
import io
class DataPool(collections.abc.MutableMapping):
"""
dictionary with sending change events
"""
# supports tranlation from dictionary keys to pytables hdf node names
# taken from: Python Ref Manual Section 2.3: Identifiers and keywords
# things are always prefixed by "dir_" or "dict_"
translation_table=""
for i in range(256):
c=chr(i)
if (c>="a" and c<="z") or \
(c>="A" and c<="Z") or \
(c>="0" and c<="9"):
translation_table+=c
else:
translation_table+="_"
class Event:
access=0
updated_value=1
new_key=2
deleted_key=3
destroy=4
def __init__(self, what, subject="", origin=None):
self.what=what
self.subject=subject
self.origin=origin
def __repr__(self):
return "<DataPool.Event origin=%s what=%d subject='%s'>"%(self.origin, self.what,self.subject)
def copy(self):
return DataPool.Event(self.what+0, self.subject+"", self.origin)
def __init__(self):
self.__mydict={}
self.__dictlock=threading.Lock()
self.__registered_listeners=[]
def __getitem__(self, name):
try:
self.__dictlock.acquire()
return self.__mydict[name]
finally:
self.__dictlock.release()
def __setitem__(self, name, value):
try:
self.__dictlock.acquire()
if name in self.__mydict:
e=DataPool.Event(DataPool.Event.updated_value,name,self)
else:
e=DataPool.Event(DataPool.Event.new_key, name,self)
self.__mydict[name]=value
finally:
self.__dictlock.release()
self.__send_event(e)
def __delitem__(self, name):
try:
self.__dictlock.acquire()
del self.__mydict[name]
finally:
self.__dictlock.release()
self.__send_event(DataPool.Event(DataPool.Event.deleted_key,name,self))
def __iter__(self):
try:
self.__dictlock.acquire()
return iter(self.__mydict)
finally:
self.__dictlock.release()
def __len__(self):
try:
self.__dictlock.acquire()
return len(self.__mydict)
finally:
self.__dictlock.release()
def keys(self):
try:
self.__dictlock.acquire()
return list(self.__mydict.keys())
finally:
self.__dictlock.release()
def __send_event(self, _event):
for listeners in self.__registered_listeners:
listeners(_event.copy())
def __del__(self):
self.__send_event(DataPool.Event(DataPool.Event.destroy))
self.__registered_listeners=None
def write_hdf5(self,hdffile,where="/",name="data_pool", complib=None, complevel=None):
if type(hdffile) is bytes:
dump_file=tables.open_file(hdffile, mode="a")
elif isinstance(hdffile,tables.File):
dump_file=hdffile
else:
raise Exception("expecting hdffile or string")
dump_group=dump_file.create_group(where, name, "DAMARIS data pool")
self.__dictlock.acquire()
dict_keys=list(self.__mydict.keys())
self.__dictlock.release()
try:
for key in dict_keys:
if key.startswith("__"):
continue
dump_dir=dump_group
# walk along the given path and create groups if necessary
namelist = key.split("/")
for part in namelist[:-1]:
dir_part="dir_"+str(part).translate(DataPool.translation_table)
if dir_part not in dump_dir:
dump_dir=dump_file.create_group(dump_dir,name=dir_part,title=part)
else:
if dump_dir._v_children[dir_part]._v_title==part:
dump_dir=dump_dir._v_children[dir_part]
else:
extension_count=0
while dir_part+"_%03d"%extension_count in dump_dir:
extension_count+=1
dump_dir=dump_file.create_group(dump_dir,
name=dir_part+"_%03d"%extension_count,
title=part)
# convert last part of key to a valid name
group_keyname="dict_"+str(namelist[-1]).translate(DataPool.translation_table)
# avoid double names by adding number extension
if group_keyname in dump_dir:
extension_count=0
while group_keyname+"_%03d"%extension_count in dump_dir:
extension_count+=1
group_keyname+="_%03d"%extension_count
self.__dictlock.acquire()
if key not in self.__mydict:
# outdated ...
self.__dictlock.release()
continue
value=self.__mydict[key]
self.__dictlock.release()
# now write data, assuming, the object is constant during write operation
if "write_to_hdf" in dir(value):
try:
value.write_to_hdf(hdffile=dump_file,
where=dump_dir,
name=group_keyname,
title=key,
complib=complib,
complevel=complevel)
except Exception as e:
print("failed to write data_pool[\"%s\"]: %s"%(key,str(e)))
traceback_file=io.StringIO()
traceback.print_tb(sys.exc_info()[2], None, traceback_file)
print("detailed traceback: %s\n"%str(e)+traceback_file.getvalue())
traceback_file=None
else:
print("don't know how to store data_pool[\"%s\"]"%key)
value=None
finally:
dump_group=None
if type(hdffile) is bytes:
dump_file.close()
dump_file=None
def register_listener(self, listening_function):
self.__registered_listeners.append(listening_function)
def unregister_listener(self, listening_function):
if listening_function in self.__registered_listeners:
self.__registered_listeners.remove(listening_function)
+188
View File
@@ -0,0 +1,188 @@
# -*- coding: iso-8859-1 -*-
#############################################################################
# #
# Name: Class Drawable #
# #
# Purpose: Base class of everything plottable #
# #
#############################################################################
class Drawable:
def __init__(self):
# Will be set correctly in one of the subclasses
self.x = []
self.y = []
self.styles = { }
self.xlabel = None
self.ylabel = None
self.title = None
self.legend = { }
self.text = {}
self.xmin = 0
self.xmax = 0
self.ymin = 0
self.ymax = 0
def get_xdata(self):
"Returns a reference to the x-Plotdata (array)"
return self.x
def set_xdata(self, pos, value):
"Sets a point in x"
try:
self.x[pos] = value
except:
raise
def get_ydata(self, channel):
"Returns the y-Plotdata of channel n (array)"
try:
return self.y[channel]
except:
raise
def set_ydata(self, channel, pos, value):
"Sets a point in y"
try:
self.y[channel][pos] = value
except:
raise
def get_number_of_channels(self):
"Returns the number of channels in y"
return len(self.y)
def get_style(self):
"Returns a reference to plot-styles (dictionary)"
return self.styles
def set_style(self, channel, value):
"Sets a channel to a certain plot-style"
if channel in self.styles:
print("Drawable Warning: Style key \"%s\" will be overwritten with \"%s\"" % (str(channel), str(value)))
self.styles[channel] = str(value)
def get_xlabel(self):
"Returns the label for the x-axis"
return self.xlabel
def set_xlabel(self, label):
"Sets the label for the x-axis"
self.xlabel = str(label)
def get_ylabel(self):
"Gets the label for the y-axis"
return self.ylabel
def set_ylabel(self, label):
"Sets the label for the y-axis"
self.ylabel = str(label)
def get_text(self, index):
"Returns labels to be plotted (List)"
if index in self.text:
return self.text[index]
else:
return None
def set_text(self, index, text):
"Sets labels to be plotted "
self.text[index] = str(text)
def get_title(self):
"Returns the title of the plot"
return self.title
def set_title(self, title):
"Sets the title of the plot"
self.title = str(title)
def get_legend(self):
"Returns the legend of the plot (Dictionary)"
return self.legend
def set_legend(self, channel, value):
"Sets the legend of the plot"
if channel in self.legend:
print("Drawable Warning: Legend key \"%s\" will be overwritten with \"%s\"" % (str(channel), str(value)))
self.legend[channel] = str(value)
def get_xmin(self):
"Returns minimun of x"
return self.x.min()
def get_xminpos(self):
"Returns smallest positive value of x"
mask = self.x > 0
return self.x[mask].min()
def set_xmin(self, xmin):
"Sets minimum of x"
self.xmin = xmin
def get_xmax(self):
"Returns maximum of x"
return self.x.max()
def set_xmax(self, xmax):
"Sets maximum of x"
self.xmax = xmax
def get_ymin(self):
"Returns minimum of y"
if isinstance(self.y, list):
return min([yarr.min() for yarr in self.y])
else:
return self.y.min()
def get_yminpos(self):
"Returns smallest positive value of y"
if isinstance(self.y, list):
ymins = []
for ys in self.y:
mask = ys > 0
ymins.append(ys[mask].min())
ymin = min(ymins)
else:
mask = self.y > 0
ymin = self.y[mask].min()
return ymin
def set_ymin(self, ymin):
"Sets minimum of y"
self.ymin = ymin
def get_ymax(self):
"Returns maximimum of y"
if isinstance(self.y, list):
return max([yarr.max() for yarr in self.y])
else:
return self.y.max()
def set_ymax(self, ymax):
"Sets maximum of y"
self.ymax = ymax
+75
View File
@@ -0,0 +1,75 @@
# -*- coding: iso-8859-1 -*-
from .Resultable import Resultable
from .Drawable import Drawable
#############################################################################
# #
# Name: Class Error_Result #
# #
# Purpose: Specialised class of Resultable #
# Contains occured error-messages from the core #
# #
#############################################################################
class Error_Result(Resultable, Drawable):
"""
Specialised class of Resultable
Contains error-messages from the core
"""
def __init__(self, error_msg = None, desc = {}, job_id = None, job_date = None):
Resultable.__init__(self)
Drawable.__init__(self)
if error_msg is not None:
self.error_message = error_msg
self.set_title("Error-Result: %s" % error_msg)
else:
self.error_message = error_msg
self.description = desc
self.job_id = job_id
self.job_date = job_date
def get_error_message(self):
return self.error_message
def set_error_message(self, error_msg):
self.set_title("Error-Result: %s" % error_msg)
self.error_message = error_msg
# No statistics
def uses_statistics(self):
return False
# Nothing to plot
def get_ydata(self):
return [0.0]
# Nothing to plot
def get_xdata(self):
return [0.0]
#overload of operators und built-ins -------------------------------------------------------
def __repr__(self):
tmp_string = "Core error-message: %s" % self.error_message
return tmp_string
def __len__(self):
return len(self.error_message)
def __str__(self):
return self.error_message
# Preventing an error when adding something to an error-result (needed for plotting error-results)
def __add__(self, other):
return self
+78
View File
@@ -0,0 +1,78 @@
# -*- coding: iso-8859-1 -*-
#############################################################################
# #
# Name: Class Errorable #
# #
# Purpose: Base class of everything what could contain a statistic error #
# #
#############################################################################
class Errorable:
"""
Base class for data objects that can have data with errors.
"""
def __init__(self):
# Will be determined in one of the subclasses
self.xerr = []
self.yerr = []
self.error_color = ""
self.bars_above = False
self.n = 0
def get_xerr(self):
"""Returns a reference to x-Error (array)"""
return self.xerr
def set_xerr(self, pos, value):
"""Sets a point in x-Error"""
try:
self.xerr[pos] = value
except:
raise
def get_yerr(self, channel):
"""Returns a list of y-Errors (list of arrays, corresponding channels)"""
try:
return self.yerr[channel]
except:
raise
def set_yerr(self, channel, pos, value):
"""Sets a point in y-Error"""
try:
self.yerr[channel][pos] = value
except:
raise
def get_error_color(self):
"""Returns the error-bar color"""
return self.error_color
def set_error_color(self, color):
"""Sets the error-bar color"""
self.error_color = color
def get_bars_above(self):
"""Gets bars-above property of errorplot"""
return self.bars_above
def set_bars_above(self, bars_above):
"""Sets bars-above property of errorplot"""
self.bars_above = bool(bars_above)
def ready_for_drawing_error(self):
"""Returns true if more than one result have been accumulated"""
return self.n >= 2
+287
View File
@@ -0,0 +1,287 @@
import threading
import math
import types
import sys
import tables
import numpy
import collections
from . import Drawable
## provide gaussian statistics for a series of measured data points
#
# AccumulatedValue provides mean and error of mean after being fed with measured data
# internally it keeps the sum, the sum of squares and the number of data points
class AccumulatedValue:
def __init__(self, mean=None, mean_err=None, n=None):
"""
one value with std. deviation
can be initialized by:
No argument: no entries
one argument: first entry
two arguments: mean and its error, n is set 2
three arguments: already existing statistics defined by mean, mean's error, n
"""
if mean is None:
self.y=0.0
self.y2=0.0
self.n=0
elif mean_err is None and n is None:
self.y=float(mean)
self.y2=self.y**2
self.n=1
elif mean_err is None:
self.n=max(1, int(n))
self.y=float(mean)*self.n
self.y2=(float(mean)**2)*self.n
elif n is None:
self.n=2
self.y=float(mean)*2
self.y2=(float(mean_err)**2+float(mean)**2)*2
else:
self.n=int(n)
self.y=float(mean)*self.n
self.y2=float(mean_err)**2*n*(n-1.0)+float(mean)**2*n
def __add__(self,y):
new_one=AccumulatedValue()
if (type(y) is types.InstanceType and isinstance(y, AccumulatedValue)):
new_one.y=self.y+y.y
new_one.y2=self.y2+y.y2
new_one.n=self.n+y.n
else:
new_one.y=self.y+float(y)
new_one.y2=self.y2+float(y)**2
new_one.n=self.n+1
return new_one
def __iadd__(self,y):
if (type(y) is types.InstanceType and isinstance(y, AccumulatedValue)):
self.y+=y.y
self.y2+=y.y2
self.n+=y.n
else:
self.y+=float(y)
self.y2+=float(y)**2
self.n+=1
return self
def copy(self):
a=AccumulatedValue()
a.y=self.y
a.y2=self.y2
a.n=self.n
return a
def mean(self):
"""
returns the mean of all added/accumulated values
"""
if self.n is None or self.n==0:
return None
else:
return self.y/self.n
def sigma(self):
"""
returns the standard deviation added/accumulated values
"""
if self.n>1:
variance=(self.y2-(self.y**2)/float(self.n))/(self.n-1.0)
if variance<0:
if variance<-1e-20:
print("variance=%g<0! assuming 0"%variance)
return 0.0
return math.sqrt(variance)
elif self.n==1:
return 0.0
else:
return None
def mean_error(self):
"""
returns the mean's error (=std.dev/sqrt(n)) of all added/accumulated values
"""
if self.n>1:
variance=(self.y2-(self.y**2)/float(self.n))/(self.n-1.0)
if variance<0:
if variance<-1e-20:
print("variance=%g<0! assuming 0"%variance)
return 0.0
return math.sqrt(variance/self.n)
elif self.n==1:
return 0.0
else:
return None
def __str__(self):
if self.n==0:
return "no value"
elif self.n==1:
return str(self.y)
else:
return "%g +/- %g (%d accumulations)"%(self.mean(),self.mean_error(),self.n)
def __repr__(self):
return str(self)
class MeasurementResult(Drawable.Drawable, collections.UserDict):
def __init__(self, quantity_name):
"""
convenient accumulation and interface to plot functions
The dictionary must not contain anything but AccumulatedValue instances
"""
Drawable.Drawable.__init__(self)
collections.UserDict.__init__(self)
self.quantity_name=quantity_name
self.lock = threading.RLock()
# get the selected item, if it does not exist, create an empty one
def __getitem__(self, key):
if key not in self:
a=AccumulatedValue()
self.data[float(key)]=a
return a
else:
return self.data[float(key)]
def __setitem__(self,key,value):
if not isinstance(value, AccumulatedValue):
value=AccumulatedValue(float(value))
return collections.UserDict.__setitem__(self,
float(key),
value)
def __add__(self, right_value):
if right_value==0:
return self.copy()
else:
raise Exception("not implemented")
def get_title(self):
return self.quantity_name
def get_xdata(self):
"""
sorted array of all dictionary entries without Accumulated Value objects with n==0
"""
keys=numpy.array([k for k in list(self.data.keys()) if not (isinstance(self.data[k], AccumulatedValue) and self.data[k].n==0)],
dtype="float64")
keys.sort()
return keys
def get_ydata(self):
return self.get_xydata()[1]
def get_xydata(self):
k=self.get_xdata()
v=numpy.array([self.data[key].mean() for key in k], dtype="float64")
return [k,v]
def get_errorplotdata(self):
k=self.get_xdata()
v=numpy.array([self.data[key].mean() for key in k], dtype="float64")
e=numpy.array([self.data[key].mean_error() for key in k], dtype="float64")
return [k,v,e]
def get_lineplotdata(self):
k=self.get_xdata()
v=numpy.array(self.y, dtype="float64")
return [k, v]
def uses_statistics(self):
"""
drawable interface method, returns True
"""
return True
def write_to_csv(self,destination=sys.stdout, delimiter=" "):
"""
writes the data to a file or to sys.stdout
destination can be a file or a filename
suitable for further processing
"""
# write sorted
the_destination=destination
if type(destination) in (str,):
the_destination=open(destination, "w")
the_destination.write("# quantity:"+str(self.quantity_name)+"\n")
the_destination.write("# x y ysigma n\n")
for x in self.get_xdata():
y=self.data[x]
if type(y) in [float, int, int]:
the_destination.write("%e%s%e%s0%s1\n"%(x, delimiter, y, delimiter, delimiter))
else:
the_destination.write("%e%s%e%s%e%s%d\n"%(x,
delimiter,
y.mean(),
delimiter,
y.mean_error(),
delimiter,
y.n))
the_destination=None
def write_to_hdf(self, hdffile, where, name, title, complib=None, complevel=None):
h5_table_format= {
"x" : tables.Float32Col(),
"y" : tables.Float32Col(),
"y_err" : tables.Float32Col(),
"n" : tables.Int64Col()
}
filter=None
if complib is not None:
if complevel is None:
complevel=9
filter=tables.Filters(complevel=complevel,complib=complib,shuffle=1)
mr_table=hdffile.create_table(where=where,name=name,
description=h5_table_format,
title=title,
filters=filter,
expectedrows=len(self))
mr_table.flavor="numpy"
mr_table.attrs.damaris_type="MeasurementResult"
self.lock.acquire()
try:
mr_table.attrs.quantity_name=self.quantity_name
row=mr_table.row
xdata=self.get_xdata()
if xdata.shape[0]!=0:
for x in self.get_xdata():
y=self.data[x]
row["x"]=x
if type(y) in [float, int, int]:
row["y"]=y
row["y_err"]=0.0
row["n"]=1
else:
row["y"]=y.mean()
row["y_err"]=y.mean_error()
row["n"]=y.n
row.append()
finally:
mr_table.flush()
self.lock.release()
def read_from_hdf(hdf_node):
"""
reads a MeasurementResult object from the hdf_node
or None if the node is not suitable
"""
if not isinstance(hdf_node, tables.Table):
return None
if hdf_node._v_attrs.damaris_type!="MeasurementResult":
return None
mr=MeasurementResult(hdf_node._v_attrs.quantity_name)
for r in hdf_node.iterrows():
mr[r["x"]]=AccumulatedValue(r["y"],r["y_err"],r["n"])
return mr
+28
View File
@@ -0,0 +1,28 @@
class Persistance :
def __init__(self, shots):
self.shots = shots
self.accu = 0
self.counter = 0
self.result_list = []
def fade(self, res):
self.counter += 1
if self.accu == 0:
self.accu=res+0
self.result_list.append(res)
if self.counter < 1:
for i,ch in enumerate(self.accu.y):
ch += res.y[i]
elif len(self.result_list) == self.shots:
self.counter = len(self.result_list)
old_result = self.result_list.pop(0)
for i,ch in enumerate(self.accu.y):
ch *= self.shots
ch -= old_result.y[i]
ch += res.y[i]
else:
for i,ch in enumerate(self.accu.y):
ch *= self.counter-1
ch += res.y[i]
self.accu /= self.counter
return self.accu
+65
View File
@@ -0,0 +1,65 @@
# -*- coding: iso-8859-1 -*-
#############################################################################
# #
# Name: Class Resultable #
# #
# Purpose: Base class of everything what could be a core-result #
# #
#############################################################################
class Resultable:
def __init__(self):
self.job_id = None
self.job_date = None
self.description = { }
def get_job_id(self):
"Returns the job-id of this result"
return self.job_id
def set_job_id(self, _id):
"Sets the job-id of this result"
self.job_id = _id
def get_job_date(self):
"Gets the date of this result"
return self.job_date
def set_job_date(self, date):
"Sets the date of this result"
self.job_date = date
def get_description_dictionary(self):
"Returns a reference to the description (Dictionary)"
return self.description
def set_description_dictionary(self, dictionary):
"Sets the entire description"
self.description = dictionary
def get_description(self, key):
"Returns the description value for a given key"
if key in self.description:
return self.description[key]
else:
print("Warning Resultable: No value for key \"%s\". Returned None" % str(key))
return None
def set_description(self, key, value):
"Adds a attribute to the description"
if key in self.description:
print("Warning: Result key \"%s\" will be overwritten with \"%s\"." % (str(key), str(value)))
self.description[key] = value
+37
View File
@@ -0,0 +1,37 @@
import numpy as N
class Signalpath:
def phase(self, degrees):
"""
rotate signal by **degrees** for phase cycling, etc.
:param degrees: rotate signal by this value
:return:
"""
if self.get_number_of_channels() != 2:
raise Exception("rotation defined only for 2 channels")
# simple case 0, 90, 180, 270 degree
reduced_angle = divmod(degrees, 90)
if abs(reduced_angle[1]) < 1e-6:
reduced_angle = reduced_angle[0] % 4
if reduced_angle == 0:
return
elif reduced_angle == 1:
self.y[1] *= -1
self.y = [self.y[1], self.y[0]]
elif reduced_angle == 2:
self.y[0] *= -1
self.y[1] *= -1
elif reduced_angle == 3:
self.y[0] *= -1
self.y = [self.y[1], self.y[0]]
else:
sin_angle = N.sin(degrees / 180.0 * N.pi)
cos_angle = N.cos(degrees / 180.0 * N.pi)
self.y = [cos_angle * self.y[0] - sin_angle * self.y[1],
sin_angle * self.y[0] + cos_angle * self.y[1]]
return self
+31
View File
@@ -0,0 +1,31 @@
# -*- coding: iso-8859-1 -*-
from .Resultable import Resultable
from .Drawable import Drawable
#############################################################################
# #
# Name: Class Temp_Result #
# #
# Purpose: Specialised class of Resultable and Drawable #
# Contains recorded temperature data #
# #
#############################################################################
class TemperatureResult(Resultable, Drawable):
"""
Specialised class of Resultable and Drawable
Contains recorded temperature data
"""
def __init__(self, x = None, y = None, desc = None, job_id = None, job_date = None):
Resultable.__init__(self)
Drawable.__init__(self)
if (x is None) and (y is None) and (desc is None) and (job_id is None) and (job_date is None):
pass
elif (x is not None) and (y is not None) and (desc is not None) and (job_id is not None) and (job_date is not None):
pass
else:
raise ValueError("Wrong usage of __init__!")
+10
View File
@@ -0,0 +1,10 @@
from damaris.data.ADC_Result import ADC_Result
from damaris.data.Accumulation import Accumulation
from damaris.data.MeasurementResult import MeasurementResult, AccumulatedValue
from damaris.data.DataPool import DataPool
from damaris.data.Error_Result import Error_Result
from damaris.data.Config_Result import Config_Result
from damaris.data.Temperature import TemperatureResult
__all__=["ADC_Result", "Accumulation", "MeasurementResult", "AccumulatedValue", "DataPool", "FFT", "Error_Result", "Config_Result", "TemperatureResult" ]
+63
View File
@@ -0,0 +1,63 @@
from scipy.optimize import fmin_powell
import numpy as N
def calculate_entropy(phi, real, imag, gamma, dwell):
"""
Calculates the entropy of the spectrum (real part).
p = phase
gamma should be adjusted such that the penalty and entropy are in the same magnitude
"""
# This is first order phasecorrection
# corr_phase = phi[0]+phi[1]*arange(0,len(signal),1.0)/len(signal) # For 0th and 1st correction
# Zero order phase correction
real_part = real*N.cos(phi)-imag*N.sin(phi)
# Either this for calculating derivatives:
# Zwei-Punkt-Formel
# real_diff = (Re[1:]-Re[:-1])/dwell
# Better this:
# Drei-Punkte-Mittelpunkt-Formel (Ränder werden nicht beachtet)
# real_diff = abs((Re[2:]-Re[:-2])/(dwell*2))
# Even better:
# Fünf-Punkte-Mittelpunkt-Formel (ohne Ränder)
real_diff = N.abs((real_part[:-4]-8*real_part[1:-3]
+8*real_part[3:-1]-2*real_part[4:])/(12*dwell))
# TODO Ränder, sind wahrscheinlich nicht kritisch
# Calculate the entropy
h = real_diff/real_diff.sum()
# Set all h with 0 to 1 (log would complain)
h[h==0]=1
entropy = N.sum(-h*N.log(h))
# My version, according the paper
#penalty = gamma*sum([val**2 for val in Re if val < 0])
# calculate penalty value: a real spectrum should have positive values
if real_part.sum() < 0:
tmp = real_part[real_part<0]
penalty = N.dot(tmp,tmp)
if gamma == 0:
gamma = entropy/penalty
penalty = N.dot(tmp,tmp)*gamma
else:
penalty = 0
#print "Entropy:",entrop,"Penalty:",penalty # Debugging
shannon = entropy+penalty
return shannon
def get_phase(result_object):
global gamma
gamma=0
real = result_object.y[0].copy()
imag = result_object.y[1].copy()
dwell = 1.0/result_object.sampling_rate
# fmin also possible
xopt = fmin_powell( func=calculate_entropy,
x0=N.array([0.0]),
args=(real, imag, gamma, dwell),
disp=0)
result_object.y[0] = real*N.cos(xopt) - imag*N.sin(xopt)
result_object.y[1] = real*N.sin(xopt) + imag*N.cos(xopt)
return result_object