plots/cartesian/XYFourierTransformCurve.cpp

0001 /*
0002     File                 : XYFourierTransformCurve.cpp
0003     Project              : LabPlot
0004     Description          : A xy-curve defined by a Fourier transform
0005     --------------------------------------------------------------------
0006     SPDX-FileCopyrightText: 2016 Stefan Gerlach <stefan.gerlach@uni.kn>
0007     SPDX-FileCopyrightText: 2017 Alexander Semke <alexander.semke@web.de>
0008     SPDX-License-Identifier: GPL-2.0-or-later
0009 */
0010
0011 /*!
0012   \class XYFourierTransformCurve
0013   \brief A xy-curve defined by a Fourier transform
0014
0015   \ingroup worksheet
0016 */
0017
0018 #include "XYFourierTransformCurve.h"
0019 #include "XYFourierTransformCurvePrivate.h"
0020 #include "backend/core/AbstractColumn.h"
0021 #include "backend/core/column/Column.h"
0022 #include "backend/gsl/errors.h"
0023 #include "backend/lib/XmlStreamReader.h"
0024 #include "backend/lib/commandtemplates.h"
0025 #include "backend/lib/macros.h"
0026
0027 extern "C" {
0028 #include "backend/nsl/nsl_sf_poly.h"
0029 }
0030
0031 #include <KLocalizedString>
0032 #include <QDebug> // qWarning()
0033 #include <QElapsedTimer>
0034 #include <QIcon>
0035 #include <QThreadPool>
0036
0037 XYFourierTransformCurve::XYFourierTransformCurve(const QString& name)
0038     : XYAnalysisCurve(name, new XYFourierTransformCurvePrivate(this), AspectType::XYFourierTransformCurve) {
0039 }
0040
0041 XYFourierTransformCurve::XYFourierTransformCurve(const QString& name, XYFourierTransformCurvePrivate* dd)
0042     : XYAnalysisCurve(name, dd, AspectType::XYFourierTransformCurve) {
0043 }
0044
0045 // no need to delete the d-pointer here - it inherits from QGraphicsItem
0046 // and is deleted during the cleanup in QGraphicsScene
0047 XYFourierTransformCurve::~XYFourierTransformCurve() = default;
0048
0049 void XYFourierTransformCurve::recalculate() {
0050     Q_D(XYFourierTransformCurve);
0051     d->recalculate();
0052 }
0053
0054 /*!
0055     Returns an icon to be used in the project explorer.
0056 */
0057 QIcon XYFourierTransformCurve::icon() const {
0058     return QIcon::fromTheme(QStringLiteral("labplot-xy-fourier-transform-curve"));
0059 }
0060
0061 // ##############################################################################
0062 // ##########################  getter methods  ##################################
0063 // ##############################################################################
0064 BASIC_SHARED_D_READER_IMPL(XYFourierTransformCurve, XYFourierTransformCurve::TransformData, transformData, transformData)
0065
0066 const XYAnalysisCurve::Result& XYFourierTransformCurve::result() const {
0067     Q_D(const XYFourierTransformCurve);
0068     return d->transformResult;
0069 }
0070
0071 // ##############################################################################
0072 // #################  setter methods and undo commands ##########################
0073 // ##############################################################################
0074 STD_SETTER_CMD_IMPL_F_S(XYFourierTransformCurve, SetTransformData, XYFourierTransformCurve::TransformData, transformData, recalculate)
0075 void XYFourierTransformCurve::setTransformData(const XYFourierTransformCurve::TransformData& transformData) {
0076     Q_D(XYFourierTransformCurve);
0077     exec(new XYFourierTransformCurveSetTransformDataCmd(d, transformData, ki18n("%1: set transform options and perform the Fourier transform")));
0078 }
0079
0080 // ##############################################################################
0081 // ######################### Private implementation #############################
0082 // ##############################################################################
0083 XYFourierTransformCurvePrivate::XYFourierTransformCurvePrivate(XYFourierTransformCurve* owner)
0084     : XYAnalysisCurvePrivate(owner)
0085     , q(owner) {
0086 }
0087
0088 // no need to delete xColumn and yColumn, they are deleted
0089 // when the parent aspect is removed
0090 XYFourierTransformCurvePrivate::~XYFourierTransformCurvePrivate() = default;
0091
0092 void XYFourierTransformCurvePrivate::resetResults() {
0093     transformResult = XYFourierTransformCurve::TransformResult();
0094 }
0095
0096 bool XYFourierTransformCurvePrivate::recalculateSpecific(const AbstractColumn* tmpXDataColumn, const AbstractColumn* tmpYDataColumn) {
0097     QElapsedTimer timer;
0098     timer.start();
0099
0100     // copy all valid data point for the transform to temporary vectors
0101     QVector<double> xdataVector;
0102     QVector<double> ydataVector;
0103     const double xmin = transformData.xRange.first();
0104     const double xmax = transformData.xRange.last();
0105
0106     int rowCount = std::min(tmpXDataColumn->rowCount(), tmpYDataColumn->rowCount());
0107     for (int row = 0; row < rowCount; ++row) {
0108         // only copy those data where _all_ values (for x and y, if given) are valid
0109         if (std::isnan(tmpXDataColumn->valueAt(row)) || std::isnan(tmpYDataColumn->valueAt(row)) || tmpXDataColumn->isMasked(row)
0110             || tmpYDataColumn->isMasked(row))
0111             continue;
0112
0113         // only when inside given range
0114         if (tmpXDataColumn->valueAt(row) >= xmin && tmpXDataColumn->valueAt(row) <= xmax) {
0115             xdataVector.append(tmpXDataColumn->valueAt(row));
0116             ydataVector.append(tmpYDataColumn->valueAt(row));
0117         }
0118     }
0119
0120     // number of data points to transform
0121     auto n = (unsigned int)ydataVector.size();
0122     if (n == 0) {
0123         transformResult.available = true;
0124         transformResult.valid = false;
0125         transformResult.status = i18n("No data points available.");
0126         return true;
0127     }
0128
0129     double* xdata = xdataVector.data();
0130     double* ydata = ydataVector.data();
0131
0132     // transform settings
0133     const nsl_sf_window_type windowType = transformData.windowType;
0134     const nsl_dft_result_type type = transformData.type;
0135     const bool twoSided = transformData.twoSided;
0136     const bool shifted = transformData.shifted;
0137     const nsl_dft_xscale xScale = transformData.xScale;
0138
0139     DEBUG("n =" << n);
0140     DEBUG("window type:" << nsl_sf_window_type_name[windowType]);
0141     DEBUG("type:" << nsl_dft_result_type_name[type]);
0142     DEBUG("scale:" << nsl_dft_xscale_name[xScale]);
0143     DEBUG("two sided:" << twoSided);
0144     DEBUG("shifted:" << shifted);
0145 #ifndef NDEBUG
0146     QDebug out = qDebug();
0147     for (unsigned int i = 0; i < n; i++)
0148         out << ydata[i];
0149 #endif
0150
0151     ///////////////////////////////////////////////////////////
0152     // transform with window
0153     gsl_set_error_handler_off();
0154     int status = nsl_dft_transform_window(ydata, 1, n, twoSided, type, windowType);
0155
0156     unsigned int N = n;
0157     if (twoSided == false)
0158         N = n / 2;
0159
0160     switch (xScale) {
0161     case nsl_dft_xscale_frequency:
0162         for (unsigned int i = 0; i < N; i++) {
0163             if (i >= n / 2 && shifted)
0164                 xdata[i] = (n - 1) / (xmax - xmin) * (i / (double)n - 1.);
0165             else
0166                 xdata[i] = (n - 1) * i / (xmax - xmin) / n;
0167         }
0168         break;
0169     case nsl_dft_xscale_index:
0170         for (unsigned int i = 0; i < N; i++) {
0171             if (i >= n / 2 && shifted)
0172                 xdata[i] = (int)i - (int)N;
0173             else
0174                 xdata[i] = i;
0175         }
0176         break;
0177     case nsl_dft_xscale_period: {
0178         double f0 = (n - 1) / (xmax - xmin) / n;
0179         for (unsigned int i = 0; i < N; i++) {
0180             double f = (n - 1) * i / (xmax - xmin) / n;
0181             xdata[i] = 1 / (f + f0);
0182         }
0183         break;
0184     }
0185     }
0186 #ifndef NDEBUG
0187     out = qDebug();
0188     for (unsigned int i = 0; i < N; i++)
0189         out << ydata[i] << '(' << xdata[i] << ')';
0190 #endif
0191
0192     xVector->resize((int)N);
0193     yVector->resize((int)N);
0194     if (shifted) {
0195         memcpy(xVector->data(), &xdata[n / 2], n / 2 * sizeof(double));
0196         memcpy(&xVector->data()[n / 2], xdata, n / 2 * sizeof(double));
0197         memcpy(yVector->data(), &ydata[n / 2], n / 2 * sizeof(double));
0198         memcpy(&yVector->data()[n / 2], ydata, n / 2 * sizeof(double));
0199     } else {
0200         memcpy(xVector->data(), xdata, N * sizeof(double));
0201         memcpy(yVector->data(), ydata, N * sizeof(double));
0202     }
0203     ///////////////////////////////////////////////////////////
0204
0205     // write the result
0206     transformResult.available = true;
0207     transformResult.valid = (status == GSL_SUCCESS);
0208     transformResult.status = gslErrorToString(status);
0209     transformResult.elapsedTime = timer.elapsed();
0210
0211     return true;
0212 }
0213
0214 // ##############################################################################
0215 // ##################  Serialization/Deserialization  ###########################
0216 // ##############################################################################
0217 //! Save as XML
0218 void XYFourierTransformCurve::save(QXmlStreamWriter* writer) const {
0219     Q_D(const XYFourierTransformCurve);
0220
0221     writer->writeStartElement(QStringLiteral("xyFourierTransformCurve"));
0222
0223     // write the base class
0224     XYAnalysisCurve::save(writer);
0225
0226     // write xy-fourier_transform-curve specific information
0227     // transform data
0228     writer->writeStartElement(QStringLiteral("transformData"));
0229     writer->writeAttribute(QStringLiteral("autoRange"), QString::number(d->transformData.autoRange));
0230     writer->writeAttribute(QStringLiteral("xRangeMin"), QString::number(d->transformData.xRange.first()));
0231     writer->writeAttribute(QStringLiteral("xRangeMax"), QString::number(d->transformData.xRange.last()));
0232     writer->writeAttribute(QStringLiteral("type"), QString::number(d->transformData.type));
0233     writer->writeAttribute(QStringLiteral("twoSided"), QString::number(d->transformData.twoSided));
0234     writer->writeAttribute(QStringLiteral("shifted"), QString::number(d->transformData.shifted));
0235     writer->writeAttribute(QStringLiteral("xScale"), QString::number(d->transformData.xScale));
0236     writer->writeAttribute(QStringLiteral("windowType"), QString::number(d->transformData.windowType));
0237     writer->writeEndElement(); // transformData
0238
0239     // transform results (generated columns)
0240     writer->writeStartElement(QStringLiteral("transformResult"));
0241     writer->writeAttribute(QStringLiteral("available"), QString::number(d->transformResult.available));
0242     writer->writeAttribute(QStringLiteral("valid"), QString::number(d->transformResult.valid));
0243     writer->writeAttribute(QStringLiteral("status"), d->transformResult.status);
0244     writer->writeAttribute(QStringLiteral("time"), QString::number(d->transformResult.elapsedTime));
0245
0246     // save calculated columns if available
0247     if (saveCalculations() && d->xColumn && d->yColumn) {
0248         d->xColumn->save(writer);
0249         d->yColumn->save(writer);
0250     }
0251     writer->writeEndElement(); //"transformResult"
0252     writer->writeEndElement(); //"xyFourierTransformCurve"
0253 }
0254
0255 //! Load from XML
0256 bool XYFourierTransformCurve::load(XmlStreamReader* reader, bool preview) {
0257     Q_D(XYFourierTransformCurve);
0258
0259     QXmlStreamAttributes attribs;
0260     QString str;
0261
0262     while (!reader->atEnd()) {
0263         reader->readNext();
0264         if (reader->isEndElement() && reader->name() == QLatin1String("xyFourierTransformCurve"))
0265             break;
0266
0267         if (!reader->isStartElement())
0268             continue;
0269
0270         if (reader->name() == QLatin1String("xyAnalysisCurve")) {
0271             if (!XYAnalysisCurve::load(reader, preview))
0272                 return false;
0273         } else if (!preview && reader->name() == QLatin1String("transformData")) {
0274             attribs = reader->attributes();
0275             READ_INT_VALUE("autoRange", transformData.autoRange, bool);
0276             READ_DOUBLE_VALUE("xRangeMin", transformData.xRange.first());
0277             READ_DOUBLE_VALUE("xRangeMax", transformData.xRange.last());
0278             READ_INT_VALUE("type", transformData.type, nsl_dft_result_type);
0279             READ_INT_VALUE("twoSided", transformData.twoSided, bool);
0280             READ_INT_VALUE("shifted", transformData.shifted, bool);
0281             READ_INT_VALUE("xScale", transformData.xScale, nsl_dft_xscale);
0282             READ_INT_VALUE("windowType", transformData.windowType, nsl_sf_window_type);
0283         } else if (!preview && reader->name() == QLatin1String("transformResult")) {
0284             attribs = reader->attributes();
0285             READ_INT_VALUE("available", transformResult.available, int);
0286             READ_INT_VALUE("valid", transformResult.valid, int);
0287             READ_STRING_VALUE("status", transformResult.status);
0288             READ_INT_VALUE("time", transformResult.elapsedTime, int);
0289         } else if (reader->name() == QLatin1String("column")) {
0290             Column* column = new Column(QString(), AbstractColumn::ColumnMode::Double);
0291             if (!column->load(reader, preview)) {
0292                 delete column;
0293                 return false;
0294             }
0295
0296             if (column->name() == QLatin1String("x"))
0297                 d->xColumn = column;
0298             else if (column->name() == QLatin1String("y"))
0299                 d->yColumn = column;
0300         } else { // unknown element
0301             reader->raiseUnknownElementWarning();
0302             if (!reader->skipToEndElement())
0303                 return false;
0304         }
0305     }
0306
0307     if (preview)
0308         return true;
0309
0310     // wait for data to be read before using the pointers
0311     QThreadPool::globalInstance()->waitForDone();
0312
0313     if (d->xColumn && d->yColumn) {
0314         d->xColumn->setHidden(true);
0315         addChild(d->xColumn);
0316
0317         d->yColumn->setHidden(true);
0318         addChild(d->yColumn);
0319
0320         d->xVector = static_cast<QVector<double>*>(d->xColumn->data());
0321         d->yVector = static_cast<QVector<double>*>(d->yColumn->data());
0322
0323         static_cast<XYCurvePrivate*>(d_ptr)->xColumn = d->xColumn;
0324         static_cast<XYCurvePrivate*>(d_ptr)->yColumn = d->yColumn;
0325
0326         recalcLogicalPoints();
0327     }
0328
0329     return true;
0330 }