Разложение функции в ряд Котельникова

[ivantozavrus]

Согласно теореме Котельникова функцию с ограниченным спектром можно разложить в ряд следующего вида:

$$f(x) = \sum_{k=-\infty}^{\infty} f(k \Delta) \frac{\sin \biggr(\frac{\pi}{\Delta} (x - k \Delta)\biggr)}{\frac{\pi}{\Delta} (x - k \Delta)},$$

где $\Delta \leq \frac{1}{2 F}, $ $F - $ интервал частот, на котором спектр функции отличен от нуля.

Задание: разложить в такой ряд функцию $sinc(\pi x)$. Для начала нужно выполнить преобразование Фурье, чтобы найти интервал $F.$ Я выполнил fft. На графике оставил только часть с положительными значениями частоты. По графику определил, что спектр становится нулевым после $f = \frac{\pi}{2.}$ Значит, $F = \pi, \quad \Delta = \frac{1}{2 \pi}.$
Но когда я выполнил расчеты и построил график для получившегося ряда, никакого точного восстановления функции не получилось. Прилагаю текст программы в Matlab и сами графики (для fft и разложения в ряд Котельникова соответственно).

Может быть, я не правильно понимаю теорему или как-то не так написал программу? Буду рад, если кто-нибудь сможет помочь.

clear

x = -5:0.01:5;
delta = 1/2/pi;

x_k = -100:delta:100;
f_k = sinc(pi*x_k);
n = length(x_k) - 1;
f = zeros(1, length(x));

for i = 1:length(x)
c = 1;
for j = -n/2:n/2
f(i) = f(i) + f_k(c) * sinc(pi/delta * (x(i) - x_k(c)));
c = c + 1;
end
end

syms t
fplot(sinc(pi*t));
grid on
hold on
plot(x, f);
legend('sinc', 'series');

[ozheredov]

ivantozavrus в сообщении #1616305 писал(а):

Задание: разложить в такой ряд функцию $sinc(\pi x)$

Извиняюсь, а это случайно не нулевой ($k=0$) член Вашей верхней формулы при $\Delta = 1$?

$$f(0) = \frac{\sin0}{0} = 1$$

[ivantozavrus]

ozheredov в сообщении #1616309 писал(а):

Извиняюсь, а это случайно не нулевой ($k=0$) член Вашей верхней формулы при $\Delta = 1$?

$$f(0) = \frac{\sin0}{0} = 1$$

Получается, что так. Но разве к функции $sinc(\pi x)$ теорема не применима? Вроде бы, кроме ограниченного спектра, никаких других требований особо и нет. То есть она так же должна раскладываться по формуле для ряда Котельникова.

[ihq.pl]

Ряд, который у вас выписан, попробуйте правильно промасштабировать

[ivantozavrus]

ihq.pl в сообщении #1616375 писал(а):

Ряд, который у вас выписан, попробуйте правильно промасштабировать

Что в данном случае значит промасштабировать?

[ihq.pl]

Изменить масштаб переменной.

[ivantozavrus]

ihq.pl в сообщении #1616380 писал(а):

Изменить масштаб переменной.

Извините, но я пока совсем не понимаю по какому принципу это нужно сделать и зачем

[ihq.pl]

ivantozavrus
Это просто звучит страшно. Можно преобразовать ваше выражение так, чтобы справа оказались обычные синки (не масштабированные), про сдвиги которых известно, что они ортонормальные. Нужно ли - другой вопрос, но этим можно воспользоваться для того, чтобы увидеть/проверить одну очевидную вещь, а заодно потренироваться)...

У вас выписана формула разложения по базису масштабированных синков. В этом базисе можно разложить любую функцию, частотный спектр которой лежит в $\Omega = [-\pi/\Delta, \pi/\Delta]$. Частотный спектр функции $sinc(x) = \sin\pi x/\pi x$ есть полоса $[-\pi, \pi]$, то есть лежит в $\Omega$, если $\Delta\leqslant 1$...

ПС: лично я интерпретировал бы это задания, как выбор хорошего $\Delta$. Ну, в смысле такого, при котором ряд получился бы простым. Например, при $\Delta = 1$ имеем $sinc(x) = sinc(x)$, неинтересно. Интересным может оказаться $\Delta= 1/2$... в таком духе, в общем;)

[ivantozavrus]

ihq.pl в сообщении #1616434 писал(а):

ivantozavrus

