Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Update frequency_correction.py #47

Open
wants to merge 2 commits into
base: master
Choose a base branch
from
Open

Update frequency_correction.py #47

Changes from 1 commit
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
fde37aa
Update frequency_correction.py
laurajanem Mar 31, 2017
d1fbbdf
Update frequency_correction module
laurajanem Apr 3, 2017
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
106 changes: 89 additions & 17 deletions suspect/processing/frequency_correction.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -13,23 +13,43 @@ def residual_water_alignment(data):
return peak_index * data.df


def spectral_registration(data, target, initial_guess=(0.0, 0.0), frequency_range=None):
def spectral_registration(data, target, initial_guess=(0.0, 0.0), frequency_range=None, time_range=None):
"""
Performs the spectral registration method to calculate the frequency and
phase shifts between the input data and the reference spectrum target. The
frequency range over which the two spectra are compared can be specified to
exclude regions where the spectra differ.

frequency range and time range over which the two spectra are compared can
be specified to exclude regions where the spectra differ, or the signal
quality is low.

Algorithm described in : Near, Jamie, et al. "Frequency and phase drift
correction of magnetic resonance spectroscopy data by spectral registration
in the time domain." Magnetic resonance in medicine 73.1 (2015): 44-50.

Parameters
----------
data : MRSData
Input signal to align

target : MRSData
Input signal data will be aligned to

initial_guess : tuple
frequency_range :

Initial guess for the optimal (frequency,phase) shifts, used as a starting
point for the scipy.optimize.leastsq function

frequency_range : tuple
(lower bounds, upper bound) of the frequency range to consider for
alignment. Specified in Hz.

time_range : tuple
(lower bound, upper bound) of the time range to use for alignment.

Returns
-------

freq_and_phase_params : tuple
Optimal (frequency shift, phase shift) parameters output by the
least squares algorithm

"""

# make sure that there are no extra dimensions in the data
Expand All @@ -41,26 +61,78 @@ def spectral_registration(data, target, initial_guess=(0.0, 0.0), frequency_rang
# case we use the spectral points between those two frequencies, or it can
# be a numpy.array of the same size as the data in which case we simply use
# that array as the weightings for the comparison

# Specify some default spectral weights (i.e. equal weighting to all region of the spectrum)
spectral_weights = numpy.ones(data.shape)
if type(frequency_range) is tuple:
spectral_weights = frequency_range[0] < data.frequency_axis() & data.frequency_axis() < frequency_range[1]
else:
spectral_weights = numpy.logical_and(data.frequency_axis() > frequency_range[0],data.frequency_axis() < frequency_range[1])
elif frequency_range is not None:
# If frequency_range is a numpy array...
spectral_weights = frequency_range

# Apply the spectral weights to the data
d = data.spectrum()*spectral_weights
t = target.spectrum()*spectral_weights

# Convert back to the time-domain MRSData object
data = d.fid()
target = t.fid()

# Check on time_axis constraints
# the supplied time range can be None, in which case the whole signal will
# used, or it can be a tuple specifying the lower and upper bounds of the
# time range to use.
if type(time_range) is tuple:
# Run the algorithm on a subset of the time domain samples only, defined by the input range
time_idx = numpy.logical_and(data.time_axis()>=time_range[0],data.time_axis()<=time_range[1])
data = data[time_idx]
target = target[time_idx]


# define a residual function for the optimizer to use
def residual(input_vector):
#transformed_data = transform_fid(data, -input_vector[0], -input_vector[1])
transformed_data = data.adjust_frequency(-input_vector[0]).adjust_phase(-input_vector[1])
residual_data = transformed_data - target
if frequency_range is not None:
spectrum = residual_data.spectrum()
weighted_spectrum = residual_data * spectral_weights
# remove zero-elements
weighted_spectrum = weighted_spectrum[weighted_spectrum != 0]
residual_data = numpy.fft.ifft(numpy.fft.ifftshift(weighted_spectrum))
return_vector = numpy.zeros(len(residual_data) * 2)
return_vector[:len(residual_data)] = residual_data.real
return_vector[len(residual_data):] = residual_data.imag
return return_vector

out = scipy.optimize.leastsq(residual, initial_guess)
return out[0][0], out[0][1]
freq_and_phase_params = (out[0][0], out[0][1])
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

this seems like an unnecessary extra line, what is the purpose of adding a variable name here and not just returning the parameters directly?

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

It is definitely unnecessary - the only reason I added it is because it helps describe the meaning of the outputs a bit, and also makes the documentation a bit easier to read and understand. i.e. in the Returns section, the explanation of the outputs is simplified by capturing them in a single variable, if that makes sense....

Copy link
Member

@bennyrowland bennyrowland Apr 1, 2017 •
edited
Loading

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

One thing to remember about Python is that it will automatically pack and unpack these tuples, so the existing line return out[0][0], out[0][1] is equivalent to return (out[0][0], out[0][1]) and similarly calling it like frequency, phase = spectral_registration() is the same as frequency_and_phase = spectral_registration(). If I wanted to make it really easy to read then I would probably suggest

frequency_shift = out[0][0]
phase_shift = out[0][1]
return frequency_shift, phase_shift

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I like that last suggestion, I'll update it to that. Basically, I usually try to write code that you can read - at the expense of it being less than efficient sometimes. I just felt that when I got down to the output statement, I wasn't sure which of those outputs corresponded to which parameter, and I had to go back up into the code to the place where they are applied to the data with the exponential to figure it out. So, I was just hoping to save someone else from having to do that too : )


