MyTetra Share
Делитесь знаниями!
Интерполяция полиномом 3 степени
Время создания: 16.09.2019 16:59
Раздел: Computer - Programming - C

//Интерполяция кубическими сплайном третьего порядка

//Вариант для вычисления средних промежуточных значений по оси времени

//Требуется построить интерполирующую функцию F(x), такую, что она принимает в указаных точках те же значения, т.е. F(x0) = y0, F(x1) = y1, ... F(xN) = yN.

//Геометрически это значит, что нужно найти кривую y = F(x) определенного типа, проходящую через систему заданных точек. В такой общей постановке задача может иметь бесчисленное множество решений или совсем не иметь решений. В случае интерполяции сплайном кривая F(x) состоит из кусочков, а именно, на каждом из отрезков [xk-1; xk] функция F(x) является кубическим полиномом

//Fk(x) = ak + bk(x-xk) + ck(x-xk)2 + dk(x-xk)3

//F = F1 на интервале [x0, x1]

//F = F2 на интервале [x1, x2]

//...

//F = FN на интервале [xN-1, xN]

//При этом, на каждом из отрезков [xk-1; xk] коэффициенты полинома ak, bk, ck, dk разные.


int DisplayPixelSmallTimeGap; //Расстояние между промежуточными позициями сетки в пикселях дисплея при прокрутке

int DisplayPixelShift; //Сдвиг временной шкалы относительно нулевой координаты времени в пикселях шкалы




float *BufProbe = new float [2048]; //Бфуер входного сгнала

float *BufProbe05 = new float [2048]; //Буфер интерполированыых промежуточных отсчетов


float *BufProbeSum = new float [2048]; //Суммарный буфер для исходных и интерполированыых промежуточных отсчетов


float P; //Число ПИ = 3.14159265358979323846;


float *x, *y, *h, *l, *delta, *lambda, *c, *d, *b; //Указатели на массивы в стиле С++ без резервирования памяти

float Hk; //

int N, k; //Внешние переменные и счетчики в циклах.


void allocmatrix() //Декларация функций. Создаем массивы для хранения коэффициентов сплайна и промежуточных данных

//(прогоночных коэффициентов), необходимых для вычисления коэффициентов полиномаю Освобождаем память.

{

//allocate memory for matrixes //Массивы для коэффициентов

x = new float[N+1];

y = new float[N+1];

h = new float[N+1];

l = new float[N+1];

delta = new float[N+1];

lambda = new float[N+1];

c = new float[N+1];

d = new float[N+1];

b = new float[N+1];

}

void freematrix() //Удаление массивов, освобождение памяти

{

delete [] x;

delete [] y;

delete [] h;

delete [] l;

delete [] delta;

delete [] lambda;

delete [] c;

delete [] d;

delete [] b;

}




//==============================================================================



//Реализация процесса интерполяции и отображения полученных результатов в событии "Нажата кнопка BitBtn1"


void __fastcall TForm1::BitBtn1Click(TObject *Sender)

