15#include "MantidHistogramData/Histogram.h"
23#include <gsl/gsl_errno.h>
36using namespace Kernel;
38using namespace DataObjects;
39using namespace HistogramData;
44 "The name of the input workspace.");
46 "The name of the output workspace.");
50 "The name of the input workspace for the imaginary part. "
51 "Leave blank if same as InputWorkspace");
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");
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");
62 "Apply an extra phase equal to this quantity "
63 "times 2*pi to the transform");
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",
89 const bool isComplex = iImag !=
EMPTY_INT();
91 const auto &xPoints =
m_inWS->points(iReal);
92 const auto nPoints =
static_cast<int>(xPoints.size());
94 std::vector<double> data(2 * nPoints);
95 const std::string transform =
getProperty(
"Transform");
99 bool addPositiveOnly =
false;
101 if (!isComplex && transform ==
"Forward") {
103 addPositiveOnly =
true;
106 m_outWS = create<HistoWorkspace>(*
m_inWS, nOut, Points(nPoints));
108 for (
int i = 0; i < nOut; ++i)
109 m_outWS->getSpectrum(i).setDetectorID(
static_cast<detid_t>(i + 1));
111 const double dx = xPoints[1] - xPoints[0];
112 double df = 1.0 / (dx * nPoints);
119 const int dys = nPoints % 2;
121 m_wavetable = gsl_fft_complex_wavetable_alloc(nPoints);
122 m_workspace = gsl_fft_complex_workspace_alloc(nPoints);
127 const bool centerShift =
true;
129 if (transform ==
"Forward") {
130 transformForward(data, nPoints, nPoints, dys, addPositiveOnly, centerShift, isComplex, iReal, iImag, df, dx);
132 transformBackward(data, nPoints, nPoints, dys, centerShift, isComplex, iReal, iImag, df);
138 if (addPositiveOnly) {
150 const bool addPositiveOnly,
const bool centerShift,
const bool isComplex,
const int iReal,
151 const int iImag,
const double df,
const double dx) {
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];
168 data[2 * i + 1] = isComplex ?
m_inImagWS->y(iImag)[j] : 0.;
185 for (
int i = 0; i < ySize; i++) {
186 int j = (ySize / 2 + i + dys) % ySize;
188 double re = data[2 * j] * dx;
189 double im = data[2 * j + 1] * dx;
194 double re1 = re * c - im * s;
195 double im1 = re * s + im * c;
202 if (addPositiveOnly) {
207 m_outWS->mutableY(2)[j] = sqrt(re * re + im * im);
215 if (xSize == ySize + 1) {
223 const bool centerShift,
const bool isComplex,
const int iReal,
const int iImag,
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.;
233 for (
int i = 0; i < ySize; i++) {
236 x -= df * (ySize / 2);
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;
244 m_outWS->mutableY(2)[i] = sqrt(re * re + im * im);
246 if (xSize == ySize + 1)
251 auto tAxis = std::make_unique<API::TextAxis>(nOut);
256 if (addPositiveOnly) {
260 tAxis->setLabel(0,
"Real Positive");
261 tAxis->setLabel(1,
"Imag Positive");
262 tAxis->setLabel(2,
"Modulus Positive");
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));
273 auto inputUnit =
m_inWS->getAxis(0)->unit();
275 std::shared_ptr<Kernel::Units::Label> lblUnit =
279 if ((inputUnit->caption() ==
"Energy" || inputUnit->caption() ==
"Energy transfer") &&
280 inputUnit->label() ==
"meV") {
281 lblUnit->setLabel(
"Time",
"ns");
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") {
291 }
else if (inputUnit->caption() ==
"d-Spacing" && inputUnit->label() ==
"Angstrom") {
293 }
else if (inputUnit->caption() ==
"q" && inputUnit->label() ==
"Angstrom^-1") {
296 m_outWS->getAxis(0)->unit() = lblUnit;
310 std::map<std::string, std::string> errors;
316 auto &
X = inWS->x(iReal);
320 errors[
"InputWorkspace"] =
"Input workspace must have at least two values";
326 errors[
"InputWorkspace"] =
"X axis must be linear (all bins have same width)";
331 auto nHist =
static_cast<int>(inWS->getNumberHistograms());
332 if (iReal >= nHist) {
333 errors[
"Real"] =
"Real out of range";
338 if (inWS->blocksize() != inImagWS->blocksize()) {
339 errors[
"Imaginary"] =
"Real and Imaginary sizes do not match";
341 nHist =
static_cast<int>(inImagWS->getNumberHistograms());
343 if (iImag >= nHist) {
344 errors[
"Imaginary"] =
"Imaginary out of range";
361 const bool acceptXRoundingErrors =
getProperty(
"AcceptXRoundingErrors");
362 const double tolerance = acceptXRoundingErrors ? 0.5 : 1e-7;
363 const double warnValue = acceptXRoundingErrors ? 0.1 : -1;
369 if (!acceptXRoundingErrors) {
373 binChecker.
setErrorType(EqualBinsChecker::ErrorType::Individual);
376 const std::string binError = binChecker.
validate();
377 return !binError.empty();
391 const size_t mid = xPoints.size() / 2;
392 shift = -xPoints[mid];
#define DECLARE_ALGORITHM(classname)
void declareProperty(std::unique_ptr< Kernel::Property > p, const std::string &doc="") override
Add a property to the list of managed properties.
TypedValue getProperty(const std::string &name) const override
Get the value of a property.
A property class for workspaces.
double getPhaseShift(const HistogramData::Points &xPoints)
Get phase shift - user supplied or auto-calculated.
std::map< std::string, std::string > validateInputs() override
Perform validation of inputs.
Mantid::API::MatrixWorkspace_sptr m_outWS
void exec() override
Executes the algorithm.
bool areBinWidthsUneven(const HistogramData::BinEdges &xBins) const
Check whether supplied values are evenly spaced.
gsl_fft_complex_wavetable * m_wavetable
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)
gsl_fft_complex_workspace * m_workspace
Mantid::API::MatrixWorkspace_const_sptr m_inWS
void init() override
Initialisation method. Declares properties to be used in algorithm.
void createUnitsLabels(double &df)
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)
Mantid::API::MatrixWorkspace_const_sptr m_inImagWS
void setupTAxis(const int nOut, const bool addPositiveOnly)
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.
int32_t detid_t
Typedef for a detector ID.
@ Input
An input workspace.
@ Output
An output workspace.