Commit 86e5e80c authored by eric pellegrini's avatar eric pellegrini
Browse files

Bug fix when computing FFT using odd number of points

parent f6aa12ab
......@@ -50,6 +50,6 @@ class GaussianInstrumentResolution(IInstrumentResolution):
mu = self._configuration["mu"]["value"]
sigma = self._configuration["sigma"]["value"]
self._frequencyWindow = (1.0/(sigma*numpy.sqrt(2.0*numpy.pi)))*numpy.exp(-0.5*((frequencies-mu)/sigma)**2)
self._timeWindow = numpy.fft.fftshift(numpy.abs(numpy.fft.ifft(self._frequencyWindow))/dt)
\ No newline at end of file
self._frequencyWindow = (numpy.sqrt(2.0*numpy.pi)/sigma)*numpy.exp(-0.5*((frequencies-mu)/sigma)**2)
self._timeWindow = numpy.fft.fftshift(numpy.fft.ifft(numpy.fft.ifftshift(self._frequencyWindow))/dt)
......@@ -51,9 +51,7 @@ class LorentzianInstrumentResolution(IInstrumentResolution):
mu = self._configuration["mu"]["value"]
sigma = self._configuration["sigma"]["value"]
fact = 0.5*sigma
self._frequencyWindow = (1.0/numpy.pi)*(fact/((frequencies-mu)**2 + fact**2))
self._timeWindow = numpy.fft.fftshift(numpy.abs(numpy.fft.ifft(self._frequencyWindow))/dt)
\ No newline at end of file
self._frequencyWindow = (2.0*sigma)/((frequencies-mu)**2 + sigma**2)
self._timeWindow = numpy.fft.fftshift(numpy.fft.ifft(numpy.fft.ifftshift(self._frequencyWindow))/dt)
\ No newline at end of file
......@@ -57,12 +57,10 @@ class PseudoVoigtInstrumentResolution(IInstrumentResolution):
muG = self._configuration["mu_gaussian"]["value"]
sigmaG = self._configuration["sigma_gaussian"]["value"]
gContribution = (1.0/(sigmaG*numpy.sqrt(2.0*numpy.pi)))*numpy.exp(-0.5*((frequencies-muG)/sigmaG)**2)
fact = 0.5*sigmaL
gaussian = (numpy.sqrt(2.0*numpy.pi)/sigmaG)*numpy.exp(-0.5*((frequencies-muG)/sigmaG)**2)
lContribution = (1.0/numpy.pi)*(fact/((frequencies-muL)**2 + fact**2))
lorentzian = (2.0*sigmaL)/((frequencies-muL)**2 + sigmaL**2)
self._frequencyWindow = eta*lContribution + (1.0-eta)*gContribution
self._timeWindow = numpy.fft.fftshift(numpy.abs(numpy.fft.ifft(self._frequencyWindow))/dt)
self._frequencyWindow = eta*lorentzian + (1.0-eta)*gaussian
self._timeWindow = numpy.fft.fftshift(numpy.fft.ifft(numpy.fft.ifftshift(self._frequencyWindow))/dt)
\ No newline at end of file
......@@ -51,6 +51,7 @@ class SquareInstrumentResolution(IInstrumentResolution):
mu = self._configuration["mu"]["value"]
sigma = self._configuration["sigma"]["value"]
self._frequencyWindow = numpy.where((numpy.abs(frequencies-mu)-sigma) > 0,0.0,1.0/(2.0*sigma))
self._frequencyWindow = 2.0*numpy.pi*numpy.where((numpy.abs(frequencies-mu)-sigma) > 0,0.0,1.0/(2.0*sigma))
self._timeWindow = numpy.fft.fftshift(numpy.abs(numpy.fft.ifft(self._frequencyWindow))/dt)
\ No newline at end of file
self._timeWindow = numpy.fft.fftshift(numpy.fft.ifft(numpy.fft.ifftshift(self._frequencyWindow))/dt)
\ No newline at end of file
......@@ -53,6 +53,6 @@ class TriangularInstrumentResolution(IInstrumentResolution):
val = numpy.abs(frequencies-mu) - sigma
self._frequencyWindow = numpy.where( val >= 0, 0.0, -val/sigma**2)
self._timeWindow = numpy.fft.fftshift(numpy.abs(numpy.fft.ifft(self._frequencyWindow))/dt)
\ No newline at end of file
self._frequencyWindow = 2.0*numpy.pi*numpy.where( val >= 0, 0.0, -val/sigma**2)
self._timeWindow = numpy.fft.fftshift(numpy.fft.ifft(numpy.fft.ifftshift(self._frequencyWindow))/dt)
......@@ -224,9 +224,18 @@ def get_spectrum(signal,window=None,timeStep=1.0,axis=0):
if window is None:
window = numpy.ones(signal.shape[axis])
window /= window[len(window)/2]
s = [numpy.newaxis]*signal.ndim
s[axis] = slice(None)
# We compute the unitary inverse fourier transform with angular frequencies as described in
# https://en.wikipedia.org/wiki/Fourier_transform#Discrete_Fourier_Transforms_and_Fast_Fourier_Transforms
# For information about the manipulation around fftshift and ifftshift
# http://www.mathworks.com/matlabcentral/newsreader/view_thread/285244
return numpy.fft.fftshift(numpy.abs(numpy.fft.fft(signal*window[s],axis=axis)),axes=axis)*timeStep
fftSignal = 0.5*numpy.fft.fftshift(numpy.fft.fft(numpy.fft.ifftshift(signal*window[s],axes=axis),axis=axis),axes=axis)*timeStep/numpy.pi
return fftSignal.real
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment