Heroic squashing of spurious furious tab bugs
This commit is contained in:
Stefan Reutter
parent
0a393b0748
commit
bb95347567
@ -6,206 +6,206 @@ import numpy as N
|
||||
import numpy.fft as F
|
||||
|
||||
class FFT:
|
||||
def __init__(self, one_result):
|
||||
# create copy of one_result and work only on the copy
|
||||
# also extract some informations
|
||||
self.the_result = one_result + 0
|
||||
self.timepoints = N.array(one_result.x)
|
||||
self.sampling_rate = one_result.get_sampling_rate()
|
||||
self.data_points = one_result.get_ydata(0).size
|
||||
self.aquisition_time = self.data_points / float(self.sampling_rate)
|
||||
self.the_result.set_xlabel('Frequency [Hz]')
|
||||
|
||||
def write_n(self, afile):
|
||||
filename = open(afile,'w')
|
||||
filename = open(afile,'a')
|
||||
#print self.the_result.get_description_dictionary()
|
||||
#filename.write('%s'%self.get_description_dictionary())
|
||||
for i in range(self.data_points):
|
||||
filename.write('%e\t%e\t%e\n'%(self.the_result.x[i], self.the_result.y[0][i], self.the_result.y[1][i]))
|
||||
filename.close()
|
||||
return self
|
||||
def __init__(self, one_result):
|
||||
# create copy of one_result and work only on the copy
|
||||
# also extract some informations
|
||||
self.the_result = one_result + 0
|
||||
self.timepoints = N.array(one_result.x)
|
||||
self.sampling_rate = one_result.get_sampling_rate()
|
||||
self.data_points = one_result.get_ydata(0).size
|
||||
self.aquisition_time = self.data_points / float(self.sampling_rate)
|
||||
self.the_result.set_xlabel('Frequency [Hz]')
|
||||
|
||||
def write_n(self, afile):
|
||||
filename = open(afile,'w')
|
||||
filename = open(afile,'a')
|
||||
#print self.the_result.get_description_dictionary()
|
||||
#filename.write('%s'%self.get_description_dictionary())
|
||||
for i in range(self.data_points):
|
||||
filename.write('%e\t%e\t%e\n'%(self.the_result.x[i], self.the_result.y[0][i], self.the_result.y[1][i]))
|
||||
filename.close()
|
||||
return self
|
||||
|
||||
def base_corr(self, cutoff=0.3, show=0):
|
||||
"""
|
||||
Subtracts the mean of the last cutoff % of the timsignal
|
||||
to get rid of the DC part in the FFT and returns the
|
||||
new data.
|
||||
If cutoff is not given, the mean of the last 30% will be
|
||||
subtracted.
|
||||
If show=1 the result is return and not the instance. This allows to plot the baseline corrected signal
|
||||
Example:
|
||||
base_corr(cutoff=0.2, show=1)
|
||||
"""
|
||||
last_points = int(cutoff*self.data_points)
|
||||
for i in range(2):
|
||||
self.the_result.y[i] = self.the_result.y[i] - self.the_result.y[i][-last_points:].mean()
|
||||
if show == 1 :
|
||||
return self.the_result
|
||||
return self
|
||||
def base_corr(self, cutoff=0.3, show=0):
|
||||
"""
|
||||
Subtracts the mean of the last cutoff % of the timsignal
|
||||
to get rid of the DC part in the FFT and returns the
|
||||
new data.
|
||||
If cutoff is not given, the mean of the last 30% will be
|
||||
subtracted.
|
||||
If show=1 the result is return and not the instance. This allows to plot the baseline corrected signal
|
||||
Example:
|
||||
base_corr(cutoff=0.2, show=1)
|
||||
"""
|
||||
last_points = int(cutoff*self.data_points)
|
||||
for i in range(2):
|
||||
self.the_result.y[i] = self.the_result.y[i] - self.the_result.y[i][-last_points:].mean()
|
||||
if show == 1 :
|
||||
return self.the_result
|
||||
return self
|
||||
|
||||
def abs_fft(self, points=None, zoom=None,write = 'off'):
|
||||
"""
|
||||
Fourier transforms the timesignal;
|
||||
points is the number of points to transform, if more points given than data points
|
||||
the rest is zero padded
|
||||
|
||||
absfft(points=4096)
|
||||
"""
|
||||
realdata = N.array(self.the_result.y[0])
|
||||
imdata = N.array(self.the_result.y[1])
|
||||
data = realdata + 1j*imdata
|
||||
fftdata = F.fftshift(F.fft(data, points))
|
||||
absfft = N.sqrt(fftdata.real**2 + fftdata.imag**2)
|
||||
# create our x axis
|
||||
n = fftdata.size
|
||||
self.the_result.x = F.fftshift(F.fftfreq(n, 1.0/self.sampling_rate))
|
||||
self.the_result.y[0] = absfft
|
||||
self.the_result.y[1] = N.zeros(n)
|
||||
if write == 'on':
|
||||
return self
|
||||
else:
|
||||
if zoom is None:
|
||||
return self.the_result
|
||||
else:
|
||||
center, width = zoom
|
||||
return self.zoom(self.the_result, center, width)
|
||||
def abs_fft(self, points=None, zoom=None,write = 'off'):
|
||||
"""
|
||||
Fourier transforms the timesignal;
|
||||
points is the number of points to transform, if more points given than data points
|
||||
the rest is zero padded
|
||||
|
||||
absfft(points=4096)
|
||||
"""
|
||||
realdata = N.array(self.the_result.y[0])
|
||||
imdata = N.array(self.the_result.y[1])
|
||||
data = realdata + 1j*imdata
|
||||
fftdata = F.fftshift(F.fft(data, points))
|
||||
absfft = N.sqrt(fftdata.real**2 + fftdata.imag**2)
|
||||
# create our x axis
|
||||
n = fftdata.size
|
||||
self.the_result.x = F.fftshift(F.fftfreq(n, 1.0/self.sampling_rate))
|
||||
self.the_result.y[0] = absfft
|
||||
self.the_result.y[1] = N.zeros(n)
|
||||
if write == 'on':
|
||||
return self
|
||||
else:
|
||||
if zoom is None:
|
||||
return self.the_result
|
||||
else:
|
||||
center, width = zoom
|
||||
return self.zoom(self.the_result, center, width)
|
||||
|
||||
|
||||
def fft(self, points=None, zoom=None, write='off'):
|
||||
realdata = N.array(self.the_result.y[0])
|
||||
imdata = N.array(self.the_result.y[1])
|
||||
data = realdata + 1j*imdata
|
||||
fftdata = F.fftshift(F.fft(data, points))
|
||||
# create our x axis
|
||||
n = fftdata.size
|
||||
self.the_result.x = F.fftshift(F.fftfreq(n, 1.0/self.sampling_rate))
|
||||
self.the_result.y[0] = fftdata.real
|
||||
self.the_result.y[1] = fftdata.imag
|
||||
if write == 'on':
|
||||
return self
|
||||
else:
|
||||
if zoom is None:
|
||||
return self.the_result
|
||||
else:
|
||||
center, width = zoom
|
||||
return self.zoom(self.the_result, center, width)
|
||||
def fft(self, points=None, zoom=None, write='off'):
|
||||
realdata = N.array(self.the_result.y[0])
|
||||
imdata = N.array(self.the_result.y[1])
|
||||
data = realdata + 1j*imdata
|
||||
fftdata = F.fftshift(F.fft(data, points))
|
||||
# create our x axis
|
||||
n = fftdata.size
|
||||
self.the_result.x = F.fftshift(F.fftfreq(n, 1.0/self.sampling_rate))
|
||||
self.the_result.y[0] = fftdata.real
|
||||
self.the_result.y[1] = fftdata.imag
|
||||
if write == 'on':
|
||||
return self
|
||||
else:
|
||||
if zoom is None:
|
||||
return self.the_result
|
||||
else:
|
||||
center, width = zoom
|
||||
return self.zoom(self.the_result, center, width)
|
||||
|
||||
def zoom(self,some_result, center="auto", width=1000):
|
||||
if center == "auto":
|
||||
i_center = int(self.the_result.y[0].argmax())
|
||||
maximum = self.the_result.y[0][i_center]
|
||||
print "Maximum at Frequency:", self.the_result.x[i_center]
|
||||
else:
|
||||
i_center = int(self.data_points/2.0+self.data_points*center/self.sampling_rate)
|
||||
#print "TODO: set width automagically"
|
||||
#if width == "auto":
|
||||
# i_width = int(self.data_points*width)
|
||||
i_width = int(self.data_points*width/self.sampling_rate)
|
||||
some_result.x=some_result.x[i_center-i_width/2:i_center+i_width/2]
|
||||
some_result.y[0]=some_result.y[0][i_center-i_width/2:i_center+i_width/2]
|
||||
some_result.y[1]=some_result.y[1][i_center-i_width/2:i_center+i_width/2]
|
||||
return some_result
|
||||
def zoom(self,some_result, center="auto", width=1000):
|
||||
if center == "auto":
|
||||
i_center = int(self.the_result.y[0].argmax())
|
||||
maximum = self.the_result.y[0][i_center]
|
||||
print "Maximum at Frequency:", self.the_result.x[i_center]
|
||||
else:
|
||||
i_center = int(self.data_points/2.0+self.data_points*center/self.sampling_rate)
|
||||
#print "TODO: set width automagically"
|
||||
#if width == "auto":
|
||||
# i_width = int(self.data_points*width)
|
||||
i_width = int(self.data_points*width/self.sampling_rate)
|
||||
some_result.x=some_result.x[i_center-i_width/2:i_center+i_width/2]
|
||||
some_result.y[0]=some_result.y[0][i_center-i_width/2:i_center+i_width/2]
|
||||
some_result.y[1]=some_result.y[1][i_center-i_width/2:i_center+i_width/2]
|
||||
return some_result
|
||||
|
||||
"""
|
||||
Apodization functions:
|
||||
* exp_window and gauss_window are S/N enhancing,
|
||||
* dexp_window and traf_window are resolution enhancing
|
||||
* standard windows [hamming, hanning, bartlett, blackman, kaiser-bessel] are also available
|
||||
self.timepoints = time points
|
||||
self.aquisition_time = aquisition time (no. samples / sampling_rate)
|
||||
line_broadening = line broadening factor (standard = 10 Hz)
|
||||
gaussian_multiplicator = Gaussian Multiplication Factor for
|
||||
the double exponential apodization
|
||||
function (standard = 0.3)
|
||||
"""
|
||||
def exp_window(self, line_broadening=10, show=0):
|
||||
apod = N.exp(-self.timepoints*line_broadening)
|
||||
for i in range(2):
|
||||
self.the_result.y[i] = self.the_result.y[i]*apod
|
||||
if show == 1 :
|
||||
return self.the_result
|
||||
return self
|
||||
|
||||
def gauss_window(self, line_broadening=10, show=0):
|
||||
apod = N.exp(-(self.timepoints*line_broadening)**2)
|
||||
for i in range(2):
|
||||
self.the_result.y[i] = self.the_result.y[i]*apod
|
||||
if show == 1 :
|
||||
return self.the_result
|
||||
return self
|
||||
|
||||
def dexp_window(self, line_broadening=10, gaussian_multiplicator=0.3, show=0):
|
||||
apod = N.exp(-(self.timepoints*line_broadening - gaussian_multiplicator*self.aquisition_time)**2)
|
||||
for i in range(2):
|
||||
self.the_result.y[i] = self.the_result.y[i]*apod
|
||||
if show == 1:
|
||||
return self.the_result
|
||||
return self
|
||||
|
||||
def traf_window(self, line_broadening=10, show=0):
|
||||
apod = (N.exp(-self.timepoints*line_broadening))**2 / ( (N.exp(-self.timepoints*line_broadening))**3
|
||||
+ (N.exp(-self.aquisition_time*line_broadening))**3 )
|
||||
for i in range(2):
|
||||
self.the_result.y[i] = self.the_result.y[i]*apod
|
||||
if show == 1:
|
||||
return self.the_result
|
||||
return self
|
||||
|
||||
def hanning_window(self, show=0):
|
||||
apod = N.hanning(self.data_points)
|
||||
for i in range(2):
|
||||
self.the_result.y[i] = self.the_result.y[i]*apod
|
||||
if show == 1:
|
||||
return self.the_result
|
||||
return self
|
||||
|
||||
def hamming_window(self, show=0):
|
||||
apod = N.hamming(self.data_points)
|
||||
for i in range(2):
|
||||
self.the_result.y[i] = self.the_result.y[i]*apod
|
||||
if show == 1:
|
||||
return self.the_result
|
||||
return self
|
||||
|
||||
def blackman_window(self, show=0):
|
||||
apod = N.blackman(self.data_points)
|
||||
for i in range(2):
|
||||
self.the_result.y[i] = self.the_result.y[i]*apod
|
||||
if show == 1:
|
||||
return self.the_result
|
||||
return self
|
||||
|
||||
def bartlett_window(self, show=0):
|
||||
apod = N.bartlett(self.data_points)
|
||||
for i in range(2):
|
||||
self.the_result.y[i] = self.the_result.y[i]*apod
|
||||
if show == 1:
|
||||
return self.the_result
|
||||
return self
|
||||
|
||||
def kaiser_window(self, beta=4, show=0, use_scipy=None):
|
||||
if use_scipy == None:
|
||||
# modified Bessel function of zero kind order from somewhere
|
||||
def I_0(x):
|
||||
i0=0
|
||||
fac = lambda n:reduce(lambda a,b:a*(b+1),range(n),1)
|
||||
for n in range(20):
|
||||
i0 += ((x/2.0)**n/(fac(n)))**2
|
||||
return i0
|
||||
|
||||
t = N.arange(self.data_points, type=N.Float) - self.data_points/2.0
|
||||
T = self.data_points
|
||||
# this is the window function array
|
||||
apod = I_0(beta*N.sqrt(1-(2*t/T)**2))/I_0(beta)
|
||||
else:
|
||||
# alternative method using scipy
|
||||
import scipy
|
||||
apod=scipy.kaiser(self.data_points, beta)
|
||||
|
||||
for i in range(2):
|
||||
self.the_result.y[i] = self.the_result.y[i]*apod
|
||||
if show == 1:
|
||||
return self.the_result
|
||||
return self
|
||||
|
||||
"""
|
||||
Apodization functions:
|
||||
* exp_window and gauss_window are S/N enhancing,
|
||||
* dexp_window and traf_window are resolution enhancing
|
||||
* standard windows [hamming, hanning, bartlett, blackman, kaiser-bessel] are also available
|
||||
self.timepoints = time points
|
||||
self.aquisition_time = aquisition time (no. samples / sampling_rate)
|
||||
line_broadening = line broadening factor (standard = 10 Hz)
|
||||
gaussian_multiplicator = Gaussian Multiplication Factor for
|
||||
the double exponential apodization
|
||||
function (standard = 0.3)
|
||||
"""
|
||||
def exp_window(self, line_broadening=10, show=0):
|
||||
apod = N.exp(-self.timepoints*line_broadening)
|
||||
for i in range(2):
|
||||
self.the_result.y[i] = self.the_result.y[i]*apod
|
||||
if show == 1 :
|
||||
return self.the_result
|
||||
return self
|
||||
|
||||
def gauss_window(self, line_broadening=10, show=0):
|
||||
apod = N.exp(-(self.timepoints*line_broadening)**2)
|
||||
for i in range(2):
|
||||
self.the_result.y[i] = self.the_result.y[i]*apod
|
||||
if show == 1 :
|
||||
return self.the_result
|
||||
return self
|
||||
|
||||
def dexp_window(self, line_broadening=10, gaussian_multiplicator=0.3, show=0):
|
||||
apod = N.exp(-(self.timepoints*line_broadening - gaussian_multiplicator*self.aquisition_time)**2)
|
||||
for i in range(2):
|
||||
self.the_result.y[i] = self.the_result.y[i]*apod
|
||||
if show == 1:
|
||||
return self.the_result
|
||||
return self
|
||||
|
||||
def traf_window(self, line_broadening=10, show=0):
|
||||
apod = (N.exp(-self.timepoints*line_broadening))**2 / ( (N.exp(-self.timepoints*line_broadening))**3
|
||||
+ (N.exp(-self.aquisition_time*line_broadening))**3 )
|
||||
for i in range(2):
|
||||
self.the_result.y[i] = self.the_result.y[i]*apod
|
||||
if show == 1:
|
||||
return self.the_result
|
||||
return self
|
||||
|
||||
def hanning_window(self, show=0):
|
||||
apod = N.hanning(self.data_points)
|
||||
for i in range(2):
|
||||
self.the_result.y[i] = self.the_result.y[i]*apod
|
||||
if show == 1:
|
||||
return self.the_result
|
||||
return self
|
||||
|
||||
def hamming_window(self, show=0):
|
||||
apod = N.hamming(self.data_points)
|
||||
for i in range(2):
|
||||
self.the_result.y[i] = self.the_result.y[i]*apod
|
||||
if show == 1:
|
||||
return self.the_result
|
||||
return self
|
||||
|
||||
def blackman_window(self, show=0):
|
||||
apod = N.blackman(self.data_points)
|
||||
for i in range(2):
|
||||
self.the_result.y[i] = self.the_result.y[i]*apod
|
||||
if show == 1:
|
||||
return self.the_result
|
||||
return self
|
||||
|
||||
def bartlett_window(self, show=0):
|
||||
apod = N.bartlett(self.data_points)
|
||||
for i in range(2):
|
||||
self.the_result.y[i] = self.the_result.y[i]*apod
|
||||
if show == 1:
|
||||
return self.the_result
|
||||
return self
|
||||
|
||||
def kaiser_window(self, beta=4, show=0, use_scipy=None):
|
||||
if use_scipy == None:
|
||||
# modified Bessel function of zero kind order from somewhere
|
||||
def I_0(x):
|
||||
i0=0
|
||||
fac = lambda n:reduce(lambda a,b:a*(b+1),range(n),1)
|
||||
for n in range(20):
|
||||
i0 += ((x/2.0)**n/(fac(n)))**2
|
||||
return i0
|
||||
|
||||
t = N.arange(self.data_points, type=N.Float) - self.data_points/2.0
|
||||
T = self.data_points
|
||||
# this is the window function array
|
||||
apod = I_0(beta*N.sqrt(1-(2*t/T)**2))/I_0(beta)
|
||||
else:
|
||||
# alternative method using scipy
|
||||
import scipy
|
||||
apod=scipy.kaiser(self.data_points, beta)
|
||||
|
||||
for i in range(2):
|
||||
self.the_result.y[i] = self.the_result.y[i]*apod
|
||||
if show == 1:
|
||||
return self.the_result
|
||||
return self
|
||||
|
||||
|
||||
@ -4,154 +4,154 @@ import autophase
|
||||
|
||||
class DamarisFFT:
|
||||
def clip(self, start=None, stop=None):
|
||||
"""
|
||||
Method for clipping data, only the timesignal between start and stop
|
||||
is returned.
|
||||
start and stop can be either time or frequency. The unit is automatically determined
|
||||
"""
|
||||
# check if start/stop order is properly
|
||||
if start > stop:
|
||||
# I could swap start/stop actually
|
||||
# TODO swap values?
|
||||
raise
|
||||
# if one uses clip as a "placeholder"
|
||||
if start==None and stop==None:
|
||||
return self
|
||||
"""
|
||||
Method for clipping data, only the timesignal between start and stop
|
||||
is returned.
|
||||
start and stop can be either time or frequency. The unit is automatically determined
|
||||
"""
|
||||
# check if start/stop order is properly
|
||||
if start > stop:
|
||||
# I could swap start/stop actually
|
||||
# TODO swap values?
|
||||
raise
|
||||
# if one uses clip as a "placeholder"
|
||||
if start==None and stop==None:
|
||||
return self
|
||||
|
||||
if start==None:
|
||||
start = 0
|
||||
if stop==None:
|
||||
stop = -1
|
||||
# check if data is fft which changes the start/stop units
|
||||
# TODO should get nicer(failsafe), i.e. flags in the object?
|
||||
if self.xlabel == "Frequency / Hz":
|
||||
isfft = True
|
||||
start = self.x.size*(0.5 + start/self.sampling_rate)
|
||||
stop = self.x.size*(0.5 + stop/self.sampling_rate)
|
||||
else:
|
||||
isfft = False
|
||||
# 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 values in the given boundaries!")
|
||||
for ch in xrange(len(self.y)):
|
||||
# clip the data for each channel
|
||||
# TODO multi records
|
||||
self.y[ch] = self.y[ch][int(start):int(stop)]
|
||||
# TODO what to do with x? Should it start from 0 or from start?
|
||||
# self.x = self.x[:int(stop)-int(start)]
|
||||
self.x = self.x[int(start):int(stop)]
|
||||
return self
|
||||
if start==None:
|
||||
start = 0
|
||||
if stop==None:
|
||||
stop = -1
|
||||
# check if data is fft which changes the start/stop units
|
||||
# TODO should get nicer(failsafe), i.e. flags in the object?
|
||||
if self.xlabel == "Frequency / Hz":
|
||||
isfft = True
|
||||
start = self.x.size*(0.5 + start/self.sampling_rate)
|
||||
stop = self.x.size*(0.5 + stop/self.sampling_rate)
|
||||
else:
|
||||
isfft = False
|
||||
# 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 values in the given boundaries!")
|
||||
for ch in xrange(len(self.y)):
|
||||
# clip the data for each channel
|
||||
# TODO multi records
|
||||
self.y[ch] = self.y[ch][int(start):int(stop)]
|
||||
# TODO what to do with x? Should it start from 0 or from start?
|
||||
# self.x = self.x[:int(stop)-int(start)]
|
||||
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.
|
||||
"""
|
||||
Correct the baseline of your data by subtracting the mean of the
|
||||
last_part fraction of your data.
|
||||
|
||||
last_part defaults to 0.1, i.e. last 10% of your data
|
||||
"""
|
||||
# TODO baselinecorrection for spectra after:
|
||||
# Heuer, A; Haeberlen, U.: J. Mag. Res.(1989) 85, Is 1, 79-94
|
||||
# Should I create an empty object?
|
||||
# I deided to do NOT a copy, but
|
||||
# rather modify the object
|
||||
n = int(self.x.size*last_part)
|
||||
for ch in xrange(len(self.y)):
|
||||
self.y[ch] -= self.y[ch][-n:].mean()
|
||||
# Skip the following due to design reasons
|
||||
# new_object.was_copied = True
|
||||
return self
|
||||
last_part defaults to 0.1, i.e. last 10% of your data
|
||||
"""
|
||||
# TODO baselinecorrection for spectra after:
|
||||
# Heuer, A; Haeberlen, U.: J. Mag. Res.(1989) 85, Is 1, 79-94
|
||||
# Should I create an empty object?
|
||||
# I deided to do NOT a copy, but
|
||||
# rather modify the object
|
||||
n = int(self.x.size*last_part)
|
||||
for ch in xrange(len(self.y)):
|
||||
self.y[ch] -= self.y[ch][-n:].mean()
|
||||
# Skip the following due to design reasons
|
||||
# new_object.was_copied = True
|
||||
return self
|
||||
|
||||
|
||||
"""
|
||||
Apodization functions:
|
||||
* exp_window and gauss_window are S/N enhancing,
|
||||
* dexp_window and traf_window are resolution enhancing
|
||||
* standard windows [hamming, hanning, bartlett, blackman, kaiser-bessel]
|
||||
are also available
|
||||
self.x = time points
|
||||
elf.aquisition_time = aquisition time (no. samples / sampling_rate)
|
||||
line_broadening = line broadening factor (standard = 10 Hz)
|
||||
gaussian_multiplicator = Gaussian Multiplication Factor for
|
||||
the double exponential apodization
|
||||
function (standard = 0.3)
|
||||
"""
|
||||
"""
|
||||
Apodization functions:
|
||||
* exp_window and gauss_window are S/N enhancing,
|
||||
* dexp_window and traf_window are resolution enhancing
|
||||
* standard windows [hamming, hanning, bartlett, blackman, kaiser-bessel]
|
||||
are also available
|
||||
self.x = time points
|
||||
elf.aquisition_time = aquisition time (no. samples / sampling_rate)
|
||||
line_broadening = line broadening factor (standard = 10 Hz)
|
||||
gaussian_multiplicator = Gaussian Multiplication Factor for
|
||||
the double exponential apodization
|
||||
function (standard = 0.3)
|
||||
"""
|
||||
def exp_window(self, line_broadening=10):
|
||||
"""
|
||||
exponential window
|
||||
"""
|
||||
"""
|
||||
exponential window
|
||||
"""
|
||||
apod = numpy.exp(-self.x*numpy.pi*line_broadening)
|
||||
for i in range(2):
|
||||
self.y[i] = self.y[i]*apod
|
||||
return self
|
||||
for i in range(2):
|
||||
self.y[i] = self.y[i]*apod
|
||||
return self
|
||||
|
||||
def gauss_window(self, line_broadening=10):
|
||||
apod = numpy.exp(-(self.x*line_broadening)**2)
|
||||
for i in range(2):
|
||||
self.y[i] = self.y[i]*apod
|
||||
return self
|
||||
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
|
||||
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
|
||||
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):
|
||||
apod = numpy.hanning(self.x.size)
|
||||
for i in range(2):
|
||||
self.y[i] = self.y[i]*apod
|
||||
return self
|
||||
apod = numpy.hanning(self.x.size)
|
||||
for i in range(2):
|
||||
self.y[i] = self.y[i]*apod
|
||||
return self
|
||||
|
||||
def hamming_window(self):
|
||||
apod = numpy.hamming(self.x.size)
|
||||
for i in range(2):
|
||||
self.y[i] = self.y[i]*apod
|
||||
return self
|
||||
apod = numpy.hamming(self.x.size)
|
||||
for i in range(2):
|
||||
self.y[i] = self.y[i]*apod
|
||||
return self
|
||||
|
||||
def blackman_window(self):
|
||||
apod = numpy.blackman(self.x.size)
|
||||
for i in range(2):
|
||||
self.y[i] = self.y[i]*apod
|
||||
return self
|
||||
apod = numpy.blackman(self.x.size)
|
||||
for i in range(2):
|
||||
self.y[i] = self.y[i]*apod
|
||||
return self
|
||||
|
||||
def bartlett_window(self):
|
||||
apod = numpy.bartlett(self.x.size)
|
||||
for i in range(2):
|
||||
self.y[i] = self.y[i]*apod
|
||||
return self
|
||||
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):
|
||||
if use_scipy == None:
|
||||
# modified Bessel function of zero kind order from somewhere
|
||||
def I_0(x):
|
||||
i0=0
|
||||
fac = lambda n:reduce(lambda a,b:a*(b+1),range(n),1)
|
||||
for n in xrange(20):
|
||||
i0 += ((x/2.0)**n/(fac(n)))**2
|
||||
return i0
|
||||
if use_scipy == None:
|
||||
# modified Bessel function of zero kind order from somewhere
|
||||
def I_0(x):
|
||||
i0=0
|
||||
fac = lambda n:reduce(lambda a,b:a*(b+1),range(n),1)
|
||||
for n in xrange(20):
|
||||
i0 += ((x/2.0)**n/(fac(n)))**2
|
||||
return i0
|
||||
|
||||
t = numpy.arange(self.x.size, type=numpy.Float) - self.x.size/2.0
|
||||
T = self.x.size
|
||||
# this is the window function array
|
||||
apod = I_0(beta*numpy.sqrt(1-(2*t/T)**2))/I_0(beta)
|
||||
else:
|
||||
# alternative method using scipy
|
||||
import scipy
|
||||
apod=scipy.kaiser(self.x.size, beta)
|
||||
t = numpy.arange(self.x.size, type=numpy.Float) - self.x.size/2.0
|
||||
T = self.x.size
|
||||
# this is the window function array
|
||||
apod = I_0(beta*numpy.sqrt(1-(2*t/T)**2))/I_0(beta)
|
||||
else:
|
||||
# alternative method using scipy
|
||||
import scipy
|
||||
apod=scipy.kaiser(self.x.size, beta)
|
||||
|
||||
for i in range(2):
|
||||
self.y[i] = self.y[i]*apod
|
||||
return self
|
||||
for i in range(2):
|
||||
self.y[i] = self.y[i]*apod
|
||||
return self
|
||||
|
||||
def autophase(self):
|
||||
"""
|
||||
@ -162,35 +162,35 @@ class DamarisFFT:
|
||||
return self
|
||||
|
||||
def fft(self, samples=None):
|
||||
"""
|
||||
Fouriertransform the timesignal inplace.
|
||||
For "zerofilling" set "samples" to a value higher than your data length.
|
||||
Shorten "samples" to truncate your data.
|
||||
samples takes only integer values
|
||||
"""
|
||||
Fouriertransform the timesignal inplace.
|
||||
For "zerofilling" set "samples" to a value higher than your data length.
|
||||
Shorten "samples" to truncate your data.
|
||||
samples takes only integer values
|
||||
|
||||
"""
|
||||
# Is this smart performance wise? Should I create an empty object?
|
||||
# Tests showed that this try except block performed 3.78ms
|
||||
# timesignal.baseline().fft()
|
||||
# with out this it needed 4.41 ms, thus this is justified :-)
|
||||
#try:
|
||||
# if self.was_copied:
|
||||
# new_object = self
|
||||
#except:
|
||||
# new_object = self+0
|
||||
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
|
||||
|
||||
"""
|
||||
# Is this smart performance wise? Should I create an empty object?
|
||||
# Tests showed that this try except block performed 3.78ms
|
||||
# timesignal.baseline().fft()
|
||||
# with out this it needed 4.41 ms, thus this is justified :-)
|
||||
#try:
|
||||
# if self.was_copied:
|
||||
# new_object = self
|
||||
#except:
|
||||
# new_object = self+0
|
||||
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):
|
||||
# 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
|
||||
# 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
|
||||
|
||||
@ -155,10 +155,10 @@ class DataPool(UserDict.DictMixin):
|
||||
complevel=complevel)
|
||||
except Exception,e:
|
||||
print "failed to write data_pool[\"%s\"]: %s"%(key,str(e))
|
||||
traceback_file=StringIO.StringIO()
|
||||
traceback.print_tb(sys.exc_info()[2], None, traceback_file)
|
||||
print "detailed traceback: %s\n"%str(e)+traceback_file.getvalue()
|
||||
traceback_file=None
|
||||
traceback_file=StringIO.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
|
||||
|
||||
@ -20,7 +20,7 @@ class AccumulatedValue:
|
||||
can be initialized by:
|
||||
No argument: no entries
|
||||
one argument: first entry
|
||||
two arguments: mean and its error, n is set 2
|
||||
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:
|
||||
@ -31,11 +31,11 @@ class AccumulatedValue:
|
||||
self.y=float(mean)
|
||||
self.y2=self.y**2
|
||||
self.n=1
|
||||
elif mean_err is None:
|
||||
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:
|
||||
elif n is None:
|
||||
self.n=2
|
||||
self.y=float(mean)*2
|
||||
self.y2=(float(mean_err)**2+float(mean)**2)*2
|
||||
@ -169,7 +169,7 @@ class MeasurementResult(Drawable.Drawable, UserDict.UserDict):
|
||||
sorted array of all dictionary entries without Accumulated Value objects with n==0
|
||||
"""
|
||||
keys=numpy.array(filter(lambda k: not (isinstance(self.data[k], AccumulatedValue) and self.data[k].n==0), self.data.keys()),
|
||||
dtype="Float64")
|
||||
dtype="Float64")
|
||||
keys.sort()
|
||||
return keys
|
||||
|
||||
|
||||
@ -1,29 +1,29 @@
|
||||
class Persistance :
|
||||
def __init__(self, shots):
|
||||
self.shots = shots
|
||||
self.accu = 0
|
||||
self.counter = 0
|
||||
self.result_list = []
|
||||
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
|
||||
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
|
||||
|
||||
@ -2,62 +2,62 @@
|
||||
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))
|
||||
"""
|
||||
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
|
||||
# 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
|
||||
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
|
||||
|
||||
@ -150,7 +150,7 @@ class BackendDriver(threading.Thread):
|
||||
# read state file
|
||||
statefile=file(self.statefilename,"r")
|
||||
statelines=statefile.readlines()
|
||||
statefile=None
|
||||
statefile=None
|
||||
statelinepattern=re.compile("<state name=\"([^\"]+)\" pid=\"([^\"]+)\" starttime=\"([^\"]+)\">")
|
||||
self.core_pid=-1
|
||||
for l in statelines:
|
||||
|
||||
@ -28,8 +28,8 @@ class ExperimentHandling(threading.Thread):
|
||||
|
||||
def run(self):
|
||||
dataspace={}
|
||||
exp_classes = __import__('damaris.experiments', dataspace, dataspace, ['Experiment'])
|
||||
for name in dir(exp_classes):
|
||||
exp_classes = __import__('damaris.experiments', dataspace, dataspace, ['Experiment'])
|
||||
for name in dir(exp_classes):
|
||||
if name[:2]=="__" and name[-2:]=="__": continue
|
||||
dataspace[name]=exp_classes.__dict__[name]
|
||||
del exp_classes
|
||||
|
||||
@ -24,18 +24,18 @@ class ExperimentWriter:
|
||||
self.inform_last_job=None
|
||||
job.job_id=self.no
|
||||
job_filename=os.path.join(self.spool,self.job_pattern%self.no)
|
||||
f=file(job_filename+".tmp","w")
|
||||
f.write(job.write_xml_string())
|
||||
f=file(job_filename+".tmp","w")
|
||||
f.write(job.write_xml_string())
|
||||
f.flush()
|
||||
f.close() # explicit close under windows necessary (don't know why)
|
||||
del f
|
||||
# this implementation tries to satisfiy msvc filehandle caching
|
||||
os.rename(job_filename+".tmp", job_filename)
|
||||
f.close() # explicit close under windows necessary (don't know why)
|
||||
del f
|
||||
# this implementation tries to satisfiy msvc filehandle caching
|
||||
os.rename(job_filename+".tmp", job_filename)
|
||||
#shutil.copyfile(job_filename+".tmp", job_filename)
|
||||
#try:
|
||||
# os.unlink(job_filename+".tmp")
|
||||
#except OSError:
|
||||
# print "could not delete temporary file %s.tmp"%job_filename
|
||||
#try:
|
||||
# os.unlink(job_filename+".tmp")
|
||||
#except OSError:
|
||||
# print "could not delete temporary file %s.tmp"%job_filename
|
||||
|
||||
self.no+=1
|
||||
|
||||
|
||||
@ -23,8 +23,8 @@ class ResultHandling(threading.Thread):
|
||||
def run(self):
|
||||
# execute it
|
||||
dataspace={}
|
||||
data_classes = __import__('damaris.data', dataspace, dataspace, ['*'])
|
||||
for name in dir(data_classes):
|
||||
data_classes = __import__('damaris.data', dataspace, dataspace, ['*'])
|
||||
for name in dir(data_classes):
|
||||
if name[:2]=="__" and name[-2:]=="__": continue
|
||||
dataspace[name]=data_classes.__dict__[name]
|
||||
del data_classes
|
||||
|
||||
Loading…
Reference in New Issue
Block a user