Помощь - Поиск - Пользователи - Календарь
Полная версия этой страницы: Окна для приближения производных sinc
Форум разработчиков электроники ELECTRONIX.ru > Cистемный уровень проектирования > Математика и Физика
_Ivana
При некоторых методах интерполяции для получения значения производят сложение близлежащих точек с определенными весовыми коэффициентами, задаваемыми таблицей оконной функции (если положение точки внутри отрезка строго определено, например 1/2 и т.п.). В качестве этих коэффициентов можно брать:

1) Значения базисных полиномов Лагранжа, рассчитанных по тому же количеству точек что содержатся в окне - тогда мы имеем "чистую" полиномиальную Лагранжевскую интерполяцию, которая становится точнее с ростом частоты дискретизации но дает большие ошибки вблизи частоты Найквиста
2) Значения функции sinc (конечное их число). Точность определяется размером окна и по моему имху равна 1/N, где N-размер окна. Но зато хорошо себя ведет вблизи частоты Найквиста.
3) Другие хитрые окна, которые наверное можно представить как произведение sinc на квазипрямоугольное окно с хитрыми загибами по краям, что имхо должно приводить к некоему компромиссу характеристик 1) и 2)

Далее: мне хочется рассчитывать производные любых порядков в любы точках (хотя для начала в самих точках дискретизации) для равномерно дискретизированной функции. Появляются возможные аналогии описанному выше:

1) Окно из производных базисных полиномов Лагранжа
2) Окно из конечного набора усеченных производных sinc (cos(x)/x - sin(x)/x^2 и т.п.)
3) Некие компромиссные окна из 1) и 2)

1) и 2) я рассчитаю самостоятельно, собственно, вопрос:
как я могу получить оптимальное окно типа 3) для производных? Может есть какое-то известное квазипрямоугольное окно, наложив которое на 2) я получу желаемое? Или если такое окно есть и оно хорошо работает для самой функции, то не факт что будет хорошо работать для производных?

Сорри за многабукаф и заранее спасибо за ответы. Если вдруг что непонятно в моем изложении - готов уточнить.
AndrewN
QUOTE (_Ivana @ Jul 7 2012, 16:01) *
Посмотрите в учебнике Бахвалова, Жидкова, Кобелькова Численные методы главы по интерполяции, аппроксимации и конечным разностям.

_Ivana
Скачал книжку, читаю, но пока не нашел ответа на свой вопрос. Приведу конкретный пример: есть ряд значений y(i), надо аппроксимировать первую производную в точке 0. Если ограничиваться аппроксимацией по 9 точкам, получаем формулу:

y'(0) = a1*(y(1)-y(-1)) + a2*(y(2)-y(-2)) + a3*(y(3)-y(-3)) + a4*(y(4)-y(-4))

1) по Лагранжу:
a1 = 4/5
a2 = -1/5
a3 = 4/105
a4 = -1/280

2) по обрезанной производной синка (если ничего не напутал, смущает знак):
a1 = -sinc'(1)
a2 = -sinc'(2)
a3 = -sinc'(3)
a4 = -sinc'(4)

Причем, очевидно, что проверку прямой линией: 2*a1 + 4*a2 + 6*a3 + 8*a4 = 1 выдерживает только Лагранж, а синк будет сильно врать. С увеличением количества точек синк будет врать все меньше (~1/N).

Вопрос - какие значения a1~a4 вы предложите? И почему такие? Или (что то же самое) - какое квазипрямоугольное окно для этих коэффициентов лучше всего наложить на значения по синку, и по каким критериям выбирать это окно?
Tiro
Цитата(_Ivana @ Jul 7 2012, 21:02) *
Скачал книжку, читаю, но пока не нашел ответа на свой вопрос. Приведу конкретный пример: есть ряд значений y(i), надо аппроксимировать первую производную в точке 0. Если ограничиваться аппроксимацией по 9 точкам, получаем формулу:

