Апроксимация сплайнами, оптимизация размещения контрольных точек Опции

[ataradov]

Эта тема является логическим продолжением вот этой, но может быть интересна и сама по себе, поэтому выношу в отдельный типик.

На картинках изображена характерная кривая (1 - общий план, 2 - приближенное начало), которую нужно представить сплайнами или любым другим способом, занчительно снижающим объем памяти, требуемый для ее хранения, но при этом допускающий сравнительно простое восстановление по точкам, не разжимя всей кривой сразу (то-есть запаковать ее ZIP-ом нельзя ).

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

В идеале нужен алгоритм, который позволит по заданному числу контрольных точек определть их положение так, чтобы СКО было минимальным. Или наоборот, по заданному СКО определить требуемое число точек и их положение.

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

Эскизы прикрепленных изображений

 

[Tanya]

На картинках изображена характерная кривая (1 - общий план, 2 - приближенное начало), которую нужно представить сплайнами или любым другим способом, занчительно снижающим объем памяти, требуемый для ее хранения, но при этом допускающий сравнительно простое восстановление по точкам, не разжимя всей кривой сразу (то-есть запаковать ее ZIP-ом нельзя ).

Мне кажется, что это композиция экспоненциального (?) спада + (-) довольно равномерно идущие затухающие пики.

[ataradov]

Мне кажется, что это композиция экспоненциального (?) спада + (-) довольно равномерно идущие затухающие пики.

Нет, сигнал сложнее, кроме того это только одна из кривых (это амплитуда перавой гармоники ноты C4), остальные хотя и похожи по форме (всегда резко растут и плавно спадают), но немного отличаются и различия эти критичны.

[Andrey_1]

Мне кажется, что это композиция экспоненциального (?) спада + (-) довольно равномерно идущие затухающие пики.

Поставьте Origin Calc и грузите туда ASCII файлы ваших функций
В проге есть нормальный сплан генератор который Вам отапроусимирует хоть средними квадратами хоть Савицким-Голеем
и еще уравнение полинома на экране покажет чем велосипед изобретать
ПС Шкалы опять в йагупопах 87-ых

Сообщение отредактировал Andrey_1 - Jan 23 2011, 14:21

[ataradov]

Поставьте Origin Calc и грузите туда ASCII файлы ваших функций

Цель - узнать алгоритм, а не пользоваться программой. Мне таких графиков несколько сотен нужно будет обработать, все руками предложите?

[Andrey_1]

Цель - узнать алгоритм, а не пользоваться программой. Мне таких графиков несколько сотен нужно будет обработать, все руками предложите?

Вам таки ехать или шашечки?
Узнаете уравнение огибающей - напишите в С или Матлабище как считать
Файлы грузить - скрипт нарисуйте
Matlab command window:
help polyfit
help polyfun

[ataradov]

Вам таки ехать или шашечки?

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

[Andrey_1]

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

Я вам предлагаю варианты решения
Выпад(к)ов и в помине нет - это юмор сэр

[АНТОН КОЗЛОВ]

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

[Guest_TSerg_*]

Для Вашего случая моделировать надо не результат, а причину.

[_Pasha]

Обратите внимание, что регрессионный полином не обязан в точности совпадать ни с одной из точек входного массива данных, в то время как ss2tf или как его там по-матлабовски... не помню - в общем, обязан совпадать в узлах интерполяции.
В итоге, Вам надо нечто гибридное - солвер для определения оптимального сегмента сплайн-аппроксимации. Внутри сегмента - обеспечение мин. отклонения, "снаружи" - минимизация кол-ва коэффициентов. Безусловным приоритетом обладает максимизация длительности сегмента.

[scifi]

В порядке бреда: коль скоро это звуки (если я не ошибаюсь), то отчего бы не применить алгоритмы из lossy codecs (MP3 и прочие)?

[Tanya]

Для Вашего случая моделировать надо не результат, а причину.

Наверное, не только в этом случае. Если верить авторскому рисунку, кроме затухания происходит перенос энергии основной моды в гармоники. И обратно?

[GetSmart]

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

Как определять точки - сильно зависит от алгоритма интерполяции/сглаживания. Для начала определитесь с методом сглаживания. А потом можно думать об алгоритме поиска точек.

[_Pasha]

Для начала определитесь с методом сглаживания.

Можно на автопилоте, например, повышая порядок полинома Лежандра для осциллирующих участков.

[ataradov]

Для Вашего случая моделировать надо не результат, а причину.