У вас выписана формула разложения по базису масштабированных синков. В этом базисе можно разложить любую функцию, частотный спектр которой лежит в $\Omega = [-\pi/\Delta, \pi/\Delta]$. Частотный спектр функции $sinc(x) = \sin\pi x/\pi x$ есть полоса $[-\pi, \pi]$, то есть лежит в $\Omega$, если $\Delta\leqslant 1$...

ПС: лично я интерпретировал бы это задания, как выбор хорошего $\Delta$. Ну, в смысле такого, при котором ряд получился бы простым. Например, при $\Delta = 1$ имеем $sinc(x) = sinc(x)$, неинтересно. Интересным может оказаться $\Delta= 1/2$... в таком духе, в общем;)

Только вот у меня получилось, что частотная полоса $sinc(\pi x) -$ это $[-\pi/2, \pi/2].$
И все же, получается что $\Delta = \frac{1}{2 \pi}$ так же под условие принадлежности спектра $\Omega = [-\pi/\Delta, \pi/\Delta]$ попадает. А разложения хорошего не получается. Где здесь подводный камень?

[svv]

Вероятно, ihq.pl намекает, что здесь у Вас ошибка:

ivantozavrus в сообщении #1616305 писал(а):

f(i) = f(i) + f_k(c) * sinc(pi/delta * (x(i) - x_k(c)));

Исправьте, и всё будет хорошо.

Comment. Существует два стандарта в определении функции $\operatorname{sinc}$. В MATLAB используется первый. Вы записали Ваш ряд через $\operatorname{sinc}$ по второму (или просто ошиблись).

Какой стандарт имел в виду автор задания в выражении $\operatorname{sinc}(\pi x)$ — загадка. С одной стороны, если подразумевалось использование MATLAB, то, наверное, первый. Да и, как сказано по ссылке, в ЦОС обычно используется первый. С другой — вид выражения как бы намекает, что $\pi$ не входит в определение $\operatorname{sinc}$ (=второй стандарт), поэтому его приходится вводить дополнительным множителем в аргументе функции. Хотя, может быть, подразумевалась именно функция $\frac{\sin\pi^2 x}{\pi^2 x}$, так что это гадание.

[ivantozavrus]

svv в сообщении #1616452 писал(а):

Вероятно, ihq.pl намекает, что здесь у Вас ошибка:

ivantozavrus в сообщении #1616305 писал(а):

f(i) = f(i) + f_k(c) * sinc(pi/delta * (x(i) - x_k(c)));

Исправьте, и всё будет хорошо.

Comment. Существует два стандарта в определении функции $\operatorname{sinc}$. В MATLAB используется первый. Вы записали Ваш ряд через $\operatorname{sinc}$ по второму (или просто ошиблись).

Какой стандарт имел в виду автор задания в выражении $\operatorname{sinc}(\pi x)$ — загадка. С одной стороны, если подразумевалось использование MATLAB, то, наверное, первый. Да и, как сказано по ссылке, в ЦОС обычно используется первый. С другой — вид выражения как бы намекает, что $\pi$ не входит в определение $\operatorname{sinc}$ (=второй стандарт), поэтому его приходится вводить дополнительным множителем в аргументе функции. Хотя, может быть, подразумевалась именно функция $\frac{\sin\pi^2 x}{\pi^2 x}$, так что это гадание.

Спасибо огромное! Теперь все получилось. Подразумевался все-таки $sinc(x)$ в первом смысле.
Получается, что я ошибся в опредлении $\Delta = \frac{1}{2 F},$ и под $F$ нужно понимать полуширину частотной полосы $f?$

[svv]

Ошибка в этом преобразовании (правая часть была закодирована в MATLAB):

$$\sum_{k=-\infty}^{\infty} f(k \Delta) \frac{\sin \biggr(\frac{\pi}{\Delta} (x - k \Delta)\biggr)}{\frac{\pi}{\Delta} (x - k \Delta)}=\sum_{k=-\infty}^{\infty} f(k \Delta) \operatorname{sinc}\left({\frac{\color{magenta}\pi}{\Delta} (x - k \Delta)}\right)$$

А именно, вместо ${\color{magenta}\pi}$ надо поставить единичку, так как по первому стандарту $\pi$ уже входит в определение $ \operatorname{sinc}$.

[svv]

