PHS3000 PET

This documents the monashspa.PHS3000.PET library that you will import into code used in the PHS3000 unit when performing experiment 4.6 PET.

monashspa.PHS3000.PET.fft(x, n=None, axis=-1, overwrite_x=False)[source]

Return discrete Fourier transform of real or complex sequence.

The returned complex array contains y(0), y(1),..., y(n-1), where

y(j) = (x * exp(-2*pi*sqrt(-1)*j*np.arange(n)/n)).sum().

Parameters:

  • x (array_like) – Array to Fourier transform.

  • n (int, optional) – Length of the Fourier transform. If n < x.shape[axis], x is truncated. If n > x.shape[axis], x is zero-padded. The default results in n = x.shape[axis].

  • axis (int, optional) – Axis along which the fft’s are computed; the default is over the last axis (i.e., axis=-1).

  • overwrite_x (bool, optional) – If True, the contents of x can be destroyed; the default is False.

Returns:

z –

with the elements:

[y(0),y(1),..,y(n/2),y(1-n/2),...,y(-1)]        if n is even
[y(0),y(1),..,y((n-1)/2),y(-(n-1)/2),...,y(-1)]  if n is odd

where:

y(j) = sum[k=0..n-1] x[k] * exp(-sqrt(-1)*j*k* 2*pi/n), j = 0..n-1

Return type:

complex ndarray

See also

ifft

Inverse FFT

rfft

FFT of a real sequence

Notes

The packing of the result is “standard”: If A = fft(a, n), then A[0] contains the zero-frequency term, A[1:n/2] contains the positive-frequency terms, and A[n/2:] contains the negative-frequency terms, in order of decreasingly negative frequency. So ,for an 8-point transform, the frequencies of the result are [0, 1, 2, 3, -4, -3, -2, -1]. To rearrange the fft output so that the zero-frequency component is centered, like [-4, -3, -2, -1, 0, 1, 2, 3], use fftshift.

Both single and double precision routines are implemented. Half precision inputs will be converted to single precision. Non-floating-point inputs will be converted to double precision. Long-double precision inputs are not supported.

This function is most efficient when n is a power of two, and least efficient when n is prime.

Note that if x is real-valued, then A[j] == A[n-j].conjugate(). If x is real-valued and n is even, then A[n/2] is real.

If the data type of x is real, a “real FFT” algorithm is automatically used, which roughly halves the computation time. To increase efficiency a little further, use rfft, which does the same calculation, but only outputs half of the symmetrical spectrum. If the data is both real and symmetrical, the dct can again double the efficiency by generating half of the spectrum from half of the signal.

Examples

>>> import numpy as np
>>> from scipy.fftpack import fft, ifft
>>> x = np.arange(5)
>>> np.allclose(fft(ifft(x)), x, atol=1e-15)  # within numerical accuracy.
True
monashspa.PHS3000.PET.fftfreq(n, d=1.0, device=None)

Return the Discrete Fourier Transform sample frequencies.

The returned float array f contains the frequency bin centers in cycles per unit of the sample spacing (with zero at the start). For instance, if the sample spacing is in seconds, then the frequency unit is cycles/second.

Given a window length n and a sample spacing d:

f = [0, 1, ...,   n/2-1,     -n/2, ..., -1] / (d*n)   if n is even
f = [0, 1, ..., (n-1)/2, -(n-1)/2, ..., -1] / (d*n)   if n is odd
Parameters:

  • n (int) – Window length.

  • d (scalar, optional) – Sample spacing (inverse of the sampling rate). Defaults to 1.

  • device (str, optional) –

    The device on which to place the created array. Default: None. For Array-API interoperability only, so must be "cpu" if passed.

    Added in version 2.0.0.

Returns:

f – Array of length n containing the sample frequencies.

Return type:

ndarray

Examples

>>> import numpy as np
>>> signal = np.array([-2, 8, 6, 4, 1, 0, 3, 5], dtype=float)
>>> fourier = np.fft.fft(signal)
>>> n = signal.size
>>> timestep = 0.1
>>> freq = np.fft.fftfreq(n, d=timestep)
>>> freq
array([ 0.  ,  1.25,  2.5 , ..., -3.75, -2.5 , -1.25])
monashspa.PHS3000.PET.ifft(x, n=None, axis=-1, overwrite_x=False)[source]

Return discrete inverse Fourier transform of real or complex sequence.

The returned complex array contains y(0), y(1),..., y(n-1), where

y(j) = (x * exp(2*pi*sqrt(-1)*j*np.arange(n)/n)).mean().

Parameters:

  • x (array_like) – Transformed data to invert.

  • n (int, optional) – Length of the inverse Fourier transform. If n < x.shape[axis], x is truncated. If n > x.shape[axis], x is zero-padded. The default results in n = x.shape[axis].

  • axis (int, optional) – Axis along which the ifft’s are computed; the default is over the last axis (i.e., axis=-1).

  • overwrite_x (bool, optional) – If True, the contents of x can be destroyed; the default is False.

Returns:

ifft – The inverse discrete Fourier transform.

Return type:

ndarray of floats

See also

fft

Forward FFT

Notes

Both single and double precision routines are implemented. Half precision inputs will be converted to single precision. Non-floating-point inputs will be converted to double precision. Long-double precision inputs are not supported.

This function is most efficient when n is a power of two, and least efficient when n is prime.

If the data type of x is real, a “real IFFT” algorithm is automatically used, which roughly halves the computation time.

Examples

