diff --git a/docs/CHANGELOG.md b/docs/CHANGELOG.md index 15d45b66..dc3d32ce 100644 --- a/docs/CHANGELOG.md +++ b/docs/CHANGELOG.md @@ -4,6 +4,7 @@ This is an UNSTABLE release and not recommended for production use (please use rc.111 instead). * jfjoch_broker: Default spot finding settings can be configured via config JSON +* jfjoch_viewer: FFT analysis of data in the dataset plot ### 1.0.0-rc.114 This is an UNSTABLE release and not recommended for production use (please use rc.111 instead). diff --git a/viewer/CMakeLists.txt b/viewer/CMakeLists.txt index 14d8039f..4587d32a 100644 --- a/viewer/CMakeLists.txt +++ b/viewer/CMakeLists.txt @@ -109,3 +109,10 @@ INSTALL( ) qt_import_plugins(jfjoch_viewer INCLUDE Qt::QXcbIntegrationPlugin) + +IF(HAS_FFTW3_H AND FFTWF_LIBRARY) + TARGET_LINK_LIBRARIES(jfjoch_viewer ${FFTWF_LIBRARY}) + MESSAGE(STATUS "FFT single-precision library found: ${FFTWF_LIBRARY}") +ELSE() + MESSAGE(WARNING "FFTW disabled") +ENDIF() \ No newline at end of file diff --git a/viewer/JFJochViewerDatasetInfo.cpp b/viewer/JFJochViewerDatasetInfo.cpp index 37cd69c8..fd36a808 100644 --- a/viewer/JFJochViewerDatasetInfo.cpp +++ b/viewer/JFJochViewerDatasetInfo.cpp @@ -129,7 +129,7 @@ void JFJochViewerDatasetInfo::UpdatePlot() { } else if (val == 8) data = dataset->b_factor; - chart_view->loadValues(data, image_number, one_over_d2, dataset->experiment.GetGoniometer()); + chart_view->loadValues(data, image_number, one_over_d2, dataset.get()); if (dataset->experiment.GetGridScan()) { stack->setCurrentWidget(grid_scan_image); grid_scan_image->loadData(data, dataset->experiment.GetGridScan().value(), one_over_d2); diff --git a/viewer/charts/JFJochDatasetInfoChartView.cpp b/viewer/charts/JFJochDatasetInfoChartView.cpp index 755c9032..48af1a8c 100644 --- a/viewer/charts/JFJochDatasetInfoChartView.cpp +++ b/viewer/charts/JFJochDatasetInfoChartView.cpp @@ -58,6 +58,18 @@ void JFJochDatasetInfoChartView::updateChart() { delete m_hoverLine; m_hoverLine = nullptr; } + +#ifdef JFJOCH_USE_FFTW + if (m_showFFT) { + buildFFTChart(); + return; + } +#endif + + buildTimeDomainChart(); +} + +void JFJochDatasetInfoChartView::buildTimeDomainChart() { if (values.size() >= binning) { // At least one full point @@ -236,6 +248,7 @@ void JFJochDatasetInfoChartView::updateChart() { } } + void JFJochDatasetInfoChartView::setBinning(int64_t val) { if (val >= 1) { binning = val; @@ -260,6 +273,14 @@ void JFJochDatasetInfoChartView::contextMenuEvent(QContextMenuEvent *event) { actXGoniometer->setChecked(m_xUseGoniometerAxis); actXGoniometer->setEnabled(goniometer_axis.has_value()); +#ifdef JFJOCH_USE_FFTW + QAction *actShowFFT = menu.addAction("Show FFT (amplitude vs Hz)"); + actShowFFT->setCheckable(true); + actShowFFT->setChecked(m_showFFT); + // Require valid sampling interval + actShowFFT->setEnabled(!values.empty() && image_time_us > 0.0); +#endif + QAction *chosen = menu.exec(event->globalPos()); if (chosen == copyXY) { QString out; @@ -278,7 +299,13 @@ void JFJochDatasetInfoChartView::contextMenuEvent(QContextMenuEvent *event) { } else if (chosen == actXGoniometer) { m_xUseGoniometerAxis = !m_xUseGoniometerAxis; updateChart(); +#ifdef JFJOCH_USE_FFTW + } else if (chosen == actShowFFT) { + m_showFFT = !m_showFFT; + updateChart(); +#endif } + } void JFJochDatasetInfoChartView::mouseMoveEvent(QMouseEvent *event) { @@ -286,6 +313,62 @@ void JFJochDatasetInfoChartView::mouseMoveEvent(QMouseEvent *event) { if (!series || values.empty()) return; +#ifdef JFJOCH_USE_FFTW + if (m_showFFT && !m_fftFrequenciesHz.empty()) { + // FFT mode: x is frequency in Hz + const QPointF chartPos = chart()->mapToValue(event->pos(), series); + double f = chartPos.x(); + if (!std::isfinite(f)) + return; + + // If we only have DC, nothing meaningful to show + if (m_fftFrequenciesHz.size() <= 1) + return; + + // Find nearest FFT bin, excluding k = 0 (DC component) + int64_t bestIdx = -1; + double bestDiff = std::numeric_limits::infinity(); + for (size_t i = 1; i < m_fftFrequenciesHz.size(); ++i) { + const double diff = std::abs(m_fftFrequenciesHz[i] - f); + if (diff < bestDiff) { + bestDiff = diff; + bestIdx = static_cast(i); + } + } + if (bestIdx < 1) + return; + + const double fBin = m_fftFrequenciesHz[static_cast(bestIdx)]; + const double amp = m_fftMagnitudes[static_cast(bestIdx)]; + + // Map x position of that bin to scene coords for vertical line + const QRectF plotArea = chart()->plotArea(); + const QPointF ptOnChart = chart()->mapToPosition(QPointF(fBin, 0.0), series); + + if (!m_hoverLine) { + m_hoverLine = new QGraphicsLineItem; + m_hoverLine->setPen(QPen(QColor(200, 0, 0, 150), 1.0)); + chart()->scene()->addItem(m_hoverLine); + } + + m_hoverLine->setLine(QLineF(ptOnChart.x(), plotArea.top(), + ptOnChart.x(), plotArea.bottom())); + + QString text = QString("f = %1 Hz, amplitude = %2") + .arg(fBin, 0, 'g', 6) + .arg(amp, 0, 'g', 6); + emit writeStatusBar(text, 6000); + + // No image loading in FFT mode + m_hoverLoadTimer->stop(); + m_hoverPendingIdx = -1; + return; + } +#endif + + if (values.empty()) + return; + // Map mouse position to chart coordinates const QPointF chartPos = chart()->mapToValue(event->pos(), series); double xVal = chartPos.x(); @@ -366,3 +449,93 @@ void JFJochDatasetInfoChartView::onHoverLoadTimeout() { } } } + +#ifdef JFJOCH_USE_FFTW +void JFJochDatasetInfoChartView::buildFFTChart() { + const size_t N = values.size(); + + if (N == 0 || !image_time_us.has_value() || image_time_us <= 0.0) { + return; + } + + // Prepare input buffer (single precision, NaN/inf treated as 0) + std::vector in(N, 0.0f); + for (size_t i = 0; i < N; ++i) { + const double v = values[i]; + in[i] = std::isfinite(v) ? static_cast(v) : 0.0f; + } + + const int n = static_cast(N); + const int nComplex = n / 2 + 1; + + std::vector out(static_cast(nComplex)); + + fftwf_plan plan = fftwf_plan_dft_r2c_1d( + n, + in.data(), + out.data(), + FFTW_ESTIMATE); + + if (!plan) { + return; + } + + fftwf_execute(plan); + fftwf_destroy_plan(plan); + + // Compute amplitude spectrum and frequencies (0 .. Nyquist) + m_fftMagnitudes.resize(static_cast(nComplex)); + m_fftFrequenciesHz.resize(static_cast(nComplex)); + + const double dt = image_time_us.value() * 1e-6; // seconds per sample + const double fs = 1.0 / dt; // sampling frequency + const double df = fs / static_cast(n); // frequency resolution + + for (int k = 0; k < nComplex; ++k) { + const double re = out[static_cast(k)][0]; + const double im = out[static_cast(k)][1]; + const double mag = std::hypot(re, im); // amplitude + + m_fftMagnitudes[static_cast(k)] = mag; + m_fftFrequenciesHz[static_cast(k)] = static_cast(k) * df; + } + + // Build chart series: X = frequency (Hz), Y = amplitude + series = new QLineSeries(this); + currentSeries = nullptr; // no "current image" marker in FFT mode + + double magMin = std::numeric_limits::infinity(); + double magMax = -std::numeric_limits::infinity(); + + for (int k = 1; k < nComplex; ++k) { + const double f = m_fftFrequenciesHz[static_cast(k)]; + const double mag = m_fftMagnitudes[static_cast(k)]; + series->append(f, mag); + if (mag < magMin) magMin = mag; + if (mag > magMax) magMax = mag; + } + + chart()->addSeries(series); + chart()->createDefaultAxes(); + + QValueAxis *axisX = qobject_cast(chart()->axisX(series)); + QValueAxis *axisY = qobject_cast(chart()->axisY(series)); + + if (axisX) { + axisX->setTitleText(QStringLiteral("Frequency (Hz)")); + axisX->setLabelsVisible(true); + } + + if (axisY) { + if (std::isfinite(magMin) && std::isfinite(magMax)) { + if (!(magMax > magMin)) { + magMax = magMin + 1.0; + } + axisY->setRange(magMin, magMax); + } + axisY->setTitleText(QStringLiteral("Amplitude")); + axisY->setLabelsVisible(true); + } +} +#endif + diff --git a/viewer/charts/JFJochDatasetInfoChartView.h b/viewer/charts/JFJochDatasetInfoChartView.h index b5850a82..fad25f07 100644 --- a/viewer/charts/JFJochDatasetInfoChartView.h +++ b/viewer/charts/JFJochDatasetInfoChartView.h @@ -3,6 +3,7 @@ #ifndef JFJOCH_JFJOCHCHARTVIEW_H #define JFJOCH_JFJOCHCHARTVIEW_H + #include #include #include @@ -10,7 +11,13 @@ #include #include #include + +#ifdef JFJOCH_USE_FFTW +#include +#endif + #include "../../common/GoniometerAxis.h" +#include "../../reader/JFJochReaderDataset.h" class JFJochDatasetInfoChartView : public QChartView { Q_OBJECT @@ -30,8 +37,26 @@ class JFJochDatasetInfoChartView : public QChartView { bool m_xUseGoniometerAxis = true; std::optional goniometer_axis; + std::optional image_time_us; + void updateChart(); + // Build regular (index or goniometer) chart + void buildTimeDomainChart(); + + // FFT view toggle + bool m_showFFT = false; + +#ifdef JFJOCH_USE_FFTW + // Cached FFT data for hover/status bar etc. + void buildFFTChart(); + std::vector m_fftMagnitudes; + std::vector m_fftFrequenciesHz; +#endif + + + + signals: void imageSelected(int64_t number); void writeStatusBar(QString string, int timeout_ms = 0); @@ -55,12 +80,21 @@ public slots: public: template - void loadValues(const std::vector &input, int64_t image, bool one_over_d2, const std::optional &goniometer = {}) { + void loadValues(const std::vector &input, + int64_t image, + bool one_over_d2, + const JFJochReaderDataset *dataset = nullptr) { m_yOneOverD = one_over_d2; values.resize(input.size()); - goniometer_axis = goniometer; + if (dataset != nullptr) { + goniometer_axis = dataset->experiment.GetGoniometer(); + image_time_us = dataset->experiment.GetImageTime().count(); + } else { + goniometer_axis = {}; + image_time_us = {}; + } for (int i = 0; i < input.size(); i++) { if (one_over_d2) {