1
0
mirror of https://github.com/f4exb/sdrangel.git synced 2024-12-23 10:05:46 -05:00

Merge pull request #919 from srcejon/noise_figure_interpolation

Noise figure updates
This commit is contained in:
Edouard Griffiths 2021-06-11 14:15:32 +02:00 committed by GitHub
commit 1b15739d29
No known key found for this signature in database
GPG Key ID: 4AEE18F83AFDEB23
16 changed files with 472 additions and 22 deletions

Binary file not shown.

Before

Width:  |  Height:  |  Size: 17 KiB

After

Width:  |  Height:  |  Size: 30 KiB

Binary file not shown.

Before

Width:  |  Height:  |  Size: 12 KiB

After

Width:  |  Height:  |  Size: 12 KiB

View File

@ -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<double> x(size);
std::vector<double> 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<double> 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<double> 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()));
}

View File

@ -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)
{
}
};

View File

@ -17,11 +17,15 @@
#include <QDebug>
#include <boost/math/interpolators/barycentric_rational.hpp>
#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<std::array<double,2>> 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<double,2> 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<double> x(size);
std::vector<double> 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<double> 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;
}

View File

@ -18,9 +18,13 @@
#ifndef INCLUDE_NOISEFIGUREENRDIALOG_H
#define INCLUDE_NOISEFIGUREENRDIALOG_H
#include <QtCharts>
#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

View File

@ -6,8 +6,8 @@
<rect>
<x>0</x>
<y>0</y>
<width>284</width>
<height>496</height>
<width>347</width>
<height>638</height>
</rect>
</property>
<property name="font">
@ -28,7 +28,7 @@
<verstretch>0</verstretch>
</sizepolicy>
</property>
<layout class="QVBoxLayout" name="verticalLayout_2">
<layout class="QVBoxLayout" name="verticalLayout_3">
<item>
<widget class="QLabel" name="enrLabel">
<property name="text">
@ -102,6 +102,130 @@
</item>
</layout>
</item>
<item>
<widget class="QWidget" name="chartContainer" native="true">
<property name="sizePolicy">
<sizepolicy hsizetype="Expanding" vsizetype="Expanding">
<horstretch>0</horstretch>
<verstretch>0</verstretch>
</sizepolicy>
</property>
<property name="minimumSize">
<size>
<width>200</width>
<height>200</height>
</size>
</property>
<property name="windowTitle">
<string>Chart</string>
</property>
<layout class="QVBoxLayout" name="verticalLayout_2" stretch="0,0,0">
<property name="spacing">
<number>2</number>
</property>
<property name="leftMargin">
<number>3</number>
</property>
<property name="topMargin">
<number>3</number>
</property>
<property name="rightMargin">
<number>3</number>
</property>
<property name="bottomMargin">
<number>3</number>
</property>
<item>
<layout class="QHBoxLayout" name="horizontalLayout">
<item>
<widget class="QLabel" name="interpolationLabel">
<property name="text">
<string>Interpolation</string>
</property>
</widget>
</item>
<item>
<widget class="QComboBox" name="interpolation">
<property name="sizePolicy">
<sizepolicy hsizetype="Expanding" vsizetype="Fixed">
<horstretch>0</horstretch>
<verstretch>0</verstretch>
</sizepolicy>
</property>
<property name="toolTip">
<string>Select interpolation method to use</string>
</property>
<item>
<property name="text">
<string>Linear</string>
</property>
</item>
<item>
<property name="text">
<string>Barycentric</string>
</property>
</item>
</widget>
</item>
</layout>
</item>
<item>
<widget class="QChartView" name="chart">
<property name="minimumSize">
<size>
<width>200</width>
<height>200</height>
</size>
</property>
<property name="toolTip">
<string>Comparison of interpolated ENR using different methods</string>
</property>
</widget>
</item>
<item>
<layout class="QHBoxLayout" name="horizontalLayout_2">
<item>
<widget class="QLabel" name="startLabel">
<property name="text">
<string>Start (MHz)</string>
</property>
</widget>
</item>
<item>
<widget class="QSpinBox" name="start">
<property name="toolTip">
<string>Start frequency to plot in MHz</string>
</property>
<property name="maximum">
<number>20000</number>
</property>
</widget>
</item>
<item>
<widget class="QLabel" name="stopLabel">
<property name="text">
<string>Stop (MHz)</string>
</property>
</widget>
</item>
<item>
<widget class="QSpinBox" name="stop">
<property name="toolTip">
<string>Stop frequency to plot in MHz</string>
</property>
<property name="maximum">
<number>20000</number>
</property>
<property name="value">
<number>3000</number>
</property>
</widget>
</item>
</layout>
</item>
</layout>
</widget>
</item>
</layout>
</widget>
</item>
@ -117,6 +241,22 @@
</item>
</layout>
</widget>
<customwidgets>
<customwidget>
<class>QChartView</class>
<extends>QGraphicsView</extends>
<header>QtCharts</header>
</customwidget>
</customwidgets>
<tabstops>
<tabstop>enr</tabstop>
<tabstop>addRow</tabstop>
<tabstop>deleteRow</tabstop>
<tabstop>interpolation</tabstop>
<tabstop>chart</tabstop>
<tabstop>start</tabstop>
<tabstop>stop</tabstop>
</tabstops>
<resources>
<include location="../../../sdrgui/resources/res.qrc"/>
<include location="../demodapt/icons.qrc"/>