y'(0) = a1*(y(1)-y(-1)) + a2*(y(2)-y(-2)) + a3*(y(3)-y(-3)) + a4*(y(4)-y(-4))

1) по Лагранжу:
a1 = 4/5
a2 = -1/5
a3 = 4/105
a4 = -1/280

2) по обрезанной производной синка (если ничего не напутал, смущает знак):
a1 = -sinc'(1)
a2 = -sinc'(2)
a3 = -sinc'(3)
a4 = -sinc'(4)

Причем, очевидно, что проверку прямой линией: 2*a1 + 4*a2 + 6*a3 + 8*a4 = 1 выдерживает только Лагранж, а синк будет сильно врать. С увеличением количества точек синк будет врать все меньше (~1/N).

Вопрос - какие значения a1~a4 вы предложите? И почему такие? Или (что то же самое) - какое квазипрямоугольное окно для этих коэффициентов лучше всего наложить на значения по синку, и по каким критериям выбирать это окно?

Тут важно, что Вы хотите получить в ответе. В книжках обычно стараются найти компромисс между порядком приближения и расширением спектра. Или даже отношение сигнал/шумоптимизируют. Если есть заранее известные характеристики сигнала, то надо на них опираться.
AndrewN
QUOTE (_Ivana @ Jul 7 2012, 22:02) *
С вычислительной точки зрения совершенно все равно какую модель дифференцировать - например, если функция аппроксимирована кубическим сплайном, то её производная - полином второй степени и считается аналитически.

Если функция задана таблицей, то считаются конечные разности, с фильтрацией высокочастотного шума, который вносит наибольший вклад в ошибку. Т.о. необходим копромисс - поскольку производная это локальное свойство, то чем меньше точек, тем лучше (но меньше чем двумя всё равно не обойтись), а для подавления шума точек желательно взять больше.

В обоих случаях считается не сама производная, а её аппроксимация и она можеть быть близка или далека от истинной производной.

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




QUOTE (Tiro @ Jul 7 2012, 23:49) *
компромисс между порядком приближения и расширением спектра
Для дифференцирования спектр нужно сузить. Но не перестараться. Если (как предел фильтра нижних частот) взять среднее арифметическое, то получим много раз гладкую функцию с производными равными нулю и очень узким спектром...

Tiro
Цитата(AndrewN @ Jul 7 2012, 23:15) *
предел фильтра нижних частот) взять среднее арифметическое, то получим много раз гладкую функцию с производными равными нулю и очень узким спектром...

Извините, но мне непонятно что такое "предел фильтра нижних частот".
_Ivana
Спасибо за направление, я все о своем, сугубо теоретическом.
Обнаружил (если нигде не ошибся), что значение производной sinc в точках натурального ряда N с меняющимся знаком совпадают со значением 1/N, например sinc'(3) = -1/3 и т.д. То есть, если взять прямую линию (к чему стремится любая функция при увеличении частоты дискретизации), то ряд получения первой производной даже не является сходящимся, то есть любое большое число окрестных точек не даст нам аппроксимации производной...
Tiro
Цитата(_Ivana @ Jul 7 2012, 23:31) *
Спасибо за направление, я все о своем, сугубо теоретическом.
Обнаружил (если нигде не ошибся), что значение производной sinc в точках натурального ряда N с меняющимся знаком совпадают со значением 1/N, например sinc'(3) = -1/3 и т.д. То есть, если взять прямую линию (к чему стремится любая функция при увеличении частоты дискретизации), то ряд получения первой производной даже не является сходящимся, то есть любое большое число окрестных точек не даст нам аппроксимации производной...

Функция это абстракция. Она вообще ни к чему не стремится. И sinc тоже абстракция. Это выражение такое, ведет себя как полагается в области определения )))
Что хочешь понять-то?
_Ivana
Я хочу понять, есть ли такое окно (и найти его), которое при относительно малом количество точек (~10-20) достаточно хорошо приближает первую (другое окно - вторую и т.д.) производную в точке. В случае полиномов Лагранжа любой степени я эту задачу решил. Теперь хочу экстраполировать на случай синка или чего-то промежуточного.

