Mantid
Loading...
Searching...
No Matches
FFTSmooth.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"
16#include "MantidHistogramData/HistogramBuilder.h"
20
21namespace Mantid::Algorithms {
22
23// Register the class into the algorithm factory
24DECLARE_ALGORITHM(FFTSmooth)
25
26using namespace Kernel;
27using namespace API;
28using namespace DataObjects;
29using namespace HistogramData;
30
33 declareProperty(std::make_unique<WorkspaceProperty<API::MatrixWorkspace>>("InputWorkspace", "", Direction::Input),
34 "The name of the input workspace.");
35 declareProperty(std::make_unique<WorkspaceProperty<API::MatrixWorkspace>>("OutputWorkspace", "", Direction::Output),
36 "The name of the output workspace.");
37
38 auto mustBePositive = std::make_shared<BoundedValidator<int>>();
39 mustBePositive->setLower(0);
40 declareProperty("WorkspaceIndex", 0, mustBePositive, "Workspace index for smoothing");
41
42 std::vector<std::string> type{"Zeroing"};
43 declareProperty("Filter", "Zeroing", std::make_shared<StringListValidator>(type), "The type of the applied filter");
44 declareProperty("Params", "", "The filter parameters");
45}
46
50 m_inWS = getProperty("InputWorkspace");
51 int spec = getProperty("WorkspaceIndex");
52
53 // Save the starting x value so it can be restored after all transforms.
54 double x0 = m_inWS->x(spec)[0];
55
56 // Symmetrize the input spectrum
57 auto dn = static_cast<int>(m_inWS->y(0).size());
58
59 HistogramBuilder builder;
60 builder.setX(m_inWS->x(0).size() + dn);
61 builder.setY(m_inWS->y(0).size() + dn);
62 builder.setDistribution(m_inWS->isDistribution());
63 API::MatrixWorkspace_sptr symmWS = create<Workspace2D>(*m_inWS, 1, builder.build());
64
65 double dx = (m_inWS->x(spec).back() - m_inWS->x(spec).front()) / static_cast<double>(m_inWS->x(spec).size() - 1);
66
67 auto &symX = symmWS->mutableX(0);
68 auto &symY = symmWS->mutableY(0);
69
70 for (int i = 0; i < dn; i++) {
71 symX[dn + i] = m_inWS->x(spec)[i];
72 symY[dn + i] = m_inWS->y(spec)[i];
73
74 symX[dn - i] = x0 - dx * i;
75 symY[dn - i] = m_inWS->y(spec)[i];
76 }
77 symmWS->mutableY(0).front() = m_inWS->y(spec).back();
78 symmWS->mutableX(0).front() = x0 - dx * dn;
79 if (m_inWS->isHistogramData())
80 symmWS->mutableX(0).back() = m_inWS->x(spec).back();
81
82 // Forward Fourier transform
83 auto fft = createChildAlgorithm("RealFFT", 0, 0.5);
84 fft->setProperty("InputWorkspace", symmWS);
85 fft->setProperty("WorkspaceIndex", 0);
86 try {
87 fft->execute();
88 } catch (...) {
89 g_log.error("Error in direct FFT algorithm");
90 throw;
91 }
92
93 m_unfilteredWS = fft->getProperty("OutputWorkspace");
94
95 // Apply the filter
96 std::string type = getProperty("Filter");
97
98 if (type == "Zeroing") {
99 std::string sn = getProperty("Params");
100 int n;
101 if (sn.empty())
102 n = 2;
103 else
104 n = std::stoi(sn);
105 if (n < 1)
106 throw std::invalid_argument("Truncation parameter must be an integer > 1");
107 zero(n);
108 }
109
110 // Backward transform
111 fft = createChildAlgorithm("RealFFT", 0.5, 1.);
112 fft->setProperty("InputWorkspace", m_filteredWS);
113 fft->setProperty("Transform", "Backward");
114 try {
115 fft->execute();
116 } catch (...) {
117 g_log.error("Error in inverse FFT algorithm");
118 throw;
119 }
120 API::MatrixWorkspace_sptr tmpWS = fft->getProperty("OutputWorkspace");
121
122 // Create output
123 builder.setX(m_inWS->x(0).size());
124 builder.setY(m_inWS->y(0).size());
125 builder.setDistribution(m_inWS->isDistribution());
126 API::MatrixWorkspace_sptr outWS = create<MatrixWorkspace>(*m_inWS, 1, builder.build());
127
128 dn = static_cast<int>(tmpWS->blocksize()) / 2;
129
130 outWS->setSharedX(0, m_inWS->sharedX(0));
131 outWS->mutableY(0).assign(tmpWS->y(0).cbegin() + dn, tmpWS->y(0).cend());
132
133 setProperty("OutputWorkspace", outWS);
134}
135
140 auto my = static_cast<int>(m_unfilteredWS->y(0).size());
141 int ny = my / n;
142
143 double f = double(ny) / my;
144
145 if (ny == 0)
146 ny = 1;
147 int nx = m_unfilteredWS->isHistogramData() ? ny + 1 : ny;
148 HistogramBuilder builder;
149 builder.setX(nx);
150 builder.setY(ny);
151 builder.setDistribution(m_unfilteredWS->isDistribution());
152 m_filteredWS = create<MatrixWorkspace>(*m_unfilteredWS, 2, builder.build());
153
154 auto &Yr = m_unfilteredWS->y(0);
155 auto &Yi = m_unfilteredWS->y(1);
156 auto &X = m_unfilteredWS->x(0);
157
158 auto &yr = m_filteredWS->mutableY(0);
159 auto &yi = m_filteredWS->mutableY(1);
160 auto &xr = m_filteredWS->mutableX(0);
161 auto &xi = m_filteredWS->mutableX(1);
162
163 yr.assign(Yr.begin(), Yr.begin() + ny);
164 yi.assign(Yi.begin(), Yi.begin() + ny);
165 xr.assign(X.begin(), X.begin() + nx);
166 xi.assign(X.begin(), X.begin() + nx);
167
168 using std::placeholders::_1;
169 std::transform(yr.begin(), yr.end(), yr.begin(), std::bind(std::multiplies<double>(), _1, f));
170 std::transform(yi.begin(), yi.end(), yi.begin(), std::bind(std::multiplies<double>(), _1, f));
171}
172
177 auto mx = static_cast<int>(m_unfilteredWS->x(0).size());
178 auto my = static_cast<int>(m_unfilteredWS->y(0).size());
179 int ny = my / n;
180
181 if (ny == 0)
182 ny = 1;
183
184 HistogramBuilder builder;
185 builder.setX(mx);
186 builder.setY(my);
187 builder.setDistribution(m_unfilteredWS->isDistribution());
188 m_filteredWS = create<MatrixWorkspace>(*m_unfilteredWS, 2, builder.build());
189
190 m_filteredWS->setSharedX(0, m_unfilteredWS->sharedX(0));
191 m_filteredWS->setSharedX(1, m_unfilteredWS->sharedX(0));
192
193 std::copy(m_unfilteredWS->y(0).cbegin(), m_unfilteredWS->y(0).begin() + ny, m_filteredWS->mutableY(0).begin());
194
195 std::copy(m_unfilteredWS->y(1).cbegin(), m_unfilteredWS->y(1).begin() + ny, m_filteredWS->mutableY(1).begin());
196}
197
198} // namespace Mantid::Algorithms
#define DECLARE_ALGORITHM(classname)
Definition: Algorithm.h:576
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
virtual std::shared_ptr< Algorithm > createChildAlgorithm(const std::string &name, const double startProgress=-1., const double endProgress=-1., const bool enableLogging=true, const int &version=-1)
Create a Child Algorithm.
Definition: Algorithm.cpp:842
Kernel::Logger & g_log
Definition: Algorithm.h:451
A property class for workspaces.
API::MatrixWorkspace_sptr m_inWS
The input workspace.
Definition: FFTSmooth.h:49
void zero(int n)
Smoothing by zeroing.
Definition: FFTSmooth.cpp:176
void exec() override
Executes the algorithm.
Definition: FFTSmooth.cpp:49
API::MatrixWorkspace_sptr m_filteredWS
Temporary workspace for keeping the filtered spectrum.
Definition: FFTSmooth.h:54
void truncate(int n)
Smoothing by truncation.
Definition: FFTSmooth.cpp:139
API::MatrixWorkspace_sptr m_unfilteredWS
Temporary workspace for keeping the unfiltered Fourier transform of the imput spectrum.
Definition: FFTSmooth.h:52
void init() override
Initialisation method. Declares properties to be used in algorithm.
Definition: FFTSmooth.cpp:32
IPropertyManager * setProperty(const std::string &name, const T &value)
Templated method to set the value of a PropertyWithValue.
void error(const std::string &msg)
Logs at error level.
Definition: Logger.cpp:77
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