View File

@ -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);

View File

@ -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:

View File

@ -696,6 +696,14 @@
<string>Excess noise ratio of the noise source in dB</string>
</property>
</column>
<column>
<property name="text">
<string>Floor (dBm)</string>
</property>
<property name="toolTip">
<string>Noise floor in dBm for 1Hz bandwidth at 290K</string>
</property>
</column>
</widget>
</item>
</layout>
@ -767,6 +775,11 @@
<string>ENR (dB)</string>
</property>
</item>
<item>
<property name="text">
<string>Noise floor (dBm)</string>
</property>
</item>
</widget>
</item>
<item>
@ -849,15 +862,20 @@
<tabstop>deltaFrequency</tabstop>
<tabstop>fftSize</tabstop>
<tabstop>fftCount</tabstop>
<tabstop>frequencySpec</tabstop>
<tabstop>start</tabstop>
<tabstop>stop</tabstop>
<tabstop>steps</tabstop>
<tabstop>step</tabstop>
<tabstop>frequencies</tabstop>
<tabstop>startStop</tabstop>
<tabstop>saveResults</tabstop>
<tabstop>clearResults</tabstop>
<tabstop>enr</tabstop>
<tabstop>control</tabstop>
<tabstop>chartSelect</tabstop>
<tabstop>openReference</tabstop>
<tabstop>clearReference</tabstop>
<tabstop>chart</tabstop>
<tabstop>results</tabstop>
</tabstops>

View File

@ -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);
}

View File

@ -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; //<! Delay in seconds before starting a measurement
QList<ENR *> m_enr;
enum Interpolation {
LINEAR,
BARYCENTRIC
} m_interpolation;
quint32 m_rgbColor;
QString m_title;

View File

@ -34,7 +34,7 @@ Determines the size (number of points) of the FFT used to measure noise power. O
<h3>5: BW</h3>
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.
<h3>6: FFTs to average</h3>
@ -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)

View File

@ -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

View File

@ -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 <http://www.gnu.org/licenses/>. //
///////////////////////////////////////////////////////////////////////////////////
#include "interpolation.h"

View File

@ -0,0 +1,100 @@
///////////////////////////////////////////////////////////////////////////////////
// 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 <http://www.gnu.org/licenses/>. //
///////////////////////////////////////////////////////////////////////////////////
#ifndef INCLUDE_INTERPOLATION_H
#define INCLUDE_INTERPOLATION_H
#include "export.h"
#include <iterator>
class SDRBASE_API Interpolation {
public:
// Linear interpolation/extrapolation
// Assumes x is sorted in increasing order
template<class Iter, class T>
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<class T>
static T extrapolate(T x0, T y0, T x1, T y1, T x)
{
return y0 + ((x-x0)/(x1-x0)) * (y1-y0);
}
// Linear interpolation
template<class T>
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