ЗЫ синк на прямой y(i) = 1 сходится к этой прямой. А первая производная нет, если я не ошибаюсь. Сейчас буду пробовать вторую, там навскидку график обнадеживает... График например есть вот здесь в самом низу
http://www.tvalx.com/Russian/MathArticles/...ncFunction.html
AndrewN
QUOTE (Tiro @ Jul 7 2012, 23:24) *
Извините, но мне непонятно что такое "предел фильтра нижних частот"
Среднее арифметическое а. Спектр а равен a*sqrt(2*pi)*delta(omega). Уже не бывает.



QUOTE (_Ivana @ Jul 7 2012, 23:46) *
Я хочу понять, есть ли такое окно (и найти его), которое при относительно малом количество точек (~10-20) достаточно хорошо приближает первую (другое окно - вторую и т.д.) производную в точке
Производная - локальная величина. Окна для вычисления производной не нужны. И sinc(*) вам абсолютно не нужен. Вы считаете производную в точке. Закройте глаза и представьте, что функция продолжается от неё и влево и вправо так далеко, что синк от дельты не отличить.



QUOTE (Tiro @ Jul 7 2012, 23:40) *
Функция это абстракция. Она вообще ни к чему не стремится
Ура, завтра в школу не пойдем - пределы отменили!
_Ivana
Цитата
Производная - локальная величина. Окна для вычисления производной не нужны. И sinc(*) вам абсолютно не нужен. Вы считаете производную в точке.

Правильно. Но значение функции в точке - тоже локальная величина. Однако при интерполяции её считают и через окно. Суть в том, что если мы представляем аппроксимирующий полином в форме базисных полиномов Лагранжа, то и производную можем рассчитать так же (я проверял, работает). А если при бесконечном увеличении точек базисные полиномы Лагранжа вырождаются в синк, то и там думал что это возможно.

ЗЫ во всяком случае, вторая производная синка в точках натурального ряда спадает быстрее чем первая, что обнадеживает (пусть хотя бы вторую удастся аппроксимировать).
Tiro
Цитата(AndrewN @ Jul 7 2012, 23:56) *
Среднее арифметическое а. Спектр а равен a*sqrt(2*pi)*delta(omega). Уже не бывает.

То есть если а 0, +1, 0, -1 и далее, а delta(omega) 0, +1, 0 , -1 то быстро через тернии к звездам?

Цитата(AndrewN @ Jul 7 2012, 23:56) *
Ура, завтра в школу не пойдем - пределы отменили!

Да, не ходите. Завтра -30 за порогом!
P.S. А как оттает, и про пределы поговорим.
AndrewN
QUOTE (Tiro @ Jul 8 2012, 01:04) *
То есть если а 0, +1, 0, -1 и далее
Нет. а это и есть среднее арифметическое. Была таблица значений функции (возможно что с шумом) на равномерной сетке x(i), i = 1,...,N, а получили a = 1/N*SUM(k=1,...,N)x(k).

Tiro
Цитата(_Ivana @ Jul 7 2012, 23:46) *
Я хочу понять, есть ли такое окно (и найти его), которое при относительно малом количество точек (~10-20) достаточно хорошо приближает первую (другое окно - вторую и т.д.) производную в точке. В случае полиномов Лагранжа любой степени я эту задачу решил. Теперь хочу экстраполировать на случай синка или чего-то промежуточного.

ЗЫ синк на прямой y(i) = 1 сходится к этой прямой. А первая производная нет, если я не ошибаюсь. Сейчас буду пробовать вторую, там навскидку график обнадеживает... График например есть вот здесь в самом низу
http://www.tvalx.com/Russian/MathArticles/...ncFunction.html