{


DisplayPixelSmallTimeGap = 8; //Расстояние между промежуточными позициями сетки в пикселях дисплея при прокрутке

DisplayPixelShift = 0; //Сдвиг временной шкалы относительно нулевой координаты времени в пикселях шкалы



k=0; //Счетчик

N = 256; //Число точек входного массива, содержащего отсчеты исходного сигнала

//Может задаваться любое число точек, 256 - это количество, выбранное

//для отладки программы

P = 3.14159265358979323846;


// BufProbe[х]; //Массив входного сигнала, полученный от внешних устройств, прочитанный

//из файла или сгенерированный программно

//Генерим тестовый сигнал

for (int i=0; i<N; i ++)

{

BufProbe[i] = 200.0*sin((P*(double)i)/10.0)*sin((P*(double)i)/178.0);


}



//Отбражаем сигнал

// Параметр PixelPerPointSignal - этот расстояние между основными отсчетами исходного сигнала в пикселях

// Параметр DisplayPixelShift - это сдвиг сигнала графика сигнала при отображении. Оба параметра рисования - особенность отладочной программы, их можно заменить просто целыми числами

//но при отладке было удобно использовать переменные, поскольку эти дальше параметры приходится менять, и проще это сделать умножив переменную на поправочный кеффициент.


Image1->Canvas->Pen->Color = clLime;

Image1->Canvas->MoveTo(0, Image1->Height/2);

Image1->Canvas->Pen->Style = psSolid;

for (int i= 0; i<2048; i++)

{

Image1->Canvas->LineTo((i*PixelPerPointSignal - DisplayPixelShift), Image1->Height/2 - BufProbe[i]);


}




allocmatrix(); //Создаем массивы прожежуточных данных


//read matrixes a and b from input BufProbe //Заполнение массивов для коэффициентов сплайна

//Всего для интерполяции необходимо 4 коеффициента

//Коеффициент y[х] - это отсчеты входного сигнала

//Интперполирующие Коеффициенты полинома (сплайнап) b[х], c[х] и d[х] вычисляются далее.


for(int i=0; i<N+1; i++)

{

x[i] = i; //Массив номеров отсчетов исходгого сигнала. Номера могут быть какими угодно, но в классическом случае,

//т.е. в нашем, это номера отсчетов массыва. По большому счету их можно генерировать

//в цикле или вычислять в процессе реализации, а не хранить в массиве.


y[i] = BufProbe[i]; //Собственно отсчеты исходного сигнала

}



// Первыми вычисляются массивы коеффициентов h[х] и l[х].

// for(k=1; k<=N; k++)

// {

// h[k] = x[k] - x[k-1];

//

// if(h[k]==0)

// {

// printf("\nError, x[%d]=x[%d]\n", k, k-1);

// return;

// }

//

// l[k] = (y[k] - y[k-1]);///h[k];

// }

// В случае интерполяции с симметричным удвоением точек коеффициент h[х] всегда равен "1" для всех N отсчетов.

//Поскольку потенциально используется для вычислений с дробными числами, необходимо указывать соответствующий тип

//или признак в виде формата для компилятора, например, число "1.0" и результаты вычислений с ним в С++ автоматически преобразуются в операции с плавающей точкой

//В данном случае коэффициент h[х] заменен на переменную Hk = 1.0; типа float.


Hk = 1.0;


//Для Hk = 1.0 вместо h[х] цикл вычисления прогоночного коэффициента l[х] сводится к виду:

//





for(k=1; k<=N; k++)

{


С = (y[k] - y[k-1]);///h[k]; // Прогоночный коэффициент l[х] для всех N отсчетов.

}


//Значения первых отсчетов прогоночных коеффициентов delta[1] и lambda[1]


delta[1] = - 0.25;

lambda[1] = 1.5*(l[2] - l[1])/2.0;


//Вычисление следующих значений прогоночных коеффициентов delta[х] и lambda[х]


for(k=3; k<=N; k++)

{

delta[k-1] = - Hk/(4.0 + delta[k-2]);

lambda[k-1] = (3.0*(l[k] - l[k-1]) - lambda[k-2])/(4.0 + delta[k-2]);



}

//Дальше вычисляем коеффициенты полинома (интерполирующей функции)



//Коеффициенты c[х]

//Крайние значения коэффициентов c[х]


c[0] = 0;

c[N] = 0;


//Остальные значения коэффициентов c[х]

for(k=N; k>=2; k--)

{


c[k-1] = delta[k-1]*c[k] + lambda[k-1];

}


//Коеффициенты d[х] и b[х]


for(k=1; k<=N; k++)

{

d[k] = (c[k] - c[k-1])/3.0;

b[k] = l[k] + (2.0*c[k] + c[k-1])/3.0;


}


//ВСЕ!!! Интерполирующий полином вычислен!



//Теперь интерполируем промежуточые осчеты в исходном сигнале BufProbe[х] (пока они находятся в массиве y[х])

// и помещаем их в специально созданный буфер BufProbe05[s];



// Поскольку самый первый промежуточный отсчет находится между нулевым и первым отсчетом входного сигнала, его значение в

// нулевой ячейке буферного массива всегда равно "0"


BufProbe05[0] = 0;



//Далее вычисляем промежуточные отсчеты интерполированного сигнала, расположенные по середине между отсчетами

//входного сигнала, т.е. находящиеся в массиве исходного сигнала и на расстоянии +0.5 от первого отсчета исходого

//сигнала в интерполируемом промежутке между отсчетами с номером к и номером к+1.

//Результат помещаем в буфер BufProbe05[k].


for(int k=1; k<=N; k++)

{

BufProbe05[k] = y[k] + b[k]*0.5 + c[k]*pow(0.5, 2) + d[k]*pow(0.5, 3);

}



//Далее графически отображаем резельтаты интерполяции для прверки корректности наших вычислений


// Все графическое отбражение здесь и далее показано на компонентах C++ Builder.



//1. Сначала отбражаем исходный сигнал. На компоненте Image1 цветом лайма выводим график исходного сигнала


Image1->Canvas->Pen->Color = clLime;

Image1->Canvas->MoveTo(0, Image1->Height/2);

Image1->Canvas->Pen->Style = psSolid;

for (int i= 0; i<2048; i++)

{

Image1->Canvas->LineTo((i*PixelPerPointSignal - DisplayPixelShift), Image1->Height/2 - BufProbe[i]);


}




//Отображаем отсчеты исходного сигнала в виде верикальных линий цветом лайма

// Параметр PixelPerPointSignal - этот расстояние между основными отсчетами исходного сигнала


//в пикселях при выбранном масштабе по оси времени в отладочной программе, которая использовалась.

//Можно заменить на любое целочисленное значение, масштабируемое программно или просто на удобное целое число.


Image1->Canvas->Pen->Color = clLime;

Image1->Canvas->MoveTo(0, Image1->Height/2);

Image1->Canvas->Pen->Style = psSolid;

for (int i= 0; i<N; i++)

{

Image1->Canvas->MoveTo(i*PixelPerPointSignal - DisplayPixelShift, Image1->Height/2);

Image1->Canvas->LineTo((i*PixelPerPointSignal - DisplayPixelShift), Image1->Height/2 - BufProbe[i]);

}




//Отображаем промежуточные интерполированные отсчеты в виде вертикальных линий красным цветом

//Оцениваем результат на корректность.



Image1->Canvas->Pen->Color = clRed;

Image1->Canvas->MoveTo(0, Image1->Height/2);


for (int i= 0; i<N; i++)

{

Image1->Canvas->MoveTo(i*PixelPerPointSignal + PixelPerPointSignal*0.5 - DisplayPixelShift, Image1->Height/2);

Image1->Canvas->LineTo((i*PixelPerPointSignal + PixelPerPointSignal*0.5 - DisplayPixelShift), Image1->Height/2 - BufProbe05[i]);

}



//Отображаем суммарный график с учетом исходных и интерполированных отсчетов сигнала желтым (clYellow) цветом

//Пока отсчеты для графического отображения берем из двух массивов - исходного BufProbe[х] и массива, содержащего

//интерполированные значения BufProbe05[х]


Image1->Canvas->Pen->Color = clYellow;

Image1->Canvas->MoveTo(0, Image1->Height/2);

Image1->Canvas->Pen->Style = psSolid;


// Выводим график в виде осциллограммы суммарного сигнала поочередно обращаясь к массивам.

//Промежуточные интерполированные точки для отображения

//на графике задаются через параметр PixelPerPointSignal


for (int i= 0; i<N; i++)

{

Image1->Canvas->LineTo(i*PixelPerPointSignal - DisplayPixelShift, Image1->Height/2 - BufProbe[i]);

Image1->Canvas->LineTo((i*PixelPerPointSignal + PixelPerPointSignal*0.5 - DisplayPixelShift), Image1->Height/2 - BufProbe05[i]);

}



//Формируем общий массив фрагмента интерполированного сигнала длиной в N точек исходного.

//Теперь длина общего массива будет в два раза больше для того же отрезка временной шкалы или области графического отбражения.

//Копируем в общий массив BufProbeSum[х] отсчеты.


for (int i= 0; i<(2*N); i = i + 2)

{

BufProbeSum[i] = BufProbe[i/2];

BufProbeSum[i+1] = BufProbe05[i/2];

}


//Отображаем график интерполированного отрезка сигнала используя общий массив BufProbeSum[х].

//Поскольку необходимо на том же участке в том же масштабе отобразить в два раза больше точек, производим

//вывод результатов с двойным шагом, поочередно ваполняя рисование для координат с четными элементами (исходные значения)

//и нечетными элементами (интерполированные значения) массива.


for (int i= 0; i<(2*N); i = i+2)

{

Image1->Canvas->LineTo(i*PixelPerPointSignal*0.5 - DisplayPixelShift, Image1->Height/2 - BufProbeSum[i]);

Image1->Canvas->LineTo(i*PixelPerPointSignal*0.5 + PixelPerPointSignal*0.5 - DisplayPixelShift, Image1->Height/2 - BufProbeSum[i+1]);

}




//Освобождаем память.

freematrix();


}






 
MyTetra Share v.0.53
Яндекс индекс цитирования