Mantid
Loading...
Searching...
No Matches
FFT.cpp
Go to the documentation of this file.
1// Mantid Repository : https://github.com/mantidproject/mantid
2//
3// Copyright © 2018 ISIS Rutherford Appleton Laboratory UKRI,
4// NScD Oak Ridge National Laboratory, European Spallation Source,
5// Institut Laue - Langevin & CSNS, Institute of High Energy Physics, CAS
6// SPDX - License - Identifier: GPL - 3.0 +
7//----------------------------------------------------------------------
8// Includes
9//----------------------------------------------------------------------
12#include "MantidAPI/TextAxis.h"
15#include "MantidHistogramData/Histogram.h"
22
23#include <gsl/gsl_errno.h>
24
25#include <algorithm>
26#include <cmath>
27#include <functional>
28#include <numeric>
29#include <sstream>
30
31namespace Mantid::Algorithms {
32
33// Register the class into the algorithm factory
35
36using namespace Kernel;
37using namespace API;
38using namespace DataObjects;
39using namespace HistogramData;
40
42void FFT::init() {
43 declareProperty(std::make_unique<WorkspaceProperty<>>("InputWorkspace", "", Direction::Input),
44 "The name of the input workspace.");
45 declareProperty(std::make_unique<WorkspaceProperty<>>("OutputWorkspace", "", Direction::Output),
46 "The name of the output workspace.");
47 // if desired, provide the imaginary part in a separate workspace.
49 std::make_unique<WorkspaceProperty<>>("InputImagWorkspace", "", Direction::Input, PropertyMode::Optional),
50 "The name of the input workspace for the imaginary part. "
51 "Leave blank if same as InputWorkspace");
52
53 auto mustBePositive = std::make_shared<BoundedValidator<int>>();
54 mustBePositive->setLower(0);
55 declareProperty("Real", 0, mustBePositive, "Spectrum number to use as real part for transform");
56 declareProperty("Imaginary", EMPTY_INT(), mustBePositive, "Spectrum number to use as imaginary part for transform");
57
58 std::vector<std::string> fft_dir{"Forward", "Backward"};
59 declareProperty("Transform", "Forward", std::make_shared<StringListValidator>(fft_dir),
60 "Direction of the transform: forward or backward");
61 declareProperty("Shift", 0.0,
62 "Apply an extra phase equal to this quantity "
63 "times 2*pi to the transform");
64 declareProperty("AutoShift", false,
65 "Automatically calculate and apply phase shift. Zero on the "
66 "X axis is assumed to be in the centre - if it is not, "
67 "setting this property will automatically correct for this.");
68 declareProperty("AcceptXRoundingErrors", false, "Continue to process the data even if X values are not evenly spaced",
70
71 // "Shift" should only be enabled if "AutoShift" is turned off
72 setPropertySettings("Shift", std::make_unique<EnabledWhenProperty>("AutoShift", IS_DEFAULT));
73}
74
80void FFT::exec() {
81 m_inWS = getProperty("InputWorkspace");
82 m_inImagWS = getProperty("InputImagWorkspace");
83
84 if (m_inImagWS == nullptr)
85 m_inImagWS = m_inWS; // workspaces are one and the same
86
87 const int iReal = getProperty("Real");
88 const int iImag = getProperty("Imaginary");
89 const bool isComplex = iImag != EMPTY_INT();
90
91 const auto &xPoints = m_inWS->points(iReal);
92 const auto nPoints = static_cast<int>(xPoints.size());
93
94 std::vector<double> data(2 * nPoints);
95 const std::string transform = getProperty("Transform");
96
97 // The number of spectra in the output workspace
98 int nOut = 3;
99 bool addPositiveOnly = false;
100 // If the input is real add 3 more spectra with positive "frequencies" only
101 if (!isComplex && transform == "Forward") {
102 nOut += 3;
103 addPositiveOnly = true;
104 }
105
106 m_outWS = create<HistoWorkspace>(*m_inWS, nOut, Points(nPoints));
107
108 for (int i = 0; i < nOut; ++i)
109 m_outWS->getSpectrum(i).setDetectorID(static_cast<detid_t>(i + 1));
110
111 const double dx = xPoints[1] - xPoints[0];
112 double df = 1.0 / (dx * nPoints);
113
114 // Output label
116
117 setupTAxis(nOut, addPositiveOnly);
118
119 const int dys = nPoints % 2;
120
121 m_wavetable = gsl_fft_complex_wavetable_alloc(nPoints);
122 m_workspace = gsl_fft_complex_workspace_alloc(nPoints);
123
124 // Hardcoded "centerShift == true" means that the zero on the x axis is
125 // assumed to be in the centre, at point with index i = ySize/2.
126 // Set to false to make zero at i = 0.
127 const bool centerShift = true;
128
129 if (transform == "Forward") {
130 transformForward(data, nPoints, nPoints, dys, addPositiveOnly, centerShift, isComplex, iReal, iImag, df, dx);
131 } else { // Backward
132 transformBackward(data, nPoints, nPoints, dys, centerShift, isComplex, iReal, iImag, df);
133 }
134
135 m_outWS->setSharedX(1, m_outWS->sharedX(0));
136 m_outWS->setSharedX(2, m_outWS->sharedX(0));
137
138 if (addPositiveOnly) {
139 m_outWS->setSharedX(m_iIm, m_outWS->sharedX(m_iRe));
140 m_outWS->setSharedX(m_iAbs, m_outWS->sharedX(m_iRe));
141 }
142
143 gsl_fft_complex_wavetable_free(m_wavetable);
144 gsl_fft_complex_workspace_free(m_workspace);
145
146 setProperty("OutputWorkspace", m_outWS);
147}
148
149void FFT::transformForward(std::vector<double> &data, const int xSize, const int ySize, const int dys,
150 const bool addPositiveOnly, const bool centerShift, const bool isComplex, const int iReal,
151 const int iImag, const double df, const double dx) {
152 /* If we translate the X-axis by -dx*ySize/2 and assume that our function is
153 * periodic
154 * along the X-axis with period equal to ySize, then dataY values must be
155 * rearranged such that
156 * dataY[i] = dataY[(ySize/2 + i + dys) % ySize]. However, we do not
157 * overwrite dataY but
158 * store the rearranged values in array 'data'.
159 * When index 'i' runs from 0 to ySize/2, data[2*i] will store dataY[j] with
160 * j running from
161 * ySize/2 to ySize. When index 'i' runs from ySize/2+1 to ySize, data[2*i]
162 * will store
163 * dataY[j] with j running from 0 to ySize.
164 */
165 for (int i = 0; i < ySize; i++) {
166 int j = centerShift ? (ySize / 2 + i) % ySize : i;
167 data[2 * i] = m_inWS->y(iReal)[j]; // even indexes filled with the real part
168 data[2 * i + 1] = isComplex ? m_inImagWS->y(iImag)[j] : 0.; // odd indexes filled with the imaginary part
169 }
170
171 double shift = getPhaseShift(m_inWS->points(iReal)); // extra phase to be applied to the transform
172
173 gsl_fft_complex_forward(data.data(), 1, ySize, m_wavetable, m_workspace);
174
175 /* The Fourier transform overwrites array 'data'. Recall that the Fourier
176 * transform is
177 * periodic along the frequency axis. Thus, 'data' takes the same values
178 * when index j runs
179 * from ySize/2 to ySize than if index j would run from -ySize/2 to 0. Thus,
180 * for negative
181 * frequencies running from -ySize/2 to 0, we use the values stored in array
182 * 'data'
183 * for index j running from ySize/2 to ySize.
184 */
185 for (int i = 0; i < ySize; i++) {
186 int j = (ySize / 2 + i + dys) % ySize;
187 m_outWS->mutableX(m_iRe)[i] = df * (-ySize / 2 + i); // zero frequency at i = ySize/2
188 double re = data[2 * j] * dx; // use j from ySize/2 to ySize for negative frequencies
189 double im = data[2 * j + 1] * dx;
190 // shift
191 {
192 double c = cos(m_outWS->x(m_iRe)[i] * shift);
193 double s = sin(m_outWS->x(m_iRe)[i] * shift);
194 double re1 = re * c - im * s;
195 double im1 = re * s + im * c;
196 re = re1;
197 im = im1;
198 }
199 m_outWS->mutableY(m_iRe)[i] = re; // real part
200 m_outWS->mutableY(m_iIm)[i] = im; // imaginary part
201 m_outWS->mutableY(m_iAbs)[i] = sqrt(re * re + im * im); // modulus
202 if (addPositiveOnly) {
203 m_outWS->dataX(0)[i] = df * i;
204 if (j < ySize / 2) {
205 m_outWS->mutableY(0)[j] = re; // real part
206 m_outWS->mutableY(1)[j] = im; // imaginary part
207 m_outWS->mutableY(2)[j] = sqrt(re * re + im * im); // modulus
208 } else {
209 m_outWS->mutableY(0)[j] = 0.; // real part
210 m_outWS->mutableY(1)[j] = 0.; // imaginary part
211 m_outWS->mutableY(2)[j] = 0.; // modulus
212 }
213 }
214 }
215 if (xSize == ySize + 1) {
216 m_outWS->mutableX(0)[ySize] = m_outWS->x(0)[ySize - 1] + df;
217 if (addPositiveOnly)
218 m_outWS->mutableX(m_iRe)[ySize] = m_outWS->x(m_iRe)[ySize - 1] + df;
219 }
220}
221
222void FFT::transformBackward(std::vector<double> &data, const int xSize, const int ySize, const int dys,
223 const bool centerShift, const bool isComplex, const int iReal, const int iImag,
224 const double df) {
225 for (int i = 0; i < ySize; i++) {
226 int j = (ySize / 2 + i) % ySize;
227 data[2 * i] = m_inWS->y(iReal)[j];
228 data[2 * i + 1] = isComplex ? m_inImagWS->y(iImag)[j] : 0.;
229 }
230
231 gsl_fft_complex_inverse(data.data(), 1, ySize, m_wavetable, m_workspace);
232
233 for (int i = 0; i < ySize; i++) {
234 double x = df * i;
235 if (centerShift) {
236 x -= df * (ySize / 2);
237 }
238 m_outWS->mutableX(0)[i] = x;
239 int j = centerShift ? (ySize / 2 + i + dys) % ySize : i;
240 double re = data[2 * j] / df;
241 double im = data[2 * j + 1] / df;
242 m_outWS->mutableY(0)[i] = re; // real part
243 m_outWS->mutableY(1)[i] = im; // imaginary part
244 m_outWS->mutableY(2)[i] = sqrt(re * re + im * im); // modulus
245 }
246 if (xSize == ySize + 1)
247 m_outWS->mutableX(0)[ySize] = m_outWS->x(0)[ySize - 1] + df;
248}
249
250void FFT::setupTAxis(const int nOut, const bool addPositiveOnly) {
251 auto tAxis = std::make_unique<API::TextAxis>(nOut);
252
253 m_iRe = 0;
254 m_iIm = 1;
255 m_iAbs = 2;
256 if (addPositiveOnly) {
257 m_iRe = 3;
258 m_iIm = 4;
259 m_iAbs = 5;
260 tAxis->setLabel(0, "Real Positive");
261 tAxis->setLabel(1, "Imag Positive");
262 tAxis->setLabel(2, "Modulus Positive");
263 }
264 tAxis->setLabel(m_iRe, "Real");
265 tAxis->setLabel(m_iIm, "Imag");
266 tAxis->setLabel(m_iAbs, "Modulus");
267 m_outWS->replaceAxis(1, std::move(tAxis));
268}
269
270void FFT::createUnitsLabels(double &df) {
271 m_outWS->getAxis(0)->unit() = UnitFactory::Instance().create("Label");
272
273 auto inputUnit = m_inWS->getAxis(0)->unit();
274 if (inputUnit) {
275 std::shared_ptr<Kernel::Units::Label> lblUnit =
276 std::dynamic_pointer_cast<Kernel::Units::Label>(UnitFactory::Instance().create("Label"));
277 if (lblUnit) {
278
279 if ((inputUnit->caption() == "Energy" || inputUnit->caption() == "Energy transfer") &&
280 inputUnit->label() == "meV") {
281 lblUnit->setLabel("Time", "ns");
282 df /= 2.418e2;
283 } else if (inputUnit->caption() == "Time" && inputUnit->label() == "s") {
284 lblUnit->setLabel("Frequency", "Hz");
285 } else if (inputUnit->caption() == "Frequency" && inputUnit->label() == "Hz") {
286 lblUnit->setLabel("Time", "s");
287 } else if (inputUnit->caption() == "Time" && inputUnit->label() == "microsecond") {
288 lblUnit->setLabel("Frequency", "MHz");
289 } else if (inputUnit->caption() == "Frequency" && inputUnit->label() == "MHz") {
290 lblUnit->setLabel("Time", Units::Symbol::Microsecond);
291 } else if (inputUnit->caption() == "d-Spacing" && inputUnit->label() == "Angstrom") {
292 lblUnit->setLabel("q", Units::Symbol::InverseAngstrom);
293 } else if (inputUnit->caption() == "q" && inputUnit->label() == "Angstrom^-1") {
294 lblUnit->setLabel("d-Spacing", Units::Symbol::Angstrom);
295 }
296 m_outWS->getAxis(0)->unit() = lblUnit;
297 }
298 }
299}
300
309std::map<std::string, std::string> FFT::validateInputs() {
310 std::map<std::string, std::string> errors;
311
312 MatrixWorkspace_const_sptr inWS = getProperty("InputWorkspace");
313 if (inWS) {
314 const int iReal = getProperty("Real");
315 const int iImag = getProperty("Imaginary");
316 auto &X = inWS->x(iReal);
317
318 // check that the workspace isn't empty
319 if (X.size() < 2) {
320 errors["InputWorkspace"] = "Input workspace must have at least two values";
321 } else {
322 // Check that the x values are evenly spaced
323 // If accepting rounding errors, just give a warning if bins are
324 // different.
325 if (areBinWidthsUneven(inWS->binEdges(iReal))) {
326 errors["InputWorkspace"] = "X axis must be linear (all bins have same width)";
327 }
328 }
329
330 // check real, imaginary spectrum numbers and workspace sizes
331 auto nHist = static_cast<int>(inWS->getNumberHistograms());
332 if (iReal >= nHist) {
333 errors["Real"] = "Real out of range";
334 }
335 if (iImag != EMPTY_INT()) {
336 MatrixWorkspace_const_sptr inImagWS = getProperty("InputImagWorkspace");
337 if (inImagWS) {
338 if (inWS->blocksize() != inImagWS->blocksize()) {
339 errors["Imaginary"] = "Real and Imaginary sizes do not match";
340 }
341 nHist = static_cast<int>(inImagWS->getNumberHistograms());
342 }
343 if (iImag >= nHist) {
344 errors["Imaginary"] = "Imaginary out of range";
345 }
346 }
347 }
348
349 return errors;
350}
351
360bool FFT::areBinWidthsUneven(const HistogramData::BinEdges &xBins) const {
361 const bool acceptXRoundingErrors = getProperty("AcceptXRoundingErrors");
362 const double tolerance = acceptXRoundingErrors ? 0.5 : 1e-7;
363 const double warnValue = acceptXRoundingErrors ? 0.1 : -1;
364
365 // TODO should be added to HistogramData as a convenience function
366 Kernel::EqualBinsChecker binChecker(xBins.rawData(), tolerance, warnValue);
367
368 // Compatibility with previous behaviour
369 if (!acceptXRoundingErrors) {
370 // Compare each bin width to the first (not the average)
371 binChecker.setReferenceBin(EqualBinsChecker::ReferenceBin::First);
372 // Use individual errors (not cumulative)
373 binChecker.setErrorType(EqualBinsChecker::ErrorType::Individual);
374 }
375
376 const std::string binError = binChecker.validate();
377 return !binError.empty();
378}
379
387double FFT::getPhaseShift(const HistogramData::Points &xPoints) {
388 double shift = 0.0;
389 const bool autoshift = getProperty("AutoShift");
390 if (autoshift) {
391 const size_t mid = xPoints.size() / 2;
392 shift = -xPoints[mid];
393 } else {
394 shift = getProperty("Shift");
395 }
396 shift *= 2 * M_PI;
397 return shift;
398}
399
400} // namespace Mantid::Algorithms
#define DECLARE_ALGORITHM(classname)
Definition: Algorithm.h:576
double tolerance
void declareProperty(std::unique_ptr< Kernel::Property > p, const std::string &doc="") override
Add a property to the list of managed properties.
Definition: Algorithm.cpp:1913
TypedValue getProperty(const std::string &name) const override
Get the value of a property.
Definition: Algorithm.cpp:2076
A property class for workspaces.
double getPhaseShift(const HistogramData::Points &xPoints)
Get phase shift - user supplied or auto-calculated.
Definition: FFT.cpp:387
std::map< std::string, std::string > validateInputs() override
Perform validation of inputs.
Definition: FFT.cpp:309
Mantid::API::MatrixWorkspace_sptr m_outWS
Definition: FFT.h:78
void exec() override
Executes the algorithm.
Definition: FFT.cpp:80
bool areBinWidthsUneven(const HistogramData::BinEdges &xBins) const
Check whether supplied values are evenly spaced.
Definition: FFT.cpp:360
gsl_fft_complex_wavetable * m_wavetable
Definition: FFT.h:79
void transformForward(std::vector< double > &data, const int xSize, const int ySize, const int dys, const bool addPositiveOnly, const bool centerShift, const bool isComplex, const int iReal, const int iImag, const double df, const double dx)
Definition: FFT.cpp:149
gsl_fft_complex_workspace * m_workspace
Definition: FFT.h:80
Mantid::API::MatrixWorkspace_const_sptr m_inWS
Definition: FFT.h:76
void init() override
Initialisation method. Declares properties to be used in algorithm.
Definition: FFT.cpp:42
void createUnitsLabels(double &df)
Definition: FFT.cpp:270
void transformBackward(std::vector< double > &data, const int xSize, const int ySize, const int dys, const bool centerShift, const bool isComplex, const int iReal, const int iImag, const double df)
Definition: FFT.cpp:222
Mantid::API::MatrixWorkspace_const_sptr m_inImagWS
Definition: FFT.h:77
void setupTAxis(const int nOut, const bool addPositiveOnly)
Definition: FFT.cpp:250
EqualBinsChecker : Checks for evenly spaced bins.
virtual std::string validate() const
Perform validation of the given X array.
virtual void setErrorType(const ErrorType &errorType)
Set whether to use cumulative errors or compare each in turn.
virtual void setReferenceBin(const ReferenceBin &refBinType)
Set whether to compare each bin to the first bin width or the average.
IPropertyManager * setProperty(const std::string &name, const T &value)
Templated method to set the value of a PropertyWithValue.
void setPropertySettings(const std::string &name, std::unique_ptr< IPropertySettings > settings)
static T & Instance()
Return a reference to the Singleton instance, creating it if it does not already exist Creation is do...
static const UnitLabel InverseAngstrom
InverseAngstrom.
static const UnitLabel Microsecond
Microsecond.
static const UnitLabel Angstrom
Angstrom.
std::shared_ptr< const MatrixWorkspace > MatrixWorkspace_const_sptr
shared pointer to the matrix workspace base class (const version)
std::unique_ptr< T > create(const P &parent, const IndexArg &indexArg, const HistArg &histArg)
This is the create() method that all the other create() methods call.
constexpr int EMPTY_INT() noexcept
Returns what we consider an "empty" integer within a property.
Definition: EmptyValues.h:25
int32_t detid_t
Typedef for a detector ID.
Definition: SpectrumInfo.h:21
@ Input
An input workspace.
Definition: Property.h:53
@ Output
An output workspace.
Definition: Property.h:54