IPKM/python-damaris
3
0
Fork 0
Code Issues 6 Pull Requests Actions Packages Activity

fixed taps/space, update documentation

Browse Source
This commit is contained in:
Markus Rosenstihl 2014-11-26 15:46:27 +00:00
parent 7aeb010a83
commit d656daef5b
1 changed files with 147 additions and 130 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

277
src/data/DamarisFFT.py
View File

@ -2,55 +2,58 @@ import numpy
import sys import sys
import autophase import autophase
class DamarisFFT:
def clip(self, start=None, stop=None):
"""
:param float start: beginning of clipping
:param float stop: end of clipping
Method for clipping data, returns only the timesignal class DamarisFFT:
between start and stop """
is returned. 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. start and stop can be either time or frequency.
The unit is automatically determined 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 # check if start/stop order is properly
if start > stop: if start > stop:
# I could swap start/stop actually start, stop = stop, start
# TODO swap values?
raise # do nothing if one uses clip as a "placeholder"
# if one uses clip as a "placeholder"
if start == None and stop == None: if start == None and stop == None:
return self return self
if start == None: if start == None:
start = 0 start = self.x[ 0 ]
if stop == None: if stop == None:
stop = -1 stop = self.x[ -1 ]
# check if data is fft which changes the start/stop units # 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": if self.xlabel == "Frequency / Hz":
isfft = True isfft = True
start = self.x.size*(0.5 + start/self.sampling_rate) start = self.x.size * (0.5 + start / self.sampling_rate)
stop = self.x.size*(0.5 + stop/self.sampling_rate) stop = self.x.size * (0.5 + stop / self.sampling_rate)
else: else:
isfft = False isfft = False
# get the corresponding indices # get the corresponding indices
start *= self.sampling_rate start *= self.sampling_rate
stop *= self.sampling_rate stop *= self.sampling_rate
# check if boundaries make sense, raise exception otherwise # check if boundaries make sense, raise exception otherwise
if numpy.abs(int(start)-int(stop))<=0: if numpy.abs( int( start ) - int( stop ) ) <= 0:
raise ValueError("start stop too close: There are no values in the given boundaries!") raise ValueError( "start stop too close: There are no samples in the given boundaries!" )
for ch in xrange(len(self.y)):
# clip the data for each channel # clip the data for each channel
# TODO multi records for ch in xrange( len( self.y ) ):
self.y[ch] = self.y[ch][int(start):int(stop)] 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( start ):int( stop ) ]
# self.x = self.x[:int(stop)-int(start)]
self.x = self.x[int(start):int(stop)]
return self return self
def baseline(self, last_part=0.1): def baseline( self, last_part=0.1 ):
""" """
Correct the baseline of your data by subtracting the mean of the Correct the baseline of your data by subtracting the mean of the
last_part fraction of your data. last_part fraction of your data.
@ -59,152 +62,166 @@ class DamarisFFT:
last_part defaults to 0.1, i.e. last 10% of your data last_part defaults to 0.1, i.e. last 10% of your data
""" """
# TODO baselinecorrection for spectra after: # TODO baseline correction for spectra after:
# Heuer, A; Haeberlen, U.: J. Mag. Res.(1989) 85, Is 1, 79-94 # Heuer, A; Haeberlen, U.: J. Mag. Res.(1989) 85, Is 1, 79-94
# Should I create an empty object? n = int( self.x.size * last_part )
# I deided to do NOT a copy, but for ch in xrange( len( self.y ) ):
# rather modify the object self.y[ ch ] -= self.y[ ch ][ -n: ].mean( )
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 return self
def exp_window( self, line_broadening=10 ):
"""
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 function Exponential window function
:param float line_broadening: Applies apodization to time signal :param float line_broadening: default 10, line broadening factor in Hz
.. math:: .. math::
\\exp\\left(-\\pi\\cdot \\textsf{line_broadening} \\cdot t\\right) \\exp\\left(-\\pi\\cdot \\textsf{line_broadening} \\cdot t\\right)
""" """
apod = numpy.exp(-self.x*numpy.pi*line_broadening) apod = numpy.exp( -self.x * numpy.pi * line_broadening )
for i in range(2): for i in range( 2 ):
self.y[i] = self.y[i]*apod self.y[ i ] = self.y[ i ] * apod
return self return self
def gauss_window(self, line_broadening=10): def gauss_window( self, line_broadening=10 ):
apod = numpy.exp(-(self.x*line_broadening)**2) """
for i in range(2): Gaussian window function
self.y[i] = self.y[i]*apod
: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 return self
def dexp_window(self, line_broadening=-10, gaussian_multiplicator=0.3): def dexp_window( self, line_broadening=-10, gaussian_multiplicator=0.3 ):
apod = numpy.exp(-(self.x*line_broadening - gaussian_multiplicator*self.x.max())**2) apod = numpy.exp( -(self.x * line_broadening - gaussian_multiplicator * self.x.max( )) ** 2 )
for i in range(2): for i in range( 2 ):
self.y[i] = self.y[i]*apod self.y[ i ] = self.y[ i ] * apod
return self return self
def traf_window(self, line_broadening=10): def traf_window( self, line_broadening=10 ):
apod = (numpy.exp(-self.x*line_broadening))**2 / ( (numpy.exp(-self.x*line_broadening))**3 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): numpy.exp( -self.x.max( ) * line_broadening )) ** 3 )
self.y[i] = self.y[i]*apod for i in range( 2 ):
self.y[ i ] = self.y[ i ] * apod
return self return self
def hanning_window(self): def hanning_window( self ):
apod = numpy.hanning(self.x.size) """
for i in range(2): Symmetric centered window (hanning)
self.y[i] = self.y[i]*apod """
apod = numpy.hanning( self.x.size )
for i in range( 2 ):
self.y[ i ] = self.y[ i ] * apod
return self return self
def hamming_window(self): def hamming_window( self ):
apod = numpy.hamming(self.x.size) """
for i in range(2): Symmetric centered window (hamming)
self.y[i] = self.y[i]*apod """
apod = numpy.hamming( self.x.size )
for i in range( 2 ):
self.y[ i ] = self.y[ i ] * apod
return self return self
def blackman_window(self): def blackman_window( self ):
apod = numpy.blackman(self.x.size) """
for i in range(2): Symmetric centered window (blackmann)
self.y[i] = self.y[i]*apod """
apod = numpy.blackman( self.x.size )
for i in range( 2 ):
self.y[ i ] = self.y[ i ] * apod
return self return self
def bartlett_window(self): def bartlett_window( self ):
apod = numpy.bartlett(self.x.size) """
for i in range(2): Symmetric centered window (bartlett)
self.y[i] = self.y[i]*apod """
apod = numpy.bartlett( self.x.size )
for i in range( 2 ):
self.y[ i ] = self.y[ i ] * apod
return self return self
def kaiser_window(self, beta=4, use_scipy=None): def kaiser_window( self, beta=4, use_scipy=None ):
"""
Symmetric centered window (kaiser)
"""
if use_scipy == None: if use_scipy == None:
# modified Bessel function of zero kind order from somewhere # modified Bessel function of zero kind order from somewhere
def I_0(x): def I_0( x ):
i0=0 i0 = 0
fac = lambda n:reduce(lambda a,b:a*(b+1),range(n),1) fac = lambda n: reduce( lambda a, b: a * (b + 1), range( n ), 1 )
for n in xrange(20): for n in xrange( 20 ):
i0 += ((x/2.0)**n/(fac(n)))**2 i0 += ((x / 2.0) ** n / (fac( n ))) ** 2
return i0 return i0
t = numpy.arange(self.x.size, type=numpy.Float) - self.x.size/2.0 t = numpy.arange( self.x.size, type=numpy.Float ) - self.x.size / 2.0
T = self.x.size T = self.x.size
# this is the window function array # this is the window function array
apod = I_0(beta*numpy.sqrt(1-(2*t/T)**2))/I_0(beta) apod = I_0( beta * numpy.sqrt( 1 - (2 * t / T) ** 2 ) ) / I_0( beta )
else: else:
# alternative method using scipy # alternative method using scipy
import scipy import scipy
apod=scipy.kaiser(self.x.size, beta)
for i in range(2): apod = scipy.kaiser( self.x.size, beta )
self.y[i] = self.y[i]*apod
return self
def autophase(self): for i in range( 2 ):
""" self.y[ i ] = self.y[ i ] * apod
works nice with a SNR above 20 dB
10 V signal height to 1V noise width
"""
autophase.get_phase(self)
return self return self
def fft(self, samples=None): def autophase( self ):
""" """
Fouriertransform the timesignal inplace. Automatically phases the data to maximize real part.
For "zerofilling" set "samples" to a value higher than your data length.
Shorten "samples" to truncate your data. Works nice with a SNR above 20 dB, i.e.
samples takes only integer values 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
""" """
# Is this smart performance wise? Should I create an empty object? fft_of_signal = numpy.fft.fft( self.y[ 0 ] + 1j * self.y[ 1 ], n=samples )
# Tests showed that this try except block performed 3.78ms fft_of_signal = numpy.fft.fftshift( fft_of_signal )
# timesignal.baseline().fft() dwell = 1.0 / self.sampling_rate
# 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 n = fft_of_signal.size
fft_frequencies = numpy.fft.fftfreq(n, dwell) fft_frequencies = numpy.fft.fftfreq( n, dwell )
self.x = numpy.fft.fftshift(fft_frequencies) self.x = numpy.fft.fftshift( fft_frequencies )
self.y[0] = fft_of_signal.real self.y[ 0 ] = fft_of_signal.real
self.y[1] = fft_of_signal.imag self.y[ 1 ] = fft_of_signal.imag
self.set_xlabel("Frequency / Hz") self.set_xlabel( "Frequency / Hz" )
return self return self
def magnitude(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 # 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[ 0 ] = numpy.sqrt( self.y[ 0 ] ** 2 + self.y[ 1 ] ** 2 )
self.y[1] *= 0 #self.y[0].copy() self.y[ 1 ] *= 0 # self.y[0].copy()
return self return self