return freq_and_phase_params

def select_target_for_spectal_registration(data, frequency_range = None):
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

you are missing an r in spectral

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Looking at this function in detail, I am not sure whether it belongs here, or perhaps needs to be renamed. This is one way of selecting a target for spectral registration, but not the only reasonable choice. I would probably expect this function to be defined at the script level, rather than the library level, to reflect the fact that it is a particular user's preference. On the other hand, it could be renamed to more precisely describe its function, e.g. find_median_magnitude_spectrum(), which suggests it could later be accompanied by other similar functions using different criteria, whereas at the moment it comes across as a bit too prescriptive, the only possible way of selecting a target. Also it might be more appropriate for this to be a function taking a set of spectra rather than FIDs given that it is selecting based on a peak?

This is just me offering my opinion on this, very happy to have a discussion about the best way to structure things. Let me know what you think.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I'm fine with changing the name to be more specific, or, perhaps adding an input to the function as is specifying which method to use (in this case the median method is the only one implemented), so more methods can be added over time without having to add entirely new functions. In that case, I'd probably still want to pass in the MRSData object, not the spectrum, in case future methods want to make use of the FID for some reason. What do you think?

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I would make an analogy with the spectral_registration() method - we chose not to call it do_frequency_correction() because there are other types of frequency correction which people might choose to use. By providing multiple different methods we suit different situations. Personally I think that having a set of methods is much better than the method as parameter option (do_frequency_correction(method="spectral_registration")) - it keeps both the code and the documentation simpler and makes it easier to pass different arguments to different functions, for example passing the spectrum in some cases and the FID in others.

"""
Identifies a candidate target signal from an ensemble that can be used for
spectral registration. The candidate is the signal whose maximum magnitude in the
frequency domain is closest to the median maximum abs values of all the
signals in the ensemble.

Parameters
----------
data : MRSData
Ensemble of candidate signals from which the target is selected.

frequency_range : tuple
Specifies the upper and lower bounds in the frequency domain to use for
selecting the target. Range is specified in Hz.

Returns
-------
target_idx : integer
Index of the selected target signal

"""
if type(frequency_range) is tuple:
spectral_weights = numpy.logical_and(data.frequency_axis(),frequency_range[0],data.frequency_axis() < frequency_range[1])
else:
spectral_weights = frequency_range

filtered_spectrum = spectral_weights*data.spectrum()
max_vals = numpy.max(numpy.abs(filtered_spectrum),axis=1)
target_idx = numpy.argmin(numpy.abs(numpy.median(max_vals)-max_vals))

return target_idx