Согласен. Это первый шаг, так как моделирование причины - это не так просто и даже для самых сложных и совершенных доступных публично моделей качество синтезируемого звука оставляет желать лучшего.

Да и почти все статьи рассматривающие такой синтез строят именно огибающие гармоник решая системы ДУ, а сам синтез делают по гармоникам, иначе никакого реал-тайма не получается.

Как определять точки - сильно зависит от алгоритма интерполяции/сглаживания. Для начала определитесь с методом сглаживания. А потом можно думать об алгоритме поиска точек.

Согласен. Рассматриваются спростые кубические сплайны и эрмитовы сплайны. Оба типа обладают интересными характеристиками и с каким получится лучше пока не ясно.

Для начала можно взять классические кубические сплайны.

Наверное, не только в этом случае. Если верить авторскому рисунку, кроме затухания происходит перенос энергии основной моды в гармоники. И обратно?

Происходит. Именно поэтому и записываются из сэмплов все гармоники до 15 кГц примерно. А при синтезе обратно собираются.

[petrov]

Почему бы просто не использовать более высокую частоту дискретизации на участке с более широкой полосой?

[ataradov]

Почему бы просто не использовать более высокую частоту дискретизации на участке с более широкой полосой?

Для хранения гармоник или сэмпла? Предлагается хранить просте передискретизированный сигнал или я не понял?

[Tanya]

Происходит. Именно поэтому и записываются из сэмплов все гармоники до 15 кГц примерно. А при синтезе обратно собираются.

Вот и попробуйте написать модель переноса... И посчитать...
Примерно так... Имеется вероятность перехода и туда и обратно. Или можно в терминах скорости. В результате образуется равновесие. При падении энергии (амплитуды (основной моды)?) ниже определенного порога каналы обмена закрываются.

[petrov]

Для хранения гармоник или сэмпла? Предлагается хранить просте передискретизированный сигнал или я не понял?

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

[ataradov]

Вот и попробуйте написать модель переноса... И посчитать...
Примерно так... Имеется вероятность перехода и туда и обратно. Или можно в терминах скорости. В результате образуется равновесие. При падении энергии (амплитуды (основной моды)?) ниже определенного порога каналы обмена закрываются.

Да, я об этом думал, но решил оставить на потом, после сжатия. Кроме того началный момент атаки скорее всего так описать не получится.

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

Тогда даже самый неоптимальный алгоритм расстановки ключевых точек в стлайнах дает большее сжатие.

[petrov]

Тогда даже самый неоптимальный алгоритм расстановки ключевых точек в стлайнах дает большее сжатие.

А где обоснование использования сплайнов? Может во время контакта с молотчком гармоники шумом модулированы и он ни в какие сплайны не вписывается? Ведь музыкальные файлы просто сплайнами не сжимают несмотря на "большее сжатие"?

[ataradov]

А где обоснование использования сплайнов?

Его нет. Я и хотел представить и посмотреть (послушать). Вероятно действительно, не стоит сразу гнаться за самым оптимальным сжатием и представить сначала сплайнами на регулярной сетке. Я просто так набросился, потму что задача казалась простой, а оказалась не очень.

Касательно сжатия музыки у меня есть только один аргумент - у огибающих очень низкая частота.

Сообщение отредактировал Taradov Alexander - Jan 24 2011, 12:43

[petrov]

у огибающих очень низкая частота.

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

[ataradov]

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

Сплайнами можно добиться побитового совпадения результата и дискретизированного до 16 бит сигнала. Тогда разницы не будет точно.

В речи например тоже огибающая медленно изменяется но для её сжатия используют линейное предсказание.

Я так понимаю линейное предсказание используют из-за ограниченности времени, так как сжатие речи идет для непосредственной передачи в линию и узнать, что опорную точку нужно поставить через несколько десятков миллисекунд просто нереально.


Я вечером попробую сжать не самым оптимальным образом и восстановить, посмотрим что получится.

[thermit]

Taradov Alexander:
Я так понимаю линейное предсказание используют из-за ограниченности времени, так как сжатие речи идет для непосредственной передачи в линию и узнать, что опорную точку нужно поставить через несколько десятков миллисекунд просто нереально.

Нет. Просто на сегодняшний день это один из наиболее эффективных способов параметризации случайных сигналов.
У меня есть подозрение, что именно lpc очень подойдет для Вашей задачи.
Что касается аппроксимации картинки из начала топика - она должна неплохо приближаться конечной суммой экспонент
вида c(k)*exp(-b(k)*n*T) k = 0 ... N-1.

