Ряд Котельникова

[Smitt&Wesson]

В общем, написал функцию по ряду Котельникова:

// Ряд Котельникова
// dt - шаг
// Y - y(t), T - Время
// M - количество отсчётов
// Возвращает - Y(t)
double HarmonicKotelnik(double &dt, double Y, double T, int M)
{
  double pi, s = 0;
  for(int k = 0; k < M; k++)
  {
    pi = M_PI * (T/dt - k);
    s += Y * k * dt * (sin(pi) / pi);
  }
  dt += dt;
  return s;
}

Здесь, программа обращения к функции (FPrint, подпрограмма вывода на консоль).

  double y = 1, dt = 0.1;
  for(int i = 0; i < 3; i++)
  {
    y = HarmonicKotelnik(dt, y, 10, 5);
    FPrint(y, false);
  }

Сомнения в том, правильно ли я вычисляю y(k dt)? Ниже формула, по которой ведутся вычисления.

[Аватар]

y(k dt) насколько понял это тот же ряд для t= k*dt, а не Y * k * dt

[p51x]

Вы спутала надо не Y * k * dt, а Y(k * dt), т.е. отсчеты

[Smitt&Wesson]

Так вот я и не врублюсь. То ли то, что в скобках просто на Y умножить т.к. выходное значение функции, оно-же входное на следующей итерации. То-ли k * dt по новой вычисляется?

[p51x]

Это отсчеты, если бы это был не дискретный случай, то там стоял бы интеграл и вычислялся бы он... Представьте: у вас отрезок разбит на m+1 частей, dt - шаг, вот и идет суммирование по этим точкам Y(0) Y(dt) Y(2dt) Y(3dt) ...

[Smitt&Wesson]

Сообщение от p51x

Это отсчеты, если бы это был не дискретный случай, то там стоял бы интеграл и вычислялся бы он... Представьте: у вас отрезок разбит на m+1 частей, dt - шаг, вот и идет суммирование по этим точкам Y(0) Y(dt) Y(2dt) Y(3dt) ...

С отсчётами всё понятно. Они в цикле вычисляются от 0 до M. Чё с этими отсчётами в коде делать? Просто с дельтами никогда не сталкивался. Что в формуле обозначает Y(...)? Дискретную математику изучал лет 35 назад. Почти всё из головы выветрилось.

[p51x]

Y(t) - функция, Y(0) Y(dt) Y(2dt) Y(3dt) ... ее значение в определенных точках

[Pavia]

Они в цикле вычисляются от 0 до M.

Отчёты не вычисляются. Они измеряются микрофона и передаются в виде массива в функцию.

Просто с дельтами никогда не сталкивался.

dt -это константа она в формуле не изменяется.
Слева у вас в формуле стоит аналоговый, т.е непрерывный сигнал.
А справа дискретный, т.е. значения в узлах сетки.

// fd - частота дискретизации в Гц.
function InterpolationKotelikov(const t:Real; const fd:Real; y:TArrayReal):Real; Overload;
var
 dt:Real;
begin
dt:=1/df;
for k:=0 to m do
  Result:=Result+y[k]*Sinc(Pi*(t/dt-k));  //Используем sinc, так как автор формулы не позаботился о делении на 0
end;

// Sin(x)/x
function  Sinc(x:Extended):Extended;
begin
If x<>0 then Result:=Sin(x)/(x) else Result:=1;
end;

Почему y[k], а не y[k*dt]? Не вижу надобности. Но если вы хотите, то можете добавить. Но тогда надо брать только кратные 2^n или заменить y[k], на виртуальную функцию.

type
 TFunc = function (t:Real):Real;

function InterpolationKotelikov(const t:Real; const fd:Real; y:TFunc):Real; Overload;
var
 dt:Real;
begin
// fd - частота дискретизации в Гц.
dt:=1/df;
for k:=0 to m do
  Result:=Result+y(k*dt)*Sinc(Pi*(t/dt-k));
end;

Учебник для справки: А.Б.Сергиенко_Цифровая обработка сигналов

[Smitt&Wesson]

Сообщение от Pavia

Отчёты не вычисляются. Они измеряются микрофона и передаются в виде массива в функцию.

dt -это константа она в формуле не изменяется.
Слева у вас в формуле стоит аналоговый, т.е непрерывный сигнал.
А справа дискретный, т.е. значения в узлах сетки.

// fd - частота дискретизации в Гц.

Почему y[k], а не y[k*dt]? Не вижу надобности. Но если вы хотите, то можете добавить. Но тогда надо брать только кратные 2^n или заменить y[k], на виртуальную функцию.

// fd - частота дискретизации в Гц.
dt:=1/df;

Учебник для справки: А.Б.Сергиенко_Цифровая обработка сигналов

Pavia, спасибо. Теперь понятно. Задаём частоту дескритизации. 1/f получаем период в секундах. Это и есть dt, которая неизменна, пока частота дескритизации не изменилась. Дальше всё понятно.

[Smitt&Wesson]

Вот доработанная функция:

// F - частота дескритизации
// M - количество отсчётов
// T - Время отсчёта
// Возвращает - Y(t)
double HarmonicKotelnik(double *Y, double T, double F, int M)
{
  double pi, dt = 1 /F, s = 0, sp;

  for(int k = 0; k < M; k++)
  {
    pi = M_PI * (T/dt - k);
    if(pi != 0)
      s += Y[k] * (sin(pi) / pi);
  }
  return s;
}

Вызов функции:
y = HarmonicKotelnik(Y, N, 1000, N);

Вот ряд, который записывается в массив Y - 0 2 4 6 8 6 4 2 0 -2 -4 -6 -8.
Частота дискретизации - 1000 Гц
Количество отсчётов и время совпадают с длинной массива Y.

Почему функция выдаёт такой результат - 5.182505423e-16 как его интерпретировать? Как спектральную плотность? Но почему -16-я степень?