>>> from scipy.fftpack import fft, ifft
>>> import numpy as np
>>> x = np.arange(5)
>>> np.allclose(ifft(fft(x)), x, atol=1e-15)  # within numerical accuracy.
True
class monashspa.PHS3000.PET.interp1d(x, y, kind='linear', axis=-1, copy=True, bounds_error=None, fill_value=nan, assume_sorted=False)[source]

Interpolate a 1-D function (legacy).

x and y are arrays of values used to approximate some function f: y = f(x). This class returns a function whose call method uses interpolation to find the value of new points.

Parameters:

  • x ((npoints, ) array_like) – A 1-D array of real values.

  • y ((..., npoints, ...) array_like) – An N-D array of real values. The length of y along the interpolation axis must be equal to the length of x. Use the axis parameter to select correct axis. Unlike other interpolators, the default interpolation axis is the last axis of y.

  • kind (str or int, optional) – Specifies the kind of interpolation as a string or as an integer specifying the order of the spline interpolator to use. The string has to be one of ‘linear’, ‘nearest’, ‘nearest-up’, ‘zero’, ‘slinear’, ‘quadratic’, ‘cubic’, ‘previous’, or ‘next’. ‘zero’, ‘slinear’, ‘quadratic’ and ‘cubic’ refer to a spline interpolation of zeroth, first, second or third order; ‘previous’ and ‘next’ simply return the previous or next value of the point; ‘nearest-up’ and ‘nearest’ differ when interpolating half-integers (e.g. 0.5, 1.5) in that ‘nearest-up’ rounds up and ‘nearest’ rounds down. Default is ‘linear’.

  • axis (int, optional) – Axis in the y array corresponding to the x-coordinate values. Unlike other interpolators, defaults to axis=-1.

  • copy (bool, optional) – If True, the class makes internal copies of x and y. If False, references to x and y are used if possible. The default is to copy.

  • bounds_error (bool, optional) – If True, a ValueError is raised any time interpolation is attempted on a value outside of the range of x (where extrapolation is necessary). If False, out of bounds values are assigned fill_value. By default, an error is raised unless fill_value="extrapolate".

  • fill_value (array-like or (array-like, array_like) or "extrapolate", optional) –

    • if a ndarray (or float), this value will be used to fill in for requested points outside of the data range. If not provided, then the default is NaN. The array-like must broadcast properly to the dimensions of the non-interpolation axes.

    • If a two-element tuple, then the first element is used as a fill value for x_new < x[0] and the second element is used for x_new > x[-1]. Anything that is not a 2-element tuple (e.g., list or ndarray, regardless of shape) is taken to be a single array-like argument meant to be used for both bounds as below, above = fill_value, fill_value. Using a two-element tuple or ndarray requires bounds_error=False.

      Added in version 0.17.0.

    • If “extrapolate”, then points outside the data range will be extrapolated.

      Added in version 0.17.0.

  • assume_sorted (bool, optional) – If False, values of x can be in any order and they are sorted first. If True, x has to be an array of monotonically increasing values.

fill_value
__call__()

See also

splrep(), splev()

UnivariateSpline()

An object-oriented wrapper of the FITPACK routines.

interp2d()

2-D interpolation

Notes

Calling interp1d with NaNs present in input values results in undefined behaviour.

Input values x and y must be convertible to float values like int or float.

If the values in x are not unique, the resulting behavior is undefined and specific to the choice of kind, i.e., changing kind will change the behavior for duplicates.

Array API Standard Support

interp1d is not in-scope for support of Python Array API Standard compatible backends other than NumPy.

See dev-arrayapi for more information.

Examples

>>> import numpy as np
>>> import matplotlib.pyplot as plt
>>> from scipy import interpolate
>>> x = np.arange(0, 10)
>>> y = np.exp(-x/3.0)
>>> f = interpolate.interp1d(x, y)

>>> xnew = np.arange(0, 9, 0.1)
>>> ynew = f(xnew)   # use interpolation function returned by `interp1d`
>>> plt.plot(x, y, 'o', xnew, ynew, '-')
>>> plt.show()
property fill_value

The fill value.

monashspa.PHS3000.PET.pet_rebuild(filepath, filter_name=None, npoints=None, call_show=True)[source]

Perform inverse radon transform on acquired PET data and plots results

Parameters:

filepath – A string containing the path to the txt file containing the PET data

Keyword Arguments:

  • filter_name –

    A string containing the name of the filter to use during reconstruction. Defaults to None (no reconstruction). Options are:

    • None: don’t reconstruct

    • ’none’: Reconstruct with no filter

    • ’ramp’: Reconstruct using the Ram-Lak filter

    • ’Shepp-Logan’: Reconstruct using the Shepp-Logan filter

    • ’cosine’: Reconstruct using the cosine filter

    • ’hamming’: Reconstruct using the hamming filter

    • ’hann’: Reconstruct using the hann filter

  • npoints – The number of points to reconstruct

  • call_show – Whether to call matplotlib.pyplot.show() at the end of the function (prior to returning). Defaults to True. Set this to False if you are not using Spyder/IPython and wish your entire script to complete before showing any plots. Note, you will need to explicitly call matplotlib.pyplot.show() if you set this to False.

Returns:

A 2D numpy array containing the coincidence counts (rows correspond to each linear stage position and columns to each rotation stage position)

monashspa.PHS3000.PET.warn(message, category=None, stacklevel=1, source=None)

Issue a warning, or maybe ignore it or raise an exception.