Тригонометрическая Интерполяция

Разложение в ряд Фурье
Ещё одним методом интерполяции является метод разложение в ряд Фурье. Суть его состоит в следующем.
Имеется четное количество значений некоторой функции.
Используя разложение Фурье она разлагается в конечную сумму
синусов и косинусов. Затем, зная указанное разложение, можно получить значение исходной функции в любой промежуточной точке. Аналитически это описывается следующим образом:

{S}_{k}(x)=\frac{{a}_{0}}{2}+\sum_{n=1}^{k}({a}_{n}*cos(n*x)+{b}_{n}*sin(n*x))

Сама программа разложения в ряд Фурье функции, заданной четным количеством своих значений приведена ниже.

Далее приводится вариант программы на языке С (компилятор BORLAND C++ 3.1)

Show »

// OPFUR02.CPP
//**********************************************************************
// Программа численного гармонического анализа.
// Тригонометрическое интерполирование.
// Фактически это разложение в ряд фурье (тригонометрическая интерполяция)
// При этом выводятся как исходные значения, так и полученные при интерполяции.
//**********************************************************************
#include <stdio.h>
#include <math.h>
#include <conio.h>

//-------------------------
// Глобальные переменные
//-------------------------
int n=12; // четное количество анализируемых точек исходной функции.
// Исходные значения анализируемой функции
float sic[12]={95,71,55,43,36,31,28,26,25,45,91,102};

//-------------------------
// Прототипы функций
//-------------------------
// Функция ОПФ и всех выводов
void opfur(void);

// Главная программа.
void main(void)
{
 clrscr();
 opfur();
}

void opfur(void)
{
 int k;
 float si,ni,d,dm,msi[6],mco[6],mm[12],s;
 float n1=2*M_PI/n;
// Получение коэффициентов Фурье и вывод их
 puts(" НОМЕР КОСИНУС СИНУС");
 puts("ГАРМОНИКИ СОСТАВЛЯЮЩАЯ СОСТАВЛЯЮЩАЯ");
 // цикл по вычисляемым парам гармоник(SIN и COS)
 for(int m=-1; ++m<=n/2-1;printf("%6d %11.3f %11.3fn",
 m,(mco[m]=d*2/n*(m?1:0.5)),(msi[m]=dm*2/n)))
 for(d=0,dm=0,k=0; k<=n-1; d+=(si=sic[k])*cos(ni=(n1*k*m)),dm+=si*sin(ni),k++);
 // Получение значений интерполяции и сравнеие их с исходными значениями
 for(int i=0;i<n;mm[i++]=s)
 for(s=0,k=0;k<n/2;s+=(mco[k]*cos(d=k*2*M_PI/n*i)+msi[k]*sin(d)),k++);
 // 0 1 2 2 1 0
 puts("n");
 puts(" номер исходная полученная");
 puts("отсчета точка точка");
 for(k=0;k<n;printf("n %3d %8.3f %6.3f",k,sic[k],mm[k]),k++);
}