Преобразование Фурье от $\operatorname{sinc}(ax)$ можно найти в справочнике, например, здесь (таблица, строка 202). При $a>0$ оно равно $\frac 1{a}\operatorname{rect}\frac {\xi}{a}$, где $\xi$ — обычная (не циклическая) частота, а $\operatorname{rect}$ — прямоугольная функция:

$$\operatorname{rect}\frac{\xi}{a} = \begin{cases} 0, & \text{if } |\xi| > \frac{a}{2} \\ \frac{1}{2}, & \text{if } |\xi| = \frac{a}{2} \\ 1, & \text{if } |\xi| < \frac{a}{2}. \end{cases}\right$$

Поскольку у нас $a=\pi$, верхняя граничная частота спектра $f_c=\frac a 2=\frac{\pi}2$. Тогда интервал дискретизации $\Delta\leqslant \frac{1}{2f_c}=\frac{1}{\pi}$.

Преподаватель хотел, чтобы Вы взяли различные значения delta и убедились, что
delta = 1/pi*0.5; % работает
delta = 1/pi*0.8; % работает
delta = 1/pi*1.0; % работает
delta = 1/pi*1.2; % не работает

[ivantozavrus]

svv в сообщении #1616468 писал(а):

Ошибка в этом преобразовании (правая часть была закодирована в MATLAB):

$$\sum_{k=-\infty}^{\infty} f(k \Delta) \frac{\sin \biggr(\frac{\pi}{\Delta} (x - k \Delta)\biggr)}{\frac{\pi}{\Delta} (x - k \Delta)}=\sum_{k=-\infty}^{\infty} f(k \Delta) \operatorname{sinc}\left({\frac{\color{magenta}\pi}{\Delta} (x - k \Delta)}\right)$$

А именно, вместо ${\color{magenta}\pi}$ надо поставить единичку, так как по первому стандарту $\pi$ уже входит в определение $ \operatorname{sinc}$.

Я имею ввиду, что полоса частот $f \in [-0.5, 0.5].$ Тогда $F = 1,$ а $\Delta \leq \frac{1}{2},$ если брать за $F$ всю полосу. Но точно восстановить функцию получается уже при $\Delta = 1.$ Значит, надо за $F$ принимать полуполосу.

[svv]

ivantozavrus в сообщении #1616479 писал(а):

Я имею ввиду, что полоса частот $f \in [-0.5, 0.5].$ Тогда $F = 1,$ а $\Delta \leq \frac{1}{2},$ если брать за $F$ всю полосу. Но точно восстановить функцию получается уже при $\Delta = 1.$ Значит, надо за $F$ принимать полуполосу.

Куда у Вас делись все $\pi$? Правильное максимальное значение интервала дискретизации $\Delta=\frac 1{\pi}$. Мне не хотелось бы использовать величину $F$, потому что она уже "скомпрометирована" (мы подбираем для неё определение, исходя из согласия с численным экспериментом, и уже есть два разных варианта). Пожалуйста, попробуйте разобраться в моих обозначениях (они такие же, как в Википедии, статья Теорема Котельникова):

ivantozavrus в сообщении #1616305 писал(а):

По графику определил, что спектр становится нулевым после $f = \frac{\pi}{2}$

Это правильно. Итак, спектр сигнала ненулевой только в полосе частот $f\in[-f_c;f_c]$, где в нашем случае $f_c= \frac{\pi}{2}$.
Тогда, по т. Котельникова, сигнал восстанавливается при
$0<\Delta \leqslant \frac {1}{2f_c}$
Подставляя сюда $f_c$, находим
$0<\Delta \leqslant \frac {1}{\pi}$
И это подтверждается численными экспериментами. При $\Delta \leqslant \frac {1}{\pi}$ сигнал восстанавливается, при $\Delta > \frac {1}{\pi}$ нет.

Кстати, посмотрите мой вариант скрипта, может, какие-то приёмчики понравятся.

Используется синтаксис Matlab M

x = -5:0.01:5;
deltamax = 1/pi; % максимальный интервал дискретизации по т.Котельникова
delta = deltamax*0.8; % попробуйте другие значения коэффициента
kmax = 100; % берем отсчёты с номерами -kmax:kmax
signal = @(x) sinc(pi*x); % anonymous function

y = zeros(size(x)); % аппроксимация
for k = -kmax:kmax
    x_k = delta*k;
    f_k = signal(x_k);
    y = y + f_k * sinc((x-x_k)/delta);
end

hold on
plot(x, signal(x));
plot(x, y);
hold off