В общем случае нет такого окна. Да и сами определения производных и проч сделаны для ограниченного набора функций. Если Ваша функция непредставима в имеющихся базисах, то, возможно, Вы найдете свое "окно". И оно будет представимо в имеющихся базисах. biggrin.gif
_Ivana
Предполагается что исходная функция дифференцируема нужное количество раз и удовлетворяет условию теоремы Котельникова. Однако, вместо урезанного синка используют же другие окна, которые были получены (из синка?) по каким-то критериям. Вот я и хочу найти такие варианты для производных.
Tiro
Цитата(_Ivana @ Jul 8 2012, 00:32) *
Предполагается что исходная функция дифференцируема нужное количество раз и удовлетворяет условию теоремы Котельникова. Однако, вместо урезанного синка используют же другие окна, которые были получены (из синка?) по каким-то критериям. Вот я и хочу найти такие варианты для производных.

Надеюсь, про связь частотно-временных представлений Вы в курсе. Теперь вопрос - зачем появились окна? Все просто. Отрезок на оси времени дает бесконечный частотный спекр. Отрезок на оси частот - бесконечен во времени. Ну такая у нас математика. Все хотят ограничить сигнал как во времени, так и по частоте.
Окна - это результат компромисса. Они имеют как частотные, так и временные представления. Для детерминированного сигнала всегда можно придумать согласованный фильтр, но его всегда можно описать в имеющихся ортогональных базисах.
_Ivana
Цитата
Окна - это результат компромисса.

Правильно. И бесконечный синк физически нереализуем, поэтому делают конечное его окно. Можно просто "обрубать" и получать краевые эффекты, а можно отклониться от чистого синка ближе к краям окна, получив другое окно.
Но я пока пробую интерполировать производную сигнала хотя бы простым ограничением производной синка.

ЗЫ сигнал в базисах можно описать, можно много чего ещё. Но на текущий локальный момент мне интересно попробовать увеличить точность аппроксимации производной.
Tiro
Цитата(_Ivana @ Jul 8 2012, 00:52) *
Правильно. И бесконечный синк физически нереализуем, поэтому делают конечное его окно. Можно просто "обрубать" и получать краевые эффекты, а можно отклониться от чистого синка ближе к краям окна, получив другое окно.
Но я пока пробую интерполировать производную сигнала хотя бы простым ограничением производной синка.

ЗЫ сигнал в базисах можно описать, можно много чего ещё. Но на текущий локальный момент мне интересно попробовать увеличить точность аппроксимации производной.

Интерполяции нужны только в случае, если невозможно точно вычислить сигнал. Все остальное - лирика.
_Ivana
За обилием лирики потерялась конкретика. Осмелюсь повториться:
Цитата
Приведу конкретный пример: есть ряд значений y(i), надо аппроксимировать первую производную в точке 0. Если ограничиваться аппроксимацией по 9 точкам, получаем формулу:

y'(0) = a1*(y(1)-y(-1)) + a2*(y(2)-y(-2)) + a3*(y(3)-y(-3)) + a4*(y(4)-y(-4))

1) по Лагранжу:
a1 = 4/5
a2 = -1/5
a3 = 4/105
a4 = -1/280

2) по обрезанной производной синка (если ничего не напутал, смущает знак):
a1 = -sinc'(1)
a2 = -sinc'(2)
a3 = -sinc'(3)
a4 = -sinc'(4)

Причем, очевидно, что проверку прямой линией: 2*a1 + 4*a2 + 6*a3 + 8*a4 = 1 выдерживает только Лагранж, а синк будет сильно врать. С увеличением количества точек синк будет врать все меньше (~1/N).

Вопрос - какие значения a1~a4 вы предложите? И почему такие?
iiv
на б-сплайнах у Вас получится
Код
a1=1
a2=-(2-sqrt(3))
a3=(2-sqrt(3))^2
a4=-(2-sqrt(3))^3
_Ivana
С последним постом параметр конкретика/лирика в этом топике скачком возрос от нуля до своего максимального значения, производная соответственно максимальна sm.gif
iiv, спасибо, попробую ваши коэффициенты на 9 точках. Но один вопрос - если взять 11, 13,.... 101 и т.д. точек - то коэффициенты будут рассчитываться по аналогии как
Код
an=+-(2-sqrt(3))^n

или не все так просто и выражение в скобках тоже будет другое? Спрашиваю потому что пока не догадываюсь как Вы получили эти коэффициенты.

Покрутил ещё Лагранжа. Обнаружил, что его коэффициенты представляют собой (для первой производной) значения производной синка в точках натурального ряда (они равны 1/n) умноженные на некие факториальные коэффициенты, которые вначале получаются около 1-цы но к концу нашего конечного окна = степени полинома резко спадают до 0, причем настолько резко, что даже логарифм их значения улетает в минус бесконечность почти вертикально. Поэтому и получается, что если мы будем брать предельный случай - только производные синка, то длина окна получается бесконечна и эти факториальные коэффициенты не ограничат чистые производные синка = 1/n и ряд не будет сходящимся.
В качестве эксперимента взял окно в 50 точек (по полиному Лагранжа 100 степени, да sm.gif ). После ~30 значений точности моего 1С-симулятора (15 знаков после запятой в десятичном формате) не хватает для этих коэффициентов - они получаются нулевые, хотя теоретически они есть все вплоть до 50-го. Если рассчитать окно по 80 точкам, то как раз 50 ненулевых коэффициентов и получится при моей точности их определения. А если рассчитать окно по 200 коэффициентам (полиному 400-й степени sm.gif ) и взять из них первые 50, то точность аппроксимации производной ограничится, после определенного порядка не будет увеличиваться при увеличении частоты дискретизации (за счет отбрасывания коэффициентов после определенного порядка), но увеличится точность вблизи частоты Найквиста.
Таким образом, имея определенную длину окна (наши 50 коэффициентов) я могу в неких пределах варьировать соотношение между абсолютной точностью аппроксимации при больших и малых частотах дискретизации, получать некие компромиссы. Но я манипулирую лишь порядком полинома и количеством оставляемых из него значений для окна. Единственно что ещё - последнее значение окна я могу рассчитать не как значение по полиному, а по правилу схождения ряда к 1-це на прямой. При этом абсолютная точность при больших частотах дискретизации несколько повышается. Но у меня есть предположения что можно обрезать ряд коэффициентов с изменением значений не одного а многих последних элементов, так чтобы погрешность была ещё меньше а окно оптимальнее. Пока не могу придумать как это сделать, буду признателен если кто подскажет.
iiv
Цитата(_Ivana @ Jul 11 2012, 14:48) *
С последним постом параметр конкретика/лирика в этом топике скачком возрос от нуля до своего максимального значения, производная соответственно максимальна sm.gif
iiv, спасибо, попробую ваши коэффициенты на 9 точках. Но один вопрос - если взять 11, 13,.... 101 и т.д. точек - то коэффициенты будут рассчитываться по аналогии как
Код
an=+-(2-sqrt(3))^n

или не все так просто и выражение в скобках тоже будет другое? Спрашиваю потому что пока не догадываюсь как Вы получили эти коэффициенты.


да, именно так. Это обычный Б-сплайн 3-ей степени. Для пятой я на раз форулы из головы написать не смогу, старый стал. А ведь было время, на экзаменах студентам такое давал, а они выводили, и приходилось на лету проверять.

Еще есть возможность нарисовать тот же Б-сплайн, но на сетке в несколько раз большей, чем снятые Ваши эксперименты. Попробуйте, пожалуйста, сами, если сложно будет, пишите, я формулы для оного перевыведу и здесь запостю.... Я когда-то постил, кажется, для des00, но уже забыл, ссылку не могу найти, и наизусть не помню, а перевыводить лениво.

С Лагранжем и синком, можно Вам один глупый совет дать, и то, и другое хорошо осциллирует. Не правильно это при обычной нашей природе. Или сплайны, или Чебышев, но не Лагранж и синком.
_Ivana
Цитата
Еще есть возможность нарисовать тот же Б-сплайн, но на сетке в несколько раз большей, чем снятые Ваши эксперименты. Попробуйте, пожалуйста, сами, если сложно будет, пишите, я формулы для оного перевыведу и здесь запостю....

Вы имеете в виду Б-сплайн построенный по.... а кстати скольким точкам? В несколько раз больше - это касается моего примера для 9 точек или для последнего для 101 точки? В любом случае вопросы сплайнов мне интересны, я попробую то что Вы говорите, только для начала хочу понять что именно имеется в виду.

Цитата
С Лагранжем и синком, можно Вам один глупый совет дать, и то, и другое хорошо осциллирует. Не правильно это при обычной нашей природе. Или сплайны, или Чебышев, но не Лагранж и синком.

1) Лагранж осциллирует в краях диапазона, например если взять точек 50, построить полином по всем точкам - в краях будут осцилляции. А если взять те же 50 точек, построить полином, но использовать его только для единственного центрального интервала, а для соседнего использовать свой полином, то имхо осцилляций не должно быть.
2) Собственно, сам вопрос аппроксимации производных у меня возник вследствие построения различных сплайнов с их использованием, в т.ч. и обобщенных сплайнов Эрмита. Например, сейчас мучаю локальный сплайн Эрмита 9-й степени с непрерывными 0,1,2,3,4 производными, при точном значении производных точность сплайна получается имхо весьма хорошая, все дело за точной аппроксимацией производных вплоть до 4-й в узлах сетки. А вот эти производные я и аппроксимирую Лагранжем, ибо больше пока не знаю как. Но, например, при моем окне в 50 точек производные аппроксимируются с более чем достаточной точностью для того же сплайна Эрмита 5-й и 7-й степени. Если интересно, могу показать экспериментальные графики точности этих сплайнов - как теоретически предельные, так и практически достигнутые sm.gif
_Ivana
Графики точности сплайна с аппроксимацией производных по ограниченному окну из 50 точек. Красный - теоретический максимум - сплайн с известными производными, зеленый - с аппроксимированными. На загиб графиков после 16 точек на период не обращайте внимания - просто у меня синус считается с такой точностью, на самом деле прямая линия продолжается вниз с тем же наклоном. Слева окно по полиному 200 порядка, справа - 400. На графике слева точность сплайна на всем участке справа не отклоняется от максимума но немного хуже слева - около частоты Найквиста. Справа точность ограничена на уровне Е-8, но лучше при малых частотах дискретизации. Степень ограничения точности сплайна можно варьировать в рамках окна из наших 50 точек путем выбора порядка полинома и, как следствие, коэффициентов окна.
_Ivana
iiv, попробовал Ваши коэффициенты.
Слева график точности кубических Эрмитовых сплайнов: синий с точными производными, красный с аппроксимацией производных по Лагранжу 8 порядка (те самые 4 коэффициента из моего примера), зеленый - аппроксимация производных 4-мя коэффициентами по Б-сплайну 3 степени (Ваши коэффициенты).
Справа график точности сплайнов Эрмита 9 степени: красный с точными производными, зеленый - с аппроксимацией через окно из 50 точек, построенное по вашим коэффициентам с дальнейшим расширением степенного ряда. Кстати по этому ряду 39-й коэффициент окна уже равен 1Е-20 а остальные не влезли в точность и равны нулю. Видно что точность аппроксимации производной не растет с увеличением длины окна. Думаю это связано с тем, что эти выкладки для Б-сплайна 3 степени и мы не можем вылезти за рамки точности этой 3 степени как бы ни расширяли окно. Повторюсь, свои окна по Лагранжу, которые работают хорошо, я считаю по полиномам 200-400 степеней.
Для просмотра полной версии этой страницы, пожалуйста, пройдите по ссылке.
Invision Power Board © 2001-2025 Invision Power Services, Inc.