C++实现整数和复数的一维线性插值interp1
导言
最近在进行Qt开发,涉及大量的matlab转C的工作,其中包括插值滤波等,这里对matlab的interp1函数进行了研究并用C++进行了重写。
经对比,结果与matlab的interp1()函数生成的插值结果一致。(如果数据不在定义域,则按照线性插值进行插值【与matlab的定义域外插值策略处理结果一致】)
函数需求分析
这里我对工作中调用matlab的interp1()函数的实例为例子进行讲解:
- 原matlab程序 :
// Ti:已有的时刻的值
// To:插入的时刻的值
// I_Id: 已有的时刻对应的y值的实部
// Q_Id:已有的时刻对应的y值的虚部
outD=interp1(Ti,I_Id,To,'linear','extrap')+1j*interp1(Ti,Q_Id,To,'linear','extrap');
-
interp1()生成的结果
-
C++版interp1()输入输出参数详情:
参数名 | 介绍 |
---|---|
Ti | 已有的点的x坐标的值 |
To | 插入的点的X坐标的值 |
I_Id | 已有的点实部的y坐标的值 |
Q_Id | 已有的点虚部的y坐标的值 |
Y | 插入的点的实部和虚部的Y坐标的值以复数的形式返回 |
compoundNumber interp1(vector<double> Ti, vector<double> To, vector<double> I_Id, vector<double> Q_Id)
/*
参数介绍
输入:(数据类型都为double类型的vector数组)
请注意:我这里的Ti和To都是排好序的
Ti: 已有的点的x坐标的值
To: 插入的点的X坐标的值
I_Id:已有的点实部的y坐标的值
Q_Id:已有的点虚部的y坐标的值
公式:Y[i] = y[i-1] + (y[i] - y[i-1]) * (X[i] - x[i-1])
输出:插入的点的实部和虚部的Y坐标的值以复数的形式返回
数据类型为自定义复数,
实部虚部都是一个double类型的vector数组
*/
- C++版interp函数生成结果 :结果表明,生成结果与matlab生成结果一致。
源码
那么废话不多说,直接上代码
#include <vector>
#include <iostream>
using namespace std;
const struct comDouble
{
vector<double> real;
vector<double> imag;
double time;
};
//复数插值
comDouble Cominterp1(vector<double> Ti, vector<double> To, vector<INT16> I_Id, vector<INT16> Q_Id) {
/*
参数介绍
输入:(数据类型都为double类型的vector数组)
请注意:我这里的Ti和To都是排好序的
Ti: 已有的点的x坐标的值
To: 插入的点的X坐标的值
I_Id:已有的点实部的y坐标的值
Q_Id:已有的点虚部的y坐标的值
公式:Y[i] = y[i-1] + (y[i] - y[i-1]) * (X[i] - x[i-1])
输出:插入的点的实部和虚部的Y坐标的值以复数的形式返回
数据类型为自定义复数,
实部虚部都是一个double类型的vector数组
*/
comDouble Y;
if (Ti.size() != I_Id.size() || I_Id.size() != Q_Id.size())
cout << "采样点个数和信号时刻数量不一致或实虚部数量不一致!" << endl;
double maxTi = Ti[Ti.size() - 1];
double minTi = Ti[0];
double maxTo = To[To.size() - 1];
double minTo = To[0];
double theK = 0;
int before = 1;
double real, imag;
for (int i = 0; i < To.size(); i++) {
待插入点在定义域左边
//if (To[i] <= minTi) {
// theK = (I_Id[1] - I_Id[0]) / (Ti[1] - Ti[0]);
// real = theK * (To[i] - Ti[0]) + I_Id[0];
// imag = theK * (To[i] - Ti[0]) + I_Id[0];
// Y.real.push_back(real);
// Y.imag.push_back(imag);
// break;
//}
//定义域右边
if (To[i] >= maxTi) {
double theRealK = (I_Id[I_Id.size() - 1] - I_Id[I_Id.size() - 2]) / (Ti[Ti.size() - 1] - Ti[Ti.size() - 2]);
double theImagK = (Q_Id[Q_Id.size() - 1] - Q_Id[Q_Id.size() - 2]) / (Ti[Ti.size() - 1] - Ti[Ti.size() - 2]);
real = theRealK * (To[i] - Ti[Ti.size() - 1]) + I_Id[I_Id.size() - 1];
imag = theImagK * (To[i] - Ti[Ti.size() - 1]) + Q_Id[Q_Id.size() - 1];
Y.real.push_back(real);
Y.imag.push_back(imag);
}
else {
for (int j = before; j < Ti.size(); j++) {
if (To[i] <= Ti[j]) {
real = (I_Id[j] - I_Id[j - 1]) / (Ti[j] - Ti[j - 1]) * (To[i] - Ti[j - 1]) + I_Id[j - 1];
imag = (Q_Id[j] - Q_Id[j - 1]) / (Ti[j] - Ti[j - 1]) * (To[i] - Ti[j - 1]) + Q_Id[j - 1];
Y.real.push_back(real);
Y.imag.push_back(imag);
before = j;
break;
}
}
}
}
return Y;
}
#include <vector>
#include <iostream>
using namespace std;
const struct comDouble
{
vector<double> real;
vector<double> imag;
double time;
};
//整数插值
vector<double> interp1(vector<double> Ti, vector<double> I_Id, vector<double> To)
{
/*
参数介绍
输入:(数据类型都为double类型的vector数组)
Ti: 已有的点的x坐标的值
To: 插入的点的X坐标的值
I_Id:已有的点实部的y坐标的值
公式:Y[i] = y[i-1] + (y[i] - y[i-1]) * (x[i] - x[i-1])
*/
vector<double> Y;
double maxTi = Ti[Ti.size() - 1];
//double minTi = Ti[0];
//double maxTo = To[To.size() - 1];
//double minTo = To[0];
double theK = 0;
double real = 0.0;
int before = 1;
for (int i = 0; i < To.size(); i++) {
待插入点在定义域左边
//if (To[i] <= minTi) {
// theK = (I_Id[1] - I_Id[0]) / (Ti[1] - Ti[0]);
// real = theK * (To[i] - Ti[0]) + I_Id[0];
// Y.push_back(real);
//}
//定义域右边
if (To[i] >= maxTi) {
theK = (I_Id[I_Id.size() - 1] - I_Id[I_Id.size() - 2]) / (Ti[Ti.size() - 1] - Ti[Ti.size() - 2]);
real = theK * (To[i] - Ti[Ti.size() - 1]) + I_Id[I_Id.size() - 1];
Y.push_back(real);
}
else {
for (int j = before; j < Ti.size(); j++) {
if (To[i] <= Ti[j]) {
real = (I_Id[j] - I_Id[j - 1]) / (Ti[j] - Ti[j - 1]) * (To[i] - Ti[j - 1]) + I_Id[j - 1];
Y.push_back(real);
before = j;
break;
}
}
}
}
return Y;
}
性能提升(vector转用动态数组,这里以复数来举例)
#include <vector>
#include <iostream>
using namespace std;
struct comDoubleC
{
double *real;
double *imag;
double time;
double Rsize;
double Isize;
comDoubleC(int reallen, int imaglen) {
this->real = new double[reallen];
this->imag = new double[imaglen];
this->time = 0;
this->Rsize = 0;
this->Isize = 0;
}
};
//复数插值
comDoubleC Cominterp1(
double Ti[],
int Tilen,
double To[],
int Tolen,
INT16 I_Id[],
INT16 Q_Id[],
comDoubleC bh_IQ) {
/*
参数介绍
输入:(数据类型都为double动态数组)
请注意:我这里的Ti和To都是排好序的
Ti: 已有的点的x坐标的值
To: 插入的点的X坐标的值
I_Id:已有的点实部的y坐标的值
Q_Id:已有的点虚部的y坐标的值
bh_IQ:结果---在函数外管理内存
公式:Y[i] = y[i-1] + (y[i] - y[i-1]) * (X[i] - x[i-1])
输出:插入的点的实部和虚部的Y坐标的值以复数的形式返回
数据类型为自定义复数,
实部虚部都是一个double类型的动态数组
*/
double real;
double imag;
double theRealK;
double theImagK;
double maxTi = Ti[Tilen - 1];
int before = 1;
for (int i = 0; i < Tolen; i++) {
//待插入点在定义域右边
if (To[i] >= maxTi) {
theRealK = (I_Id[Tilen - 1] - I_Id[Tilen - 2]) / (Ti[Tilen - 1] - Ti[Tilen - 2]);
theImagK = (Q_Id[Tilen - 1] - Q_Id[Tilen - 2]) / (Ti[Tilen - 1] - Ti[Tilen - 2]);
real = theRealK * (To[i] - Ti[Tilen - 1]) + I_Id[Tilen - 1];
imag = theImagK * (To[i] - Ti[Tilen - 1]) + Q_Id[Tilen - 1];
bh_IQ.real[i] = real;
bh_IQ.imag[i] = imag;
}
else {
for (int j = before; j < Tilen; j++) {
if (To[i] <= Ti[j]) {
real = (I_Id[j] - I_Id[j - 1]) / (Ti[j] - Ti[j - 1]) * (To[i] - Ti[j - 1]) + I_Id[j - 1];
imag = (Q_Id[j] - Q_Id[j - 1]) / (Ti[j] - Ti[j - 1]) * (To[i] - Ti[j - 1]) + Q_Id[j - 1];
bh_IQ.real[i] = real;
bh_IQ.imag[i] = imag;
before = j;
break;
}
}
}
}
return bh_IQ;
}
注意
该文章仅个人学习使用,欢迎大家一起交流学习