[GetSmart]

Что касается аппроксимации картинки из начала топика - она должна неплохо приближаться конечной суммой экспонент
вида c(k)*exp(-b(k)*n*T) k = 0 ... N-1.

Если это Фурье-коэффициенты, то имхо ключевыми точками сплайнов можно сделать короче.

[Tanya]

Да, я об этом думал, но решил оставить на потом, после сжатия. Кроме того началный момент атаки скорее всего так описать не получится.

Почему же... Комбинация возрастающего до полочки + затухание.
Я Вас не пойму. Если есть (мечтать не вредно...) адекватное описание (модель - дифференциальное уравнение) с некоторым количеством параметров, то это и есть сжатие... А аппроксимировать можно 1000+ разными способами. А в случае существенной нелинейности толку от этого...
Было бы интересно посмотреть, как ведет себя все это при разных начальных энергиях.

[thermit]

GetSmart:
Если это Фурье-коэффициенты, то имхо ключевыми точками сплайнов можно сделать короче.

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

[ataradov]

Если есть (мечтать не вредно...) адекватное описание (модель - дифференциальное уравнение) с некоторым количеством параметров, то это и есть сжатие...

Ну да, только адекватной модели нет.

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

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

Я прмерно представляю на что нужно смотреть, просто времени это требует много и не все сразу

[Andrey_1]

Наверное, не только в этом случае. Если верить авторскому рисунку, кроме затухания происходит перенос энергии основной моды в гармоники. И обратно?

Это очень умная мысль - как раз такой "перенос" и возникает когда нет желания выполнять частотный анализ в режиме ZOOM - алаверды Бобу Рэндаллу Частотный Анализ 1989 год издание фирмы Брюль и Къер

[igorchem]

В идеале нужен алгоритм, который позволит по заданному числу контрольных точек определть их положение так, чтобы СКО было минимальным. Или наоборот, по заданному СКО определить требуемое число точек и их положение.

Если Вам необходимо найти оптимальные распределения точек для сплайн аппроксимации k-ой степени в l2 норме, можно поступить так:
вначале аппроксимируете Вашу функцию сплайном k+2 степени s_{k+2}(x).
Далее ищите разбиение x_0 < ... < x_i < ... < x_n так, чтобы
\int_{x_i}^{x_{i+1}} \left(\frac{\delta^{k+1} s_{k+2}(x)}{\delta x^{k+1}}\right)^{\frac1{k+1}} dx было одинакого на всех отрезках. Если запутаетесь или не сможете сами, постораюсь формулы выложить, но, очень надеюсь, что сами справитесь.

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

Моя идея - запишите эти функции (n штук) вначале на очень мелкой сетке с кусочно-постоянными или кусочно-линейными базисными функциями, чтобы шаг такой сетки был существенно меньше минимальной осцилляции. У Вас будет n векторов, образуйте из них матрицу. Сделайте этой матрице сингулярное разложение и выбросьте все сингулярные вектора, сингулярные значения которых будут меньше точности оцифровки вашого сигнала. Примените то, что я описал к оставшимся левым сингулярным векторам. Есть большая уверенность, что после всего этого Вы сильно сожмете эти данные, распаковка их будет съедать только несколько операций на точку, да и шум Вы тоже погасите.

Сообщение отредактировал igorchem - Mar 3 2011, 09:37

[Andrey_1]

Если Вам необходимо найти оптимальные распределения точек для сплайн аппроксимации k-ой степени в l2 норме, можно поступить так:
вначале аппроксимируете Вашу функцию сплайном k+2 степени s_{k+2}(x).
Далее ищите разбиение x_0 < ... < x_i < ... < x_n так, чтобы
\int_{x_i}^{x_{i+1}} \left(\frac{\delta^{k+1} s_{k+2}(x)}{\delta x^{k+1}}\right)^{\frac1{k+1}} dx было одинакого на всех отрезках. Если запутаетесь или не сможете сами, постораюсь формулы выложить, но, очень надеюсь, что сами справитесь.

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

Моя идея - запишите эти функции (n штук) вначале на очень мелкой кусочно-постоянной или кусочно линейной сетке, чтобы шаг такой сетки был существенно меньше минимальной осцилляции. У Вас будет n векторов, образуйте из них матрицу. Сделайте этой матрице сингулярное разложение и выбросьте все сингулярные вектора, сингулярные значения которых будут меньше точности оцифровки вашого сигнала. Примените то, что я описал к оставшимся левым сингулярным векторам. Есть большая уверенность, что после всего этого Вы сильно сожмете эти данные, распаковка их будет съедать только несколько операций на точку, да и шум Вы тоже погасите.

