13#include <gsl/gsl_errno.h>
14#include <gsl/gsl_fft_halfcomplex.h>
15#include <gsl/gsl_fft_real.h>
17#define REAL(z, i) ((z)[2 * (i)])
18#define IMAG(z, i) ((z)[2 * (i) + 1])
29#include "MantidHistogramData/LinearGenerator.h"
36using namespace Kernel;
42 "The name of the input workspace.");
44 "The name of the output workspace. It contains three "
45 "spectra: the real, the imaginary parts of the transform and "
48 auto mustBePositive = std::make_shared<BoundedValidator<int>>();
49 mustBePositive->setLower(0);
51 "The index of the spectrum in the input workspace to transform.");
53 std::vector<std::string> fft_dir{
"Forward",
"Backward"};
54 declareProperty(
"Transform",
"Forward", std::make_shared<StringListValidator>(fft_dir),
55 R
"(The direction of the transform: "Forward" or "Backward".)");
57 "Ignores the requirement that X bins be linear and of the same size. "
58 "FFT result will not be valid for the X axis, and should be ignored.");
70 int spec = (transform ==
"Forward") ?
getProperty(
"WorkspaceIndex") : 0;
72 const auto &
X = inWS->x(spec);
73 auto ySize =
static_cast<int>(inWS->blocksize());
76 throw std::invalid_argument(
"Property WorkspaceIndex is out of range");
79 double dx = (
X.back() -
X.front()) /
static_cast<double>(
X.size() - 1);
81 for (
size_t i = 0; i <
X.size() - 2; i++)
82 if (std::abs(dx -
X[i + 1] +
X[i]) / dx > 1e-7)
83 throw std::invalid_argument(
"X axis must be linear (all bins have same "
84 "width). This can be ignored if "
85 "IgnoreXBins is set to true.");
90 double df = 1.0 / (dx * ySize);
92 if (transform ==
"Forward") {
93 int yOutSize = ySize / 2 + 1;
94 int xOutSize = inWS->isHistogramData() ? yOutSize + 1 : yOutSize;
95 bool odd = ySize % 2 != 0;
98 auto tAxis = std::make_unique<API::TextAxis>(3);
99 tAxis->setLabel(0,
"Real");
100 tAxis->setLabel(1,
"Imag");
101 tAxis->setLabel(2,
"Modulus");
102 outWS->replaceAxis(1, std::move(tAxis));
104 gsl_fft_real_workspace *
workspace = gsl_fft_real_workspace_alloc(ySize);
105 std::vector<double> data(2 * ySize);
107 auto &yData = inWS->mutableY(spec);
108 for (
int i = 0; i < ySize; i++) {
112 gsl_fft_real_wavetable *
wavetable = gsl_fft_real_wavetable_alloc(ySize);
117 auto &
x = outWS->mutableX(0);
118 auto &y1 = outWS->mutableY(0);
119 auto &y2 = outWS->mutableY(1);
120 auto &y3 = outWS->mutableY(2);
121 for (
int i = 0; i < yOutSize; i++) {
124 double re = i != 0 ? data[j - 1] : data[0];
125 double im = (i != 0 && (odd || i != yOutSize - 1)) ? data[j] : 0;
128 y3[i] = dx * sqrt(re * re + im * im);
130 if (inWS->isHistogramData()) {
131 outWS->mutableX(0)[yOutSize] = outWS->mutableX(0)[yOutSize - 1] + df;
133 outWS->setSharedX(1, outWS->sharedX(0));
134 outWS->setSharedX(2, outWS->sharedX(0));
138 if (inWS->getNumberHistograms() < 2)
139 throw std::runtime_error(
"The input workspace must have at least 2 spectra.");
141 int yOutSize = (ySize - 1) * 2;
142 if (inWS->y(1).back() != 0.0)
144 int xOutSize = inWS->isHistogramData() ? yOutSize + 1 : yOutSize;
145 bool odd = yOutSize % 2 != 0;
147 df = 1.0 / (dx * (yOutSize));
150 auto tAxis = std::make_unique<API::TextAxis>(1);
151 tAxis->setLabel(0,
"Real");
152 outWS->replaceAxis(1, std::move(tAxis));
154 gsl_fft_real_workspace *
workspace = gsl_fft_real_workspace_alloc(yOutSize);
156 auto &xData = outWS->mutableX(0);
157 auto &yData = outWS->mutableY(0);
158 auto &y0 = inWS->mutableY(0);
159 auto &y1 = inWS->mutableY(1);
160 for (
int i = 0; i < ySize; i++) {
164 yData[j - 1] = y0[i];
165 if (odd || i != ySize - 1) {
173 gsl_fft_halfcomplex_wavetable *
wavetable = gsl_fft_halfcomplex_wavetable_alloc(yOutSize);
177 gsl_fft_halfcomplex_wavetable_free(
wavetable);
180 std::generate(xData.begin(), xData.end(), HistogramData::LinearGenerator(0, df));
183 if (outWS->isHistogramData())
184 outWS->mutableX(0)[yOutSize] = outWS->mutableX(0)[yOutSize - 1] + df;
#define DECLARE_ALGORITHM(classname)
gsl_fft_real_wavetable * wavetable
IPeaksWorkspace_sptr workspace
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.
void exec() override
Executes the algorithm.
void init() override
Initialisation method. Declares properties to be used in algorithm.
IPropertyManager * setProperty(const std::string &name, const T &value)
Templated method to set the value of a PropertyWithValue.
static T & Instance()
Return a reference to the Singleton instance, creating it if it does not already exist Creation is do...
std::shared_ptr< MatrixWorkspace > MatrixWorkspace_sptr
shared pointer to the matrix workspace base class
@ Input
An input workspace.
@ Output
An output workspace.