Mantid
Loading...
Searching...
No Matches
RealFFT.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 +
12
13#include <gsl/gsl_errno.h>
14#include <gsl/gsl_fft_halfcomplex.h>
15#include <gsl/gsl_fft_real.h>
16
17#define REAL(z, i) ((z)[2 * (i)])
18#define IMAG(z, i) ((z)[2 * (i) + 1])
19
20#include <algorithm>
21#include <cmath>
22#include <functional>
23#include <numeric>
24#include <sstream>
25
28
29#include "MantidHistogramData/LinearGenerator.h"
30
31namespace Mantid::Algorithms {
32
33// Register the class into the algorithm factory
34DECLARE_ALGORITHM(RealFFT)
35
36using namespace Kernel;
37using namespace API;
38
41 declareProperty(std::make_unique<WorkspaceProperty<API::MatrixWorkspace>>("InputWorkspace", "", Direction::Input),
42 "The name of the input workspace.");
43 declareProperty(std::make_unique<WorkspaceProperty<API::MatrixWorkspace>>("OutputWorkspace", "", Direction::Output),
44 "The name of the output workspace. It contains three "
45 "spectra: the real, the imaginary parts of the transform and "
46 "their modulus.");
47
48 auto mustBePositive = std::make_shared<BoundedValidator<int>>();
49 mustBePositive->setLower(0);
50 declareProperty("WorkspaceIndex", 0, mustBePositive,
51 "The index of the spectrum in the input workspace to transform.");
52
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".)");
56 declareProperty("IgnoreXBins", false,
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.");
59}
60
66 API::MatrixWorkspace_sptr inWS = getProperty("InputWorkspace");
67 std::string transform = getProperty("Transform");
68 bool IgnoreXBins = getProperty("IgnoreXBins");
69
70 int spec = (transform == "Forward") ? getProperty("WorkspaceIndex") : 0;
71
72 const auto &X = inWS->x(spec);
73 auto ySize = static_cast<int>(inWS->blocksize());
74
75 if (spec >= ySize)
76 throw std::invalid_argument("Property WorkspaceIndex is out of range");
77
78 // Check that the x values are evenly spaced
79 double dx = (X.back() - X.front()) / static_cast<double>(X.size() - 1);
80 if (!IgnoreXBins) {
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.");
86 }
87
89
90 double df = 1.0 / (dx * ySize);
91
92 if (transform == "Forward") {
93 int yOutSize = ySize / 2 + 1;
94 int xOutSize = inWS->isHistogramData() ? yOutSize + 1 : yOutSize;
95 bool odd = ySize % 2 != 0;
96
97 outWS = WorkspaceFactory::Instance().create(inWS, 3, xOutSize, yOutSize);
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));
103
104 gsl_fft_real_workspace *workspace = gsl_fft_real_workspace_alloc(ySize);
105 std::vector<double> data(2 * ySize);
106
107 auto &yData = inWS->mutableY(spec);
108 for (int i = 0; i < ySize; i++) {
109 data[i] = yData[i];
110 }
111
112 gsl_fft_real_wavetable *wavetable = gsl_fft_real_wavetable_alloc(ySize);
113 gsl_fft_real_transform(data.data(), 1, ySize, wavetable, workspace);
114 gsl_fft_real_wavetable_free(wavetable);
115 gsl_fft_real_workspace_free(workspace);
116
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++) {
122 int j = i * 2;
123 x[i] = df * i;
124 double re = i != 0 ? data[j - 1] : data[0];
125 double im = (i != 0 && (odd || i != yOutSize - 1)) ? data[j] : 0;
126 y1[i] = re * dx; // real part
127 y2[i] = im * dx; // imaginary part
128 y3[i] = dx * sqrt(re * re + im * im); // modulus
129 }
130 if (inWS->isHistogramData()) {
131 outWS->mutableX(0)[yOutSize] = outWS->mutableX(0)[yOutSize - 1] + df;
132 }
133 outWS->setSharedX(1, outWS->sharedX(0));
134 outWS->setSharedX(2, outWS->sharedX(0));
135 } else // Backward
136 {
137
138 if (inWS->getNumberHistograms() < 2)
139 throw std::runtime_error("The input workspace must have at least 2 spectra.");
140
141 int yOutSize = (ySize - 1) * 2;
142 if (inWS->y(1).back() != 0.0)
143 yOutSize++;
144 int xOutSize = inWS->isHistogramData() ? yOutSize + 1 : yOutSize;
145 bool odd = yOutSize % 2 != 0;
146
147 df = 1.0 / (dx * (yOutSize));
148
149 outWS = WorkspaceFactory::Instance().create(inWS, 1, xOutSize, yOutSize);
150 auto tAxis = std::make_unique<API::TextAxis>(1);
151 tAxis->setLabel(0, "Real");
152 outWS->replaceAxis(1, std::move(tAxis));
153
154 gsl_fft_real_workspace *workspace = gsl_fft_real_workspace_alloc(yOutSize);
155
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++) {
161 int j = i * 2;
162 xData[i] = df * i;
163 if (i != 0) {
164 yData[j - 1] = y0[i];
165 if (odd || i != ySize - 1) {
166 yData[j] = y1[i];
167 }
168 } else {
169 yData[0] = y0[0];
170 }
171 }
172
173 gsl_fft_halfcomplex_wavetable *wavetable = gsl_fft_halfcomplex_wavetable_alloc(yOutSize);
174
175 // &(yData[0]) because gsl func wants non const double data[]
176 gsl_fft_halfcomplex_inverse(&(yData[0]), 1, yOutSize, wavetable, workspace);
177 gsl_fft_halfcomplex_wavetable_free(wavetable);
178 gsl_fft_real_workspace_free(workspace);
179
180 std::generate(xData.begin(), xData.end(), HistogramData::LinearGenerator(0, df));
181 yData /= df;
182
183 if (outWS->isHistogramData())
184 outWS->mutableX(0)[yOutSize] = outWS->mutableX(0)[yOutSize - 1] + df;
185 }
186
187 setProperty("OutputWorkspace", outWS);
188}
189
190} // namespace Mantid::Algorithms
#define DECLARE_ALGORITHM(classname)
Definition: Algorithm.h:576
gsl_fft_real_wavetable * wavetable
IPeaksWorkspace_sptr workspace
Definition: IndexPeaks.cpp:114
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.
void exec() override
Executes the algorithm.
Definition: RealFFT.cpp:65
void init() override
Initialisation method. Declares properties to be used in algorithm.
Definition: RealFFT.cpp:40
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.
Definition: Property.h:53
@ Output
An output workspace.
Definition: Property.h:54