help polyfit
POLYFIT Fit polynomial to data.
P = POLYFIT(X,Y,N) finds the coefficients of a polynomial P(X) of
degree N that fits the data Y best in a least-squares sense. P is a
row vector of length N+1 containing the polynomial coefficients in
descending powers, P(1)*X^N + P(2)*X^(N-1) +...+ P(N)*X + P(N+1).

[P,S] = POLYFIT(X,Y,N) returns the polynomial coefficients P and a
structure S for use with POLYVAL to obtain error estimates for
predictions. S contains fields for the triangular factor ® from a QR
decomposition of the Vandermonde matrix of X, the degrees of freedom
(df), and the norm of the residuals (normr). If the data Y are random,
an estimate of the covariance matrix of P is (Rinv*Rinv')*normr^2/df,
where Rinv is the inverse of R.

[P,S,MU] = POLYFIT(X,Y,N) finds the coefficients of a polynomial in
XHAT = (X-MU(1))/MU(2) where MU(1) = MEAN(X) and MU(2) = STD(X). This
centering and scaling transformation improves the numerical properties
of both the polynomial and the fitting algorithm.

Warning messages result if N is >= length(X), if X has repeated, or
nearly repeated, points, or if X might need centering and scaling.

Class support for inputs X,Y:
float: double, single

See also poly, polyval, roots, lscov.

[igorchem]

help polyfit
POLYFIT Fit polynomial to data.
...

Уважаемый Андрей,

очень хочу поинтересоваться смыслом Вашего ответа на мой топик, мне не понятно несколько моментов:

1. зачем копировать текст от матворкса без ссылки - Ваша копия получилась с не очень красивой визуализацией математических формул, поэтому, на мой вгляд, было бы правильннее дать ссылку на полифит, например так:
http://www.mathworks.com/help/techdoc/ref/polyfit.html

2. как полифит согласуется с вопросом топикстартера?

3. как полифит может быть ответом на мой текст об оптимальном разбиении сплайн-аппроксимации?

4. Вы хоть раз пробовали вписывать такие функции как у топикстартера в полиномы высоких степеней сами, и видели, какая ошибка хоть в l1, хоть в l2 или l_{\infty} получается не говоря об устойчивости такого вписывания? Если нет, то не советуйте то, что не знаете, а если да, очень рад буду от Вас увидеть такие результаты (мне даже достаточно ссылки на Вашу публикацию в нормальном цитируемом журнале), с радостью просвещусь Вашими знаниями!

С уважением

ИИ

[ataradov]

очень хочу поинтересоваться смыслом Вашего ответа на мой топик, мне не понятно несколько моментов:

Не обращайте внимания, человек то-ли самоутверждается, то-ли посты копит для становления своим.

По теме: в формуле с помощью \delta - это обычная производная так странно записана или это какая-то особая производная?

Выбор числа k: при степенях полинома больше 10 получается уже не очень хорошо, много осциляций.

Каков физический смысл подинтегрального выражения?

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

[igorchem]

Написал Вам ответ и прикрепил свой первый пост в файле

http://www.elegant-mathematics.com/electronix.pdf

через пару дней - уберу

[ataradov]

Написал Вам ответ и прикрепил свой первый пост в файле

Спасибо, попробую разобраться.

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

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

Соотетственно затухающую часть можно не сжимать, а порсто заново восстанавливать при помощи простого банка фильтров. Получается естественно не 1:1 исходный графиг, но и фортепиано все звучат немного по-разнму, так что такое решение даже лучше простого сжатия.

Сложный харакер огибающей на приведенном графике обусловлен биениями хора струн (нота С4, 3 струны). Как именно в такой модели физично сделать биения я пока не знаю, но работаю над этим.

[igorchem]

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

Соотетственно затухающую часть можно не сжимать, а порсто заново восстанавливать при помощи простого банка фильтров. Получается естественно не 1:1 исходный графиг, но и фортепиано все звучат немного по-разнму, так что такое решение даже лучше простого сжатия.

Понятно, и, кстати, очень разумно. Иначе надо было бы огород городить или заумные пакеты пользовать.

У Вас действительно стоит задача много (100 и более) схожих семплов саппроксимировать? Если да, попоробуйте сингулярное разложение. Реально поможет! Если есть трудности - пишите в личку, помогу