From fde37aa1cfdbd9921839929dcab1efd4e045bfa1 Mon Sep 17 00:00:00 2001 From: laurajanem Date: Fri, 31 Mar 2017 16:37:31 -0400 Subject: [PATCH 1/2] Update frequency_correction.py Many changes made Added module for automatic identification of a target to use for spectral registration Fixed bug in the way the subset of the frequency range was identified Added additional input to spectral registration module to enable limiting the time range used for alignment Quality of spectral registration is significantly improved over previous version --- suspect/processing/frequency_correction.py | 106 +++++++++++++++++---- 1 file changed, 89 insertions(+), 17 deletions(-) diff --git a/suspect/processing/frequency_correction.py b/suspect/processing/frequency_correction.py index 7830507..56d5b46 100644 --- a/suspect/processing/frequency_correction.py +++ b/suspect/processing/frequency_correction.py @@ -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 @@ -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]) + + return freq_and_phase_params + +def select_target_for_spectal_registration(data, frequency_range = None): + """ + 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 + From d1fbbdf03ae9381e15e0abb8e8c3bfc92a7a4da8 Mon Sep 17 00:00:00 2001 From: laurajanem Date: Mon, 3 Apr 2017 10:41:03 -0400 Subject: [PATCH 2/2] Update frequency_correction module Changed output structure for spectral registration function. Changed name of function for identifying target for spectral registration to be more explicit. Changed input to target identification function to MRSSpectrum object instead of MRSData. --- suspect/processing/frequency_correction.py | 28 +++++++++++++--------- 1 file changed, 17 insertions(+), 11 deletions(-) diff --git a/suspect/processing/frequency_correction.py b/suspect/processing/frequency_correction.py index 56d5b46..3b5fb67 100644 --- a/suspect/processing/frequency_correction.py +++ b/suspect/processing/frequency_correction.py @@ -46,9 +46,11 @@ def spectral_registration(data, target, initial_guess=(0.0, 0.0), frequency_rang Returns ------- - freq_and_phase_params : tuple - Optimal (frequency shift, phase shift) parameters output by the - least squares algorithm + frequency_shift : float + Optimal frequency shift output by the least squares algorithm, in Hz. + + phase_shift : float + Optimal phase shift output by the least squares algorithm, in radians. """ @@ -99,11 +101,12 @@ def residual(input_vector): return return_vector out = scipy.optimize.leastsq(residual, initial_guess) - freq_and_phase_params = (out[0][0], out[0][1]) + frequency_shift = out[0][0] + phase_shift = out[0][1] - return freq_and_phase_params + return frequency_shift, phase_shift -def select_target_for_spectal_registration(data, frequency_range = None): +def select_target_for_spectral_reg_median_mag(spectral_data, frequency_range = None): """ 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 @@ -112,8 +115,9 @@ def select_target_for_spectal_registration(data, frequency_range = None): Parameters ---------- - data : MRSData - Ensemble of candidate signals from which the target is selected. + spectral_data : MRSSpectrum + Ensemble of candidate signals from which the target is selected, in + frequency domain frequency_range : tuple Specifies the upper and lower bounds in the frequency domain to use for @@ -125,14 +129,16 @@ def select_target_for_spectal_registration(data, frequency_range = None): Index of the selected target signal """ + spectral_weights = numpy.ones(spectral_data.shape) 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 = numpy.logical_and(spectral_data.frequency_axis(),frequency_range[0],spectral_data.frequency_axis() < frequency_range[1]) + elif frequency_range is not None: spectral_weights = frequency_range - filtered_spectrum = spectral_weights*data.spectrum() + filtered_spectrum = spectral_weights*spectral_data 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 +