Reputation: 91
I've written some code that enables me to take data from DICOM files and separate the data into the individual FID signals.
from pylab import *
import dicom
import numpy as np
import scipy as sp
plan=dicom.read_file("1.3.46.670589.11.38085.5.22.3.1.4792.2013050818105496124")
all_points = array(plan.SpectroscopyData)
cmplx_data = all_points[0::2] -1j*all_points[1::2]
frames = int(plan.NumberOfFrames)
fid_pts = len(cmplx_data)/frames
del_t = plan.AcquisitionDuration / (frames * fid_pts)
fid_list = []
for fidN in arange(frames):
offset = fidN * fid_pts
current_fid = cmplx_data[offset:offset+fid_pts]
fid_list.append(current_fid)
I'd now like to quantify some of this data so, after applying a Fourier transform and shift I've tried to use Scipy's quad function and encountered the following error:
spec = fftshift(fft(fid_list[0])))
sp.integrate.quad(spec, 660.0, 700.0)
error Traceback (most recent call last)
/home/dominicc/Experiments/In Vitro/Glu, Cr/Phantom 1: 10mM Cr, 5mM Glu/WIP_SV_PRESS_ME_128TEs_5ms_spacing_1828/<ipython-input-107-17cb50e45927> in <module>()
----> 1 sp.integrate.quad(fid, 660.0, 700.0)
/usr/lib/python2.7/dist-packages/scipy/integrate/quadpack.pyc in quad(func, a, b, args, full_output, epsabs, epsrel, limit, points, weight, wvar, wopts, maxp1, limlst)
243 if type(args) != type(()): args = (args,)
244 if (weight is None):
--> 245 retval = _quad(func,a,b,args,full_output,epsabs,epsrel,limit,points)
246 else:
247 retval = _quad_weight(func,a,b,args,full_output,epsabs,epsrel,limlst,limit,maxp1,weight,wvar,wopts)
/usr/lib/python2.7/dist-packages/scipy/integrate/quadpack.pyc in _quad(func, a, b, args, full_output, epsabs, epsrel, limit, points)
307 if points is None:
308 if infbounds == 0:
--> 309 return _quadpack._qagse(func,a,b,args,full_output,epsabs,epsrel,limit)
310 else:
311 return _quadpack._qagie(func,bound,infbounds,args,full_output,epsabs,epsrel,limit)
error: First argument must be a callable function.
Could somebody please suggest a way to make this object callable? I'm still not particularly clear after reading What is a "callable" in Python?. After reading that I tried
def fid():
spec = fftshift(fft(fid_list[0]))
return spec
Which returned the same error.
Any help would be greatly appreciated, thanks.
Upvotes: 2
Views: 15581
Reputation: 1
As a beginner being interest in FFT of QM wave nature for years, I'd like to contribute a simple coding to the coming Python beginners who might be asking the same on Stack Exchange.
The error syntax seeing is related to how I am attempting to integrate the FFT_psik (i.e. an array of FFT results) over the frequency domain k. The issue here is that np.abs(FFT_psik)**2 is an array (not a callable scalar function
) , and when I am attempting to multiply an array by array inside the integrand, things get messy for the quad function mathematically.
I had ever encountered the same syntax problem, can be outlined here:
def integrand(k_value):
return k_value**2 * np.abs(FFT_psik)**2 # FFT_psik is an array of FFT of psi
FFT_psik_2nd_moment_quad, error_quad = quad(integrand, k.min(), k.max(), full_output=True)
print(f"FFT_psik_2nd_moment_quad: {FFT_psik_2nd_moment_quad}")
print(f"Estimated error from quad: {error_quad}")
To resolve it, one needs to interpolate np.abs(FFT_psik)**2 so that given any k value, quad can compute the scalar function value for any k values. The interp1d
function from scipy.interpolate is a convenient way to fix the problem.
First, make np.abs(FFT_psik)**2 callable scalar function
using interp1d
.
Next, rewrite my quad integrand as a scalar function
and use the interpolation to integrate with quad. Here's an example showing how it works fine with resolving this:
#!/usr/bin/python
import numpy as np
from numpy.fft import fft, fftshift
from scipy.integrate import simps, quad
from scipy.interpolate import interp1d
def psi_x(x, α):
A = (2*α/np.pi)**(1/4)
return A * np.exp(-α * x**2)
# Define x, k arrays and constant α
α = 3
x1 = -20 # Adjust based on the desired Δk size
x2 = 20 # Adjust based on the desired Δk size
N = 1000 # Increase for better accuracy in x-domain
x = np.linspace(x1, x2, N) # Adjust the range and points as needed
Δx = x[1] - x[0]
k = fftshift(np.fft.fftfreq(x.size, Δx)) # frequency grid array in K-space
Δk = k[1] - k[0] # the desired Δk size is smaller enough to catch the varying of FFT spectrum ROI
# Calculate the wave function in x-space and its FFT spectrum in k-space
psi = psi_x(x, α)
FFT_psik = fftshift(fft(psi))
### Calculate the 2nd moments using 'np.sum' method and 'simps' integration
# FFT_psik_2nd_moment = np.sum(k*k*np.abs(FFT_psik)**2) / np.sum(np.abs(FFT_psik)**2) # no good for my use case, error is unpredictable large
FFT_psik_2nd_moment = np.sum(k**2 * np.abs(FFT_psik)**2) * Δk # reliable one
FFT_psik_2nd_momentS = simps(k**2 * np.abs(FFT_psik)**2, k)
print("Second Moment Estimation Alternatives: \n")
print(f"FFT_psik_2nd_moment of np.sum(^2) : {FFT_psik_2nd_moment}")
print(f"Using 'simps' integration method : {FFT_psik_2nd_momentS} \n")
### Using Gaussian quadrature to compute the 2nd moment of FFT_psik. One needs to define inegrand as a scalar function (not array) before calling 'quad'
# 1. Interpolate |FFT_psik|^2 making it callable 'scalar function'
FFT_psik_inter = interp1d(k, np.abs(FFT_psik)**2, kind='cubic', fill_value="extrapolate")
# 2. Define the integrand for 'quad'
def integrand(k_val):
return k_val**2 * FFT_psik_inter(k_val)
# 3. Integrate using the saclar integrand() with quad
FFT_psik_2nd_moment_quad, error_quad = quad(integrand, k.min(), k.max())
print("Using Gaussian quad numerical integration: \n")
print(f"FFT_psik_2nd_moment by quad method: {FFT_psik_2nd_moment_quad}")
print(f"Estimated error from Gaussian quad: {error_quad}") # the side benefit of using quad is providing the reliable 'estimated error'
The actual outputs in my Python environment are as follows:
Second Moment Estimation Alternatives:
FFT_psik_2nd_moment of np.sum(^2) : 47.399363716988454
Using 'simps' integration method : 47.39936371698599
Using Gaussian quad numerical integration:
FFT_psik_2nd_moment by quad method: 47.39936294876803
Estimated error from Gaussian quad: 6.517835867177056e-07
By interpolating, one may effectively create a callable scalar function
which returns the squared magnitude of FFT_psik for arbitrary k values. This resolves the issues when using quad.
Realizing this might be arriving quite late for docar
who asked, being late for even after a decade, hope that this could be beneficial for newcomers to Python mathematics, as it has been for me.
Upvotes: 0
Reputation: 94595
Since you are integrating a function defined only on a grid (i.e., an array of numbers), you need to use one of the routines from "Integration, given fixed samples".
In fact, the first argument of quad()
must be a function, something you can call (a "callable"). For instance, you can do:
>>> from scipy.integrate import quad
>>> def f(x):
... return x**2
...
>>> quad(f, 0, 1)
(0.33333333333333337, 3.700743415417189e-15)
This can be done instead by sampling f()
:
>>> from scipy.integrate import simps
>>> x_grid = numpy.linspace(0, 1, 101) # 101 numbers between 0 and 1
>>> simps(x_grid**2, dx=x_grid[1]-x_grid[0]) # x_grid**2 is an *array*. dx=0.01 between x_grid values
0.33333333333333337
Your situation is like in this second example: you integrate a sampled function, so it is natural to use one or cumtrapz()
, simps()
or romb()
.
Upvotes: 5