diff --git a/doc/img/NoiseFigure_plugin_enr.png b/doc/img/NoiseFigure_plugin_enr.png index 94b75e2ed..ee2c83dff 100644 Binary files a/doc/img/NoiseFigure_plugin_enr.png and b/doc/img/NoiseFigure_plugin_enr.png differ diff --git a/doc/img/NoiseFigure_plugin_results.png b/doc/img/NoiseFigure_plugin_results.png index 0c300ae6c..7d6778d55 100644 Binary files a/doc/img/NoiseFigure_plugin_results.png and b/doc/img/NoiseFigure_plugin_results.png differ diff --git a/plugins/channelrx/noisefigure/noisefigure.cpp b/plugins/channelrx/noisefigure/noisefigure.cpp index 138b97527..5fc75cb16 100644 --- a/plugins/channelrx/noisefigure/noisefigure.cpp +++ b/plugins/channelrx/noisefigure/noisefigure.cpp @@ -37,6 +37,7 @@ #include "device/deviceapi.h" #include "feature/feature.h" #include "util/db.h" +#include "util/interpolation.h" #include "channel/channelwebapiutils.h" #include "maincore.h" @@ -161,10 +162,10 @@ void NoiseFigure::processVISA(QStringList commands) // Calculate ENR at specified frequency double NoiseFigure::calcENR(double frequency) { + double enr; int size = m_settings.m_enr.size(); if (size >= 2) { - int order = size > 3 ? 3 : size - 1; std::vector x(size); std::vector y(size); for (int i = 0; i < size; i++) @@ -172,19 +173,28 @@ double NoiseFigure::calcENR(double frequency) x[i] = m_settings.m_enr[i]->m_frequency; y[i] = m_settings.m_enr[i]->m_enr; } - boost::math::barycentric_rational interpolant(std::move(x), std::move(y), order); - double enr = interpolant(frequency); - qDebug() << "ENR at " << frequency << " interpolated to " << enr; - return enr; + if (m_settings.m_interpolation == NoiseFigureSettings::LINEAR) + { + enr = Interpolation::linear(x.begin(), x.end(), y.begin(), frequency); + } + else + { + int order = size - 1; + boost::math::barycentric_rational interpolant(std::move(x), std::move(y), order); + enr = interpolant(frequency); + } } else if (size == 1) { - return m_settings.m_enr[0]->m_enr; + enr = m_settings.m_enr[0]->m_enr; } else { - return 0.0; + // Shouldn't get here + enr = 0.0; } + qDebug() << "ENR at " << frequency << " interpolated to " << enr; + return enr; } // FSM for running measurements over multiple frequencies @@ -295,15 +305,19 @@ void NoiseFigure::nextState() case COMPLETE: // Calculate noise figure and temperature using Y-factor method + double t = 290.0; + double k = 1.38064852e-23; + double bw = 1; double y = m_onPower - m_offPower; double enr = calcENR(m_measurementFrequency/1e6); double nf = 10.0*log10(pow(10.0, enr/10.0)/(pow(10.0, y/10.0)-1.0)); - double temp = 290.0*(pow(10.0, nf/10.0)-1.0); + double temp = t*(pow(10.0, nf/10.0)-1.0); + double floor = 10.0*log10(1000.0*k*t) + nf + 10*log10(bw); // Send result to GUI if (getMessageQueueToGUI()) { - MsgNFMeasurement *msg = MsgNFMeasurement::create(m_measurementFrequency/1e6, nf, temp, y, enr); + MsgNFMeasurement *msg = MsgNFMeasurement::create(m_measurementFrequency/1e6, nf, temp, y, enr, floor); getMessageQueueToGUI()->push(msg); } @@ -407,7 +421,7 @@ bool NoiseFigure::handleMessage(const Message& cmd) { if (m_state == IDLE) { - if (!m_settings.m_visaDevice.isEmpty()) + if (!m_settings.m_visaDevice.isEmpty()) { if (openVISADevice()) { QTimer::singleShot(0, this, SLOT(nextState())); @@ -415,7 +429,7 @@ bool NoiseFigure::handleMessage(const Message& cmd) getMessageQueueToGUI()->push(MsgFinished::create(QString("Failed to open VISA device %1").arg(m_settings.m_visaDevice))); } } - else + else { QTimer::singleShot(0, this, SLOT(nextState())); } diff --git a/plugins/channelrx/noisefigure/noisefigure.h b/plugins/channelrx/noisefigure/noisefigure.h index b174bf526..9a1deba49 100644 --- a/plugins/channelrx/noisefigure/noisefigure.h +++ b/plugins/channelrx/noisefigure/noisefigure.h @@ -95,10 +95,11 @@ public: double getTemp() const { return m_temp; } double getY() const { return m_y; } double getENR() const { return m_enr; } + double getFloor() const { return m_floor; } - static MsgNFMeasurement* create(double frequency, double nf, double temp, double y, double enr) + static MsgNFMeasurement* create(double frequency, double nf, double temp, double y, double enr, double floor) { - return new MsgNFMeasurement(frequency, nf, temp, y, enr); + return new MsgNFMeasurement(frequency, nf, temp, y, enr, floor); } private: @@ -107,14 +108,16 @@ public: double m_temp; // In Kelvin double m_y; // In dB double m_enr; // In dB + double m_floor; // In dBm - MsgNFMeasurement(double frequency, double nf, double temp, double y, double enr) : + MsgNFMeasurement(double frequency, double nf, double temp, double y, double enr, double floor) : Message(), m_frequency(frequency), m_nf(nf), m_temp(temp), m_y(y), - m_enr(enr) + m_enr(enr), + m_floor(floor) { } }; diff --git a/plugins/channelrx/noisefigure/noisefigureenrdialog.cpp b/plugins/channelrx/noisefigure/noisefigureenrdialog.cpp index c979698cc..1200fe989 100644 --- a/plugins/channelrx/noisefigure/noisefigureenrdialog.cpp +++ b/plugins/channelrx/noisefigure/noisefigureenrdialog.cpp @@ -17,11 +17,15 @@ #include +#include + #include "noisefigureenrdialog.h" +#include "util/interpolation.h" NoiseFigureENRDialog::NoiseFigureENRDialog(NoiseFigureSettings *settings, QWidget* parent) : QDialog(parent), m_settings(settings), + m_chart(nullptr), ui(new Ui::NoiseFigureENRDialog) { ui->setupUi(this); @@ -29,6 +33,8 @@ NoiseFigureENRDialog::NoiseFigureENRDialog(NoiseFigureSettings *settings, QWidge for (int i = 0; i < m_settings->m_enr.size(); i++) { addRow( m_settings->m_enr[i]->m_frequency, m_settings->m_enr[i]->m_enr); } + ui->interpolation->setCurrentIndex((int)m_settings->m_interpolation); + plotChart(); } NoiseFigureENRDialog::~NoiseFigureENRDialog() @@ -52,11 +58,13 @@ void NoiseFigureENRDialog::accept() NoiseFigureSettings::ENR *enr = new NoiseFigureSettings::ENR(freqValue, enrValue); m_settings->m_enr.append(enr); } + m_settings->m_interpolation = (NoiseFigureSettings::Interpolation)ui->interpolation->currentIndex(); } void NoiseFigureENRDialog::addRow(double freq, double enr) { ui->enr->setSortingEnabled(false); + ui->enr->blockSignals(true); int row = ui->enr->rowCount(); ui->enr->setRowCount(row + 1); QTableWidgetItem *freqItem = new QTableWidgetItem(); @@ -65,6 +73,7 @@ void NoiseFigureENRDialog::addRow(double freq, double enr) ui->enr->setItem(row, ENR_COL_ENR, enrItem); freqItem->setData(Qt::DisplayRole, freq); enrItem->setData(Qt::DisplayRole, enr); + ui->enr->blockSignals(false); ui->enr->setSortingEnabled(true); } @@ -81,4 +90,125 @@ void NoiseFigureENRDialog::on_deleteRow_clicked() int row = indexList.at(0).row(); ui->enr->removeRow(row); } + plotChart(); +} + +void NoiseFigureENRDialog::on_enr_cellChanged(int row, int column) +{ + (void) row; + (void) column; + plotChart(); +} + +void NoiseFigureENRDialog::on_start_valueChanged(int value) +{ + (void) value; + plotChart(); +} + +void NoiseFigureENRDialog::on_stop_valueChanged(int value) +{ + (void) value; + plotChart(); +} + +void NoiseFigureENRDialog::plotChart() +{ + QChart *oldChart = m_chart; + + m_chart = new QChart(); + + m_chart->layout()->setContentsMargins(0, 0, 0, 0); + m_chart->setMargins(QMargins(1, 1, 1, 1)); + m_chart->setTheme(QChart::ChartThemeDark); + + int size = ui->enr->rowCount(); + + QLineSeries *linearSeries = new QLineSeries(); + QLineSeries *barySeries = new QLineSeries(); + double maxENR = 0.0; + double minENR = 1000.0; + + if (size >= 2) + { + // Sort in to ascending frequency + std::vector> points; + for (int i = 0; i < size; i++) + { + QTableWidgetItem *freqItem = ui->enr->item(i, ENR_COL_FREQ); + QTableWidgetItem *enrItem = ui->enr->item(i, ENR_COL_ENR); + double freq = freqItem->data(Qt::DisplayRole).toDouble(); + double enr = enrItem->data(Qt::DisplayRole).toDouble(); + std::array p = {freq, enr}; + points.push_back(p); + } + sort(points.begin(), points.end()); + + // Copy to vectors for interpolation routines. Is there a better way? + std::vector x(size); + std::vector y(size); + for (int i = 0; i < size; i++) + { + x[i] = points[i][0]; + y[i] = points[i][1]; + } + int order = size - 1; + boost::math::barycentric_rational interpolant(std::move(x), std::move(y), order); + + x.resize(size); + y.resize(size); + for (int i = 0; i < size; i++) + { + x[i] = points[i][0]; + y[i] = points[i][1]; + } + + // Plot interpolated ENR using both methods for comparison + linearSeries->setName("Linear"); + barySeries->setName("Barycentric"); + double startFreq = ui->start->value(); + double stopFreq = ui->stop->value(); + double step = (stopFreq - startFreq) / 50.0; + for (double freq = startFreq; freq < stopFreq; freq += step) + { + double enr; + + // Linear interpolation + enr = Interpolation::linear(x.begin(), x.end(), y.begin(), freq); + linearSeries->append(freq, enr); + minENR = std::min(minENR, enr); + maxENR = std::max(maxENR, enr); + + // Barycentric interpolation + enr = interpolant(freq); + barySeries->append(freq, enr); + minENR = std::min(minENR, enr); + maxENR = std::max(maxENR, enr); + } + } + + QValueAxis *xAxis = new QValueAxis(); + QValueAxis *yAxis = new QValueAxis(); + + m_chart->addAxis(xAxis, Qt::AlignBottom); + m_chart->addAxis(yAxis, Qt::AlignLeft); + + xAxis->setTitleText("Frequency (MHz)"); + yAxis->setTitleText("ENR (dB)"); + + m_chart->addSeries(linearSeries); + + m_chart->addSeries(barySeries); + + linearSeries->attachAxis(xAxis); + linearSeries->attachAxis(yAxis); + + barySeries->attachAxis(xAxis); + barySeries->attachAxis(yAxis); + + yAxis->setRange(std::floor(minENR * 10.0)/10.0, std::ceil(maxENR * 10.0)/10.0); + + ui->chart->setChart(m_chart); + + delete oldChart; } diff --git a/plugins/channelrx/noisefigure/noisefigureenrdialog.h b/plugins/channelrx/noisefigure/noisefigureenrdialog.h index 84d157b31..7a2fae08b 100644 --- a/plugins/channelrx/noisefigure/noisefigureenrdialog.h +++ b/plugins/channelrx/noisefigure/noisefigureenrdialog.h @@ -18,9 +18,13 @@ #ifndef INCLUDE_NOISEFIGUREENRDIALOG_H #define INCLUDE_NOISEFIGUREENRDIALOG_H +#include + #include "ui_noisefigureenrdialog.h" #include "noisefiguresettings.h" +using namespace QtCharts; + class NoiseFigureENRDialog : public QDialog { Q_OBJECT @@ -39,10 +43,15 @@ private slots: void accept(); void on_addRow_clicked(); void on_deleteRow_clicked(); + void on_enr_cellChanged(int row, int column); + void on_start_valueChanged(int value); + void on_stop_valueChanged(int value); private: + QChart *m_chart; Ui::NoiseFigureENRDialog* ui; void addRow(double freq, double enr); + void plotChart(); }; #endif // INCLUDE_NOISEFIGUREENRDIALOG_H diff --git a/plugins/channelrx/noisefigure/noisefigureenrdialog.ui b/plugins/channelrx/noisefigure/noisefigureenrdialog.ui index a0986ffe2..09c0a5e25 100644 --- a/plugins/channelrx/noisefigure/noisefigureenrdialog.ui +++ b/plugins/channelrx/noisefigure/noisefigureenrdialog.ui @@ -6,8 +6,8 @@ 0 0 - 284 - 496 + 347 + 638 @@ -28,7 +28,7 @@ 0 - + @@ -102,6 +102,130 @@ + + + + + 0 + 0 + + + + + 200 + 200 + + + + Chart + + + + 2 + + + 3 + + + 3 + + + 3 + + + 3 + + + + + + + Interpolation + + + + + + + + 0 + 0 + + + + Select interpolation method to use + + + + Linear + + + + + Barycentric + + + + + + + + + + + 200 + 200 + + + + Comparison of interpolated ENR using different methods + + + + + + + + + Start (MHz) + + + + + + + Start frequency to plot in MHz + + + 20000 + + + + + + + Stop (MHz) + + + + + + + Stop frequency to plot in MHz + + + 20000 + + + 3000 + + + + + + + + @@ -117,6 +241,22 @@ + + + QChartView + QGraphicsView +
QtCharts
+
+
+ + enr + addRow + deleteRow + interpolation + chart + start + stop + diff --git a/plugins/channelrx/noisefigure/noisefiguregui.cpp b/plugins/channelrx/noisefigure/noisefiguregui.cpp index ba3578768..fd9b7fbe5 100644 --- a/plugins/channelrx/noisefigure/noisefiguregui.cpp +++ b/plugins/channelrx/noisefigure/noisefiguregui.cpp @@ -71,6 +71,7 @@ void NoiseFigureGUI::resizeTable() ui->results->setItem(row, RESULTS_COL_TEMP, new QTableWidgetItem("10000")); ui->results->setItem(row, RESULTS_COL_Y, new QTableWidgetItem("10.00")); ui->results->setItem(row, RESULTS_COL_ENR, new QTableWidgetItem("10.00")); + ui->results->setItem(row, RESULTS_COL_FLOOR, new QTableWidgetItem("-174.00")); ui->results->resizeColumnsToContents(); ui->results->removeRow(row); } @@ -86,17 +87,20 @@ void NoiseFigureGUI::measurementReceived(NoiseFigure::MsgNFMeasurement& report) QTableWidgetItem *tempItem = new QTableWidgetItem(); QTableWidgetItem *yItem = new QTableWidgetItem(); QTableWidgetItem *enrItem = new QTableWidgetItem(); + QTableWidgetItem *floorItem = new QTableWidgetItem(); ui->results->setItem(row, RESULTS_COL_FREQ, freqItem); ui->results->setItem(row, RESULTS_COL_NF, nfItem); ui->results->setItem(row, RESULTS_COL_TEMP, tempItem); ui->results->setItem(row, RESULTS_COL_Y, yItem); ui->results->setItem(row, RESULTS_COL_ENR, enrItem); + ui->results->setItem(row, RESULTS_COL_FLOOR, floorItem); freqItem->setData(Qt::DisplayRole, report.getFrequency()); nfItem->setData(Qt::DisplayRole, report.getNF()); tempItem->setData(Qt::DisplayRole, report.getTemp()); yItem->setData(Qt::DisplayRole, report.getY()); enrItem->setData(Qt::DisplayRole, report.getENR()); + floorItem->setData(Qt::DisplayRole, report.getFloor()); ui->results->setSortingEnabled(true); plotChart(); @@ -647,6 +651,7 @@ NoiseFigureGUI::NoiseFigureGUI(PluginAPI* pluginAPI, DeviceUISet *deviceUISet, B ui->results->setItemDelegateForColumn(RESULTS_COL_TEMP, new DecimalDelegate(0)); ui->results->setItemDelegateForColumn(RESULTS_COL_Y, new DecimalDelegate(2)); ui->results->setItemDelegateForColumn(RESULTS_COL_ENR, new DecimalDelegate(2)); + ui->results->setItemDelegateForColumn(RESULTS_COL_FLOOR, new DecimalDelegate(1)); displaySettings(); applySettings(true); diff --git a/plugins/channelrx/noisefigure/noisefiguregui.h b/plugins/channelrx/noisefigure/noisefiguregui.h index 337e39f32..f614aa611 100644 --- a/plugins/channelrx/noisefigure/noisefiguregui.h +++ b/plugins/channelrx/noisefigure/noisefiguregui.h @@ -106,7 +106,8 @@ private: RESULTS_COL_NF, RESULTS_COL_TEMP, RESULTS_COL_Y, - RESULTS_COL_ENR + RESULTS_COL_ENR, + RESULTS_COL_FLOOR }; private slots: diff --git a/plugins/channelrx/noisefigure/noisefiguregui.ui b/plugins/channelrx/noisefigure/noisefiguregui.ui index 4898fad03..7e1b04288 100644 --- a/plugins/channelrx/noisefigure/noisefiguregui.ui +++ b/plugins/channelrx/noisefigure/noisefiguregui.ui @@ -696,6 +696,14 @@ Excess noise ratio of the noise source in dB + + + Floor (dBm) + + + Noise floor in dBm for 1Hz bandwidth at 290K + + @@ -767,6 +775,11 @@ ENR (dB) + + + Noise floor (dBm) + + @@ -849,15 +862,20 @@ deltaFrequency fftSize fftCount + frequencySpec start stop steps + step + frequencies startStop saveResults clearResults enr control chartSelect + openReference + clearReference chart results diff --git a/plugins/channelrx/noisefigure/noisefiguresettings.cpp b/plugins/channelrx/noisefigure/noisefiguresettings.cpp index c85850181..35d3ac513 100644 --- a/plugins/channelrx/noisefigure/noisefiguresettings.cpp +++ b/plugins/channelrx/noisefigure/noisefiguresettings.cpp @@ -54,6 +54,7 @@ void NoiseFigureSettings::resetToDefaults() m_powerDelay = 0.5; qDeleteAll(m_enr); m_enr << new ENR(1000.0, 15.0); + m_interpolation = LINEAR; m_rgbColor = QColor(0, 100, 200).rgb(); m_title = "Noise Figure"; m_streamIndex = 0; @@ -106,6 +107,8 @@ QByteArray NoiseFigureSettings::serialize() const s.writeU32(24, m_reverseAPIDeviceIndex); s.writeU32(25, m_reverseAPIChannelIndex); + s.writeS32(26, (int)m_interpolation); + for (int i = 0; i < NOISEFIGURE_COLUMNS; i++) { s.writeS32(100 + i, m_resultsColumnIndexes[i]); } @@ -174,6 +177,8 @@ bool NoiseFigureSettings::deserialize(const QByteArray& data) d.readU32(25, &utmp, 0); m_reverseAPIChannelIndex = utmp > 99 ? 99 : utmp; + d.readS32(26, (int*)&m_interpolation, LINEAR); + for (int i = 0; i < NOISEFIGURE_COLUMNS; i++) { d.readS32(100 + i, &m_resultsColumnIndexes[i], i); } diff --git a/plugins/channelrx/noisefigure/noisefiguresettings.h b/plugins/channelrx/noisefigure/noisefiguresettings.h index a7eea22c3..a569c7e17 100644 --- a/plugins/channelrx/noisefigure/noisefiguresettings.h +++ b/plugins/channelrx/noisefigure/noisefiguresettings.h @@ -28,7 +28,7 @@ class Serializable; // Number of columns in the table -#define NOISEFIGURE_COLUMNS 5 +#define NOISEFIGURE_COLUMNS 6 struct NoiseFigureSettings { @@ -69,6 +69,10 @@ struct NoiseFigureSettings double m_powerDelay; // m_enr; + enum Interpolation { + LINEAR, + BARYCENTRIC + } m_interpolation; quint32 m_rgbColor; QString m_title; diff --git a/plugins/channelrx/noisefigure/readme.md b/plugins/channelrx/noisefigure/readme.md index 5aae88682..066b6cbce 100644 --- a/plugins/channelrx/noisefigure/readme.md +++ b/plugins/channelrx/noisefigure/readme.md @@ -34,7 +34,7 @@ Determines the size (number of points) of the FFT used to measure noise power. O

5: BW

-Displays the measurement bandwidth in kHz as determined by (2). +Displays the measurement bandwidth in Hz as determined by the FFT size (4) and the sample rate.

6: FFTs to average

@@ -64,7 +64,7 @@ Clears the current results from the table and chart. Opens the ENR dialog to allow entering the Excess Noise Ratios (ENRs) for noise source. ENR specifies the difference in noise source power output in dB from when the source is powered off compared to when it is powered on. This typically varies with frequency, so values should be entered for each calibrated frequency. When a measurement is attempted at a frequency for which a value is not specified, an interpolated value will be used. -Barycentric Rational Interpolation is used. +You can choose between linear and barycentric rational interpolation, and the difference between the two is shown in the chart. ![Noise figure ENR dialog](../../../doc/img/NoiseFigure_plugin_enr.png) @@ -87,6 +87,7 @@ Displays measurement results. * T - Calculated noise temperature in Kelvin with a reference temperature of 290K. * Y - Measured Y factor in dB (difference in measured power when the noise source is on and off). * ENR - Excess noise factor of the noise source in dB. +* Floor - Noise floor in dBm assuming 1Hz bandwidth at 290K. ![Noise figure results table](../../../doc/img/NoiseFigure_plugin_results.png) diff --git a/sdrbase/CMakeLists.txt b/sdrbase/CMakeLists.txt index c07266312..01c14871c 100644 --- a/sdrbase/CMakeLists.txt +++ b/sdrbase/CMakeLists.txt @@ -191,6 +191,7 @@ set(sdrbase_SOURCES util/fits.cpp util/golay2312.cpp util/httpdownloadmanager.cpp + util/interpolation.cpp util/lfsr.cpp util/maidenhead.cpp util/message.cpp @@ -396,6 +397,7 @@ set(sdrbase_HEADERS util/httpdownloadmanager.h util/incrementalarray.h util/incrementalvector.h + util/interpolation.h util/lfsr.h util/maidenhead.h util/message.h diff --git a/sdrbase/util/interpolation.cpp b/sdrbase/util/interpolation.cpp new file mode 100644 index 000000000..69b194daa --- /dev/null +++ b/sdrbase/util/interpolation.cpp @@ -0,0 +1,18 @@ +/////////////////////////////////////////////////////////////////////////////////// +// Copyright (C) 2021 Jon Beniston, M7RCE // +// // +// This program is free software; you can redistribute it and/or modify // +// it under the terms of the GNU General Public License as published by // +// the Free Software Foundation as version 3 of the License, or // +// (at your option) any later version. // +// // +// This program is distributed in the hope that it will be useful, // +// but WITHOUT ANY WARRANTY; without even the implied warranty of // +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the // +// GNU General Public License V3 for more details. // +// // +// You should have received a copy of the GNU General Public License // +// along with this program. If not, see . // +/////////////////////////////////////////////////////////////////////////////////// + +#include "interpolation.h" diff --git a/sdrbase/util/interpolation.h b/sdrbase/util/interpolation.h new file mode 100644 index 000000000..0d0e8285d --- /dev/null +++ b/sdrbase/util/interpolation.h @@ -0,0 +1,99 @@ +/////////////////////////////////////////////////////////////////////////////////// +// Copyright (C) 2021 Jon Beniston, M7RCE // +// // +// This program is free software; you can redistribute it and/or modify // +// it under the terms of the GNU General Public License as published by // +// the Free Software Foundation as version 3 of the License, or // +// (at your option) any later version. // +// // +// This program is distributed in the hope that it will be useful, // +// but WITHOUT ANY WARRANTY; without even the implied warranty of // +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the // +// GNU General Public License V3 for more details. // +// // +// You should have received a copy of the GNU General Public License // +// along with this program. If not, see . // +/////////////////////////////////////////////////////////////////////////////////// + +#ifndef INCLUDE_INTERPOLATION_H +#define INCLUDE_INTERPOLATION_H + +#include "export.h" + +class SDRBASE_API Interpolation { +public: + + // Linear interpolation/extrapolation + // Assumes x is sorted in increasing order + template + static T linear(Iter xBegin, Iter xEnd, Iter yBegin, T x) + { + // Find first point in x that target is bigger than + int i = 0; + while (xBegin != xEnd) + { + if (x < *xBegin) { + break; + } + ++xBegin; + i++; + } + T x0, x1, y0, y1; + if (i == 0) { + // Extrapolate left + x0 = *xBegin; + ++xBegin; + x1 = *xBegin; + y0 = *yBegin; + ++yBegin; + y1 = *yBegin; + + return extrapolate(x0, y0, x1, y1, x); + } else if (xBegin > xEnd) { + // Extrapolate right + Iter xCur = xBegin; + std::advance(xCur, i - 2); + x0 = *xCur; + ++xCur; + x1 = *xCur; + + Iter yCur = yBegin; + std::advance(yCur, i - 2); + y0 = *yCur; + ++yCur; + y1 = *yCur; + + return extrapolate(x0, y0, x1, y1, x); + } else { + // Interpolate + x1 = *xBegin; + --xBegin; + x0 = *xBegin; + + Iter yCur = yBegin; + std::advance(yCur, i - 1); + y0 = *yCur; + ++yCur; + y1 = *yCur; + + return interpolate(x0, y0, x1, y1, x); + } + } + + // Linear extrapolation + template + static T extrapolate(T x0, T y0, T x1, T y1, T x) + { + return y0 + ((x-x0)/(x1-x0)) * (y1-y0); + } + + // Linear interpolation + template + static T interpolate(T x0, T y0, T x1, T y1, T x) + { + return (y0*(x1-x) + y1*(x-x0)) / (x1-x0); + } + +}; + +#endif // INCLUDE_INTERPOLATION_H