#
# powerspectrumseries.py - Classes used by the PowerSpectrumPanel.
#
# Author: Paul McCarthy <pauldmccarthy@gmail.com>
#
"""This module provides :class:`.DataSeries` sub-classes which are used
by the :class:`.PowerSpectrumPanel` for plotting power spectra.
The following classes are provided:
.. autosummary::
:nosignatures:
PowerSpectrumSeries
VoxelPowerSpectrumSeries
ComplexPowerSpectrumSeries
ImaginaryPowerSpectrumSeries
MagnitudePowerSpectrumSeries
PhasePowerSpectrumSeries
MelodicPowerSpectrumSeries
"""
import logging
import numpy as np
import numpy.fft as fft
import fsl.data.image as fslimage
import fsl.data.melodicimage as fslmelimage
import fsleyes_props as props
import fsleyes.colourmaps as fslcm
import fsleyes.strings as strings
from . import dataseries
log = logging.getLogger(__name__)
[docs]def calcPowerSpectrum(data):
"""Calculates a power spectrum for the given one-dimensional data array.
:arg data: Numpy array containing the time series data
:returns: If ``data`` contains real values, the magnitude of the power
spectrum is returned. If ``data`` contains complex values,
the complex power spectrum is returned.
"""
# Fourier transform on complex data
if np.issubdtype(data.dtype, np.complexfloating):
data = fft.fft(data)
data = fft.fftshift(data)
# Fourier transform on real data - we
# calculate and return the magnitude.
# We also drop the first (zero-frequency)
# term (see the rfft docs) as it is
# kind of useless for display purposes
else:
data = fft.rfft(data)[1:]
data = magnitude(data)
return data
[docs]def calcFrequencies(nsamples, sampleTime, dtype):
"""Calculates the frequencies of the power spectrum for the given
data.
:arg nsamples: Number of samples in the input time series data
:arg sampleTime: Time between each data point
:arg dtype: Data type - the calculation differs depending on
whether the data is real or complex.
:returns: A ``numpy`` array containing the frequencies of the
power spectrum for ``data``
"""
if np.issubdtype(dtype, np.complexfloating):
xdata = fft.fftfreq(nsamples, sampleTime)
xdata = fft.fftshift(xdata)
else:
# Drop the zero-frequency term
# (see calcPowerSpectrum)
xdata = fft.rfftfreq(nsamples, sampleTime)[1:]
return xdata
[docs]def magnitude(data):
"""Returns the magnitude of the given complex data. """
real = data.real
imag = data.imag
return np.sqrt(real ** 2 + imag ** 2)
[docs]def phase(data):
"""Returns the phase of the given complex data. """
real = data.real
imag = data.imag
return np.arctan2(imag, real)
[docs]def normalise(data, dmin=None, dmax=None, nmin=-1, nmax=1):
"""Returns ``data``, rescaled to the range [nmin, nmax].
If dmin and dmax are provided, the data is normalised with respect to them,
rather than being normalised by the data minimum/maximum.
"""
if dmin is None: dmin = data.min()
if dmax is None: dmax = data.max()
return nmin + (nmax - nmin) * (data - dmin) / (dmax - dmin)
[docs]def phaseCorrection(spectrum, freqs, p0, p1):
"""Applies phase correction to the given complex power spectrum.
:arg spectrum: Complex-valued power spectrum
:arg freqs: Spectrum frequency bins
:arg p0: Zero order phase correction term
:arg p1: First order phase correction term
:returns: The corrected power spectrum.
"""
exp = 1j * 2 * np.pi * (p0 / 360 + freqs * p1)
return np.exp(exp) * spectrum
[docs]class PowerSpectrumSeries:
"""The ``PowerSpectrumSeries`` encapsulates a power spectrum data series
from an overlay. The ``PowerSpectrumSeries`` class is a base mixin class
for all other classes in this module.
"""
varNorm = props.Boolean(default=False)
"""If ``True``, the fourier-transformed data is normalised to the range
[0, 1] before plotting.
.. note:: The :class:`ComplexPowerSpectrumSeries` applies normalisation
differently.
"""
@property
def sampleTime(self):
"""Returns the time between time series samples for the overlay
data. """
if isinstance(self.overlay, fslmelimage.MelodicImage):
return self.overlay.tr
elif isinstance(self.overlay, fslimage.Image):
return self.overlay.pixdim[3]
else:
return 1
[docs]class VoxelPowerSpectrumSeries(dataseries.VoxelDataSeries,
PowerSpectrumSeries):
"""The ``VoxelPowerSpectrumSeries`` class encapsulates the power spectrum
of a single voxel from a 4D :class:`.Image` overlay. The voxel is dictated
by the :attr:`.DisplayContext.location` property.
"""
[docs] def __init__(self, *args, **kwargs):
"""Create a ``VoxelPowerSpectrumSeries``. All arguments are passed
to the :meth:`VoxelDataSeries.__init__` method. A :exc:`ValueError`
is raised if the overlay is not a 4D :class:`.Image`.
"""
dataseries.VoxelDataSeries.__init__(self, *args, **kwargs)
if self.overlay.ndim < 4:
raise ValueError('Overlay is not a 4D image')
[docs] def currentVoxelData(self, location):
"""Overrides :meth:`.VoxelDataSeries.currentVoxelData`. Retrieves
the data at the specified location, then performs a fourier transform
on it and returnes the result.
"""
data = dataseries.VoxelDataSeries.currentVoxelData(self, location)
data = calcPowerSpectrum(data)
return data
[docs] def getData(self):
"""Returns the ``(xdata, ydata)`` to be plotted from the current voxel
location.
"""
overlay = self.overlay
ydata = self.dataAtCurrentVoxel()
xdata = calcFrequencies(overlay.shape[3],
self.sampleTime,
overlay.dtype)
if self.varNorm:
ydata = normalise(ydata)
return xdata, ydata
[docs]class ComplexPowerSpectrumSeries(VoxelPowerSpectrumSeries):
"""This class is the frequency-spectrum equivalent of the
:class:`.ComplexTimeSeries` class - see it for more details.
"""
plotReal = props.Boolean(default=True)
plotImaginary = props.Boolean(default=False)
plotMagnitude = props.Boolean(default=False)
plotPhase = props.Boolean(default=False)
zeroOrderPhaseCorrection = props.Real(default=0)
"""Apply zero order phase correction to the power spectrum of the complex
data.
"""
firstOrderPhaseCorrection = props.Real(default=0)
"""Apply first order phase correction to the power spectrum of the complex
data.
"""
[docs] def __init__(self, overlay, overlayList, displayCtx, plotCanvas):
"""Create a ``ComplexPowerSpectrumSeries``. All arguments are
passed through to the :class:`VoxelPowerSpectrumSeries` constructor.
"""
VoxelPowerSpectrumSeries.__init__(
self, overlay, overlayList, displayCtx, plotCanvas)
# Separate DataSeries for the imaginary/
# magnitude/phase signals, returned by
# the extraSeries method
self.__imagps = ImaginaryPowerSpectrumSeries(
self, overlay, overlayList, displayCtx, plotCanvas)
self.__magps = MagnitudePowerSpectrumSeries(
self, overlay, overlayList, displayCtx, plotCanvas)
self.__phaseps = PhasePowerSpectrumSeries(
self, overlay, overlayList, displayCtx, plotCanvas)
for ps in (self.__imagps, self.__magps, self.__phaseps):
ps.colour = fslcm.randomDarkColour()
ps.bindProps('alpha', self)
ps.bindProps('lineWidth', self)
ps.bindProps('lineStyle', self)
[docs] def makeLabelBase(self):
"""Returns a string to be used as the label prefix for this
``ComplexPowerSpectrumSeries`` instance, and for the imaginary,
magnitude, and phase child series.
"""
return VoxelPowerSpectrumSeries.makeLabel(self)
[docs] def makeLabel(self):
"""Returns a label to use for this data series. """
return '{} ({})'.format(self.makeLabelBase(), strings.labels[self])
[docs] def getData(self, component='real'):
"""If :attr:`plotReal` is true, returns the real component of the power
spectrum of the data at the current voxel. Otherwise returns ``(None,
None)``.
Every time this method is called, the power spectrum is retrieved (see
the :class:`VoxelPowerSpectrumSeries` class), phase correction is
applied if set, andthe data is normalised, if set. A tuple containing
the ``(xdata, ydata)`` is returned, with ``ydata`` containing the
requested ``component`` ( ``'real'``, ``'imaginary'``,
``'magnitude'``, or ``'phase'``).
This method is called by the :class:`ImaginarySpectrumPowerSeries`,
:class:`MagnitudeSpectrumPowerSeries`, and
:class:`PhasePowerSpectrumPowerSeries` instances that are associated
with this data series.
"""
if ((component == 'real') and (not self.plotReal)) or \
((component == 'imaginary') and (not self.plotImaginary)) or \
((component == 'magnitude') and (not self.plotMagnitude)) or \
((component == 'phase') and (not self.plotPhase)):
return None, None
# See VoxelPowerSpectrumSeries - the data
# is already fourier-transformed
ydata = self.dataAtCurrentVoxel()
if ydata is None:
return None, None
# All of the calculations below are repeated
# for each real/imag/mag/phase series that
# gets plotted. But keeping the code together
# and clean is currently more important than
# performance, as there is not really any
# performance hit.
overlay = self.overlay
xdata = calcFrequencies(overlay.shape[3],
self.sampleTime,
overlay.dtype)
if self.zeroOrderPhaseCorrection != 0 or \
self.firstOrderPhaseCorrection != 0:
ydata = phaseCorrection(ydata,
xdata,
self.zeroOrderPhaseCorrection,
self.firstOrderPhaseCorrection)
# Normalise magnitude, real, imaginary
# components with respect to magnitude.
# Normalise phase independently.
if self.varNorm:
mag = magnitude(ydata)
mr = mag.min(), mag.max()
if component == 'phase': ydata = normalise(phase(ydata))
elif component == 'magnitude': ydata = normalise(mag)
elif component == 'real': ydata = normalise(ydata.real, *mr)
elif component == 'imaginary': ydata = normalise(ydata.imag, *mr)
elif component == 'real': ydata = ydata.real
elif component == 'imaginary': ydata = ydata.imag
elif component == 'magnitude': ydata = magnitude(ydata)
elif component == 'phase': ydata = phase(ydata)
return xdata, ydata
[docs]class ImaginaryPowerSpectrumSeries(dataseries.DataSeries):
"""An ``ImaginaryPowerSpectrumSeries`` represents the power spectrum of the
imaginary component of a complex-valued image.
``ImaginaryPowerSpectrumSeries`` instances are created by
:class:`ComplexPowerSpectrumSeries` instances.
"""
[docs] def __init__(self, parent, *args, **kwargs):
"""Create an ``ImaginaryPowerSpectrumSeries``.
:arg parent: The :class:`ComplexPowerSpectrumSeries` which owns this
``ImaginaryPowerSpectrumSeries``.
All other arguments are passed through to the :class:`DataSeries`
constructor.
"""
dataseries.DataSeries.__init__(self, *args, **kwargs)
self.__parent = parent
[docs] def makeLabel(self):
"""Returns a string representation of this
``ImaginaryPowerSpectrumSeries`` instance.
"""
return '{} ({})'.format(self.__parent.makeLabelBase(),
strings.labels[self])
[docs] def getData(self):
"""Returns the imaginary component of the power spectrum. """
return self.__parent.getData('imaginary')
[docs]class MagnitudePowerSpectrumSeries(dataseries.DataSeries):
"""An ``MagnitudePowerSpectrumSeries`` represents the magnitude of a
complex-valued image. ``MagnitudePowerSpectrumSeries`` instances are
created by :class:`ComplexPowerSpectrumSeries` instances.
"""
[docs] def __init__(self, parent, *args, **kwargs):
"""Create an ``ImaginaryPowerSpectrumSeries``.
:arg parent: The :class:`ComplexPowerSpectrumSeries` which owns this
``ImaginaryPowerSpectrumSeries``.
All other arguments are passed through to the :class:`DataSeries`
constructor.
"""
dataseries.DataSeries.__init__(self, *args, **kwargs)
self.__parent = parent
[docs] def makeLabel(self):
"""Returns a string representation of this
``MagnitudePowerSpectrumSeries`` instance.
"""
return '{} ({})'.format(self.__parent.makeLabelBase(),
strings.labels[self])
[docs] def getData(self):
"""Returns the magnitude of the complex power spectrum. """
return self.__parent.getData('magnitude')
[docs]class PhasePowerSpectrumSeries(dataseries.DataSeries):
"""An ``PhasePowerSpectrumSeries`` represents the phase of a complex-valued
image. ``PhasePowerSpectrumSeries`` instances are created by
:class:`ComplexPowerSpectrumSeries` instances.
"""
[docs] def __init__(self, parent, *args, **kwargs):
"""Create an ``ImaginaryPowerSpectrumSeries``.
:arg parent: The :class:`ComplexPowerSpectrumSeries` which owns this
``ImaginaryPowerSpectrumSeries``.
All other arguments are passed through to the :class:`DataSeries`
constructor.
"""
dataseries.DataSeries.__init__(self, *args, **kwargs)
self.__parent = parent
[docs] def makeLabel(self):
"""Returns a string representation of this ``PhasePowerSpectrumSeries``
instance.
"""
return '{} ({})'.format(self.__parent.makeLabelBase(),
strings.labels[self])
[docs] def getData(self):
"""Returns the phase of the complex power spectrum. """
return self.__parent.getData('phase')
[docs]class MelodicPowerSpectrumSeries(dataseries.DataSeries,
PowerSpectrumSeries):
"""The ``MelodicPowerSpectrumSeries`` class encapsulates the power spectrum
of the time course for a single component of a :class:`.MelodicImage`. The
component is dictated by the :attr:`.NiftiOpts.volume` property.
"""
[docs] def __init__(self, *args, **kwargs):
"""Create a ``MelodicPowerSpectrumSeries``. All arguments are passed
through to the :meth:`PowerSpectrumSeries.__init__` method.
"""
dataseries.DataSeries.__init__(self, *args, **kwargs)
if not isinstance(self.overlay, fslmelimage.MelodicImage):
raise ValueError('Overlay is not a MelodicImage')
self.varNorm = False
self.disableProperty('varNorm')
[docs] def makeLabel(self):
"""Returns a label that can be used for this
``MelodicPowerSpectrumSeries``.
"""
display = self.displayCtx.getDisplay(self.overlay)
opts = display.opts
component = opts.volume
return '{} [component {}]'.format(display.name, component + 1)
[docs] def getData(self):
"""Returns the power spectrum for the current component of the
:class:`.MelodicImage`, as defined by the :attr:`.NiftiOpts.volume`
property.
"""
overlay = self.overlay
opts = self.displayCtx.getOpts(overlay)
component = opts.volume
ydata = overlay.getComponentPowerSpectrum(component)
nsamples = len(overlay.getComponentTimeSeries(component))
xdata = calcFrequencies(nsamples, overlay.tr, np.float32)
return xdata, ydata
[docs]class MeshPowerSpectrumSeries(dataseries.DataSeries,
PowerSpectrumSeries):
"""A ``MeshPowerSpectrumSeries`` object encapsulates the power spectrum for
the data from a :class:`.Mesh` overlay which has some time series
vertex data associated with it. See the :attr:`.MeshOpts.vertexData`
property.
"""
[docs] def __init__(self, *args, **kwargs):
"""Create a ``MeshPowerSpectrumSeries`` instance. All arguments are
passed through to :meth:`.DataSeries.__init__`.
"""
dataseries.DataSeries.__init__(self, *args, **kwargs)
[docs] def makeLabel(self):
"""Returns a label to use for this ``MeshPowerSpectrumSeries`` on the
legend.
"""
display = self.displayCtx.getDisplay(self.overlay)
if self.__haveData():
opts = display.opts
vidx = opts.getVertex()
return '{} [{}]'.format(display.name, vidx)
else:
return display.name
def __haveData(self):
"""Returns ``True`` if there is currently time series data to show
for this ``MeshPowerSpectrumSeries``, ``False`` otherwise.
"""
opts = self.displayCtx.getOpts(self.overlay)
vidx = opts.getVertex()
vd = opts.getVertexData()
return vidx is not None and vd is not None and vd.shape[1] > 1
[docs] def getData(self):
"""Returns the power spectrum of the data at the current location for
the :class:`.Mesh`, or ``None, None`` if there is no data.
"""
if not self.__haveData():
return None, None
opts = self.displayCtx.getOpts(self.overlay)
vidx = opts.getVertex()
if vidx is None:
return None, None
vd = opts.getVertexData()
data = vd[vidx, :]
xdata = calcFrequencies( len(data), self.sampleTime, data.dtype)
ydata = calcPowerSpectrum(data)
return xdata, ydata