Придумал алгоритм интерполяции. Протестируем результаты?, вызов от дилетанта Опции

[_Ivana]

Пару дней назад пришлось вникнуть в азы ЦОС на примере ресемплинга. Почитал пару-тройку страничек в инете, про двойной апсемплинг, оконные функции, полиномы Лагранжа-Фарроу и т.п. И придумал свой алгоритм Без таблиц, на чистой математике. По объему операций существенно меньше Фарроу. Реализовал его, как обычно, сначала на 1С Посмотрел графики. Потом перевел в целочисленную математику и с помощью коллеги наваял программку - ресемплер wav файлов. С заданием частоты ресемплинга и выбором варианта - своего или Фарроу.

Собственно, предлагаю: заинтересованные лица выкладывают вавки, я их ресемплю с нужной частотой двумя вариантами и выкладываю на всеобщее скачивание/заслушивание. Желающие делятся своими мнениями по поводу

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

Графические примеры работы алгоритмов: точки 50, 70, 20, края диапазона добиты нулями. График - ось абсцисс вниз, ординат - вправо.
Фарроу:

CODE

      *--- 0
     *
    *
   *
   *
  *
  *
   *
   *
    *
      *--- 0
          *
              *
                   *
                        *
                             *
                                   *
                                        *
                                             *
                                                   *
                                                        *--- 50
                                                            *
                                                               *
                                                                  *
                                                                     *
                                                                        *
                                                                          *
                                                                           *
                                                                            *
                                                                            *
                                                                            *--- 70
                                                                        *
                                                                    *
                                                               *
                                                          *
                                                     *
                                                *
                                          *
                                    *
                               *
                          *--- 20
                      *
                   *
                 *
              *
            *
           *
         *
        *
      *
      *--- 0
     *
     *
    *
    *
    *
    *
     *
     *
     *

Мой алгоритм:

CODE

      *--- 0
     *
     *
    *
   *
  *
  *
  *
  *
   *
      *--- 0
         *
             *
                  *
                       *
                             *
                                   *
                                         *
                                               *
                                                   *
                                                        *--- 50
                                                           *
                                                               *
                                                                  *
                                                                     *
                                                                        *
                                                                          *
                                                                            *
                                                                            *
                                                                            *
                                                                            *--- 70
                                                                         *
                                                                      *
                                                                 *
                                                           *
                                                     *
                                               *
                                         *
                                   *
                              *
                          *--- 20
                      *
                   *
                 *
              *
            *
           *
         *
        *
       *
      *--- 0
     *
    *
    *
    *
    *
     *
     *
     *
     *

За качество обоих графиков извиняйте - они в псевдографике с грубым "разрешением", на самом деле они более красивые

Сообщение отредактировал _Ivana - Apr 7 2012, 12:31

[ViKo]

Не логичнее выложить сам алгоритм?

[_Ivana]

Я бы не хотел его обнародовать. По крайней мере, пока не увижу его описание где-либо

[litv]

заслушивание на слух это - для музыкантов. Вы возьмите спектр от Вашего примера( до и после ресемплирования) и от Фароу например .

[_Ivana]

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

ОФФ - увидел "Воронеж" подписи и решил себе тоже написать

Сообщение отредактировал _Ivana - Apr 7 2012, 13:37

[des00]

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

ИМХО на глаз и на слух это субъективная оценка, надо сравнить системы (для интерполяторов лагранжа, в том числе по структуре Фарроу, накладываются определенные ограничения применения), построить частотные характеристики интерполятора при разных задержках, оценить ошибки интерполяции и т.д. Без алгоритма, по звуковому файлу это сделать крайне затруднительно.

Да и фарроу в звуковой аппаратуре лучше заменить банком фир фильтров (если вычислительный ресурс позволяет).

[_Ivana]

des00 спасибо за задание направления дальнейшего анализа. Чувствую, придется самому разбираться с перечисленными вами методами анализа Я наивно полагал что вавка может дать информацию. Хотя если кого заинтересует - предложение о тестовом сравнительном ресемплинге в силе.

И конечно Фарроу не самый лучший вариант. Я его беру только как сравнительный пример: мой алгоритм в принципе другой, но содержит меньше операций чем Фарроу. Если он окажется лучше последнего по результатам, мне будет приятно

[des00]

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

найдите в сети документ " Maximally Flat FD FIR Filter: Lagrange Interpolation" (первая строка гугла ) по нему сможете понять какими характеристиками описывается интерполятор лагранжа.

Если он окажется лучше последнего по результатам, мне будет приятно

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

[_Ivana]

Спасибо. Если кому есть что добавить по методам оценки алгоритмов, буду только рад. Я разберусь с методами оценки и анализа интерполяторов, протестирую свой алгоритм и не буду скрывать результаты

[thermit]

_Ivana:
По объему операций существенно меньше Фарроу.

1. Фарроу - структура интерполятора лагранжа с минимальным числом умножителей.

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

Сообщение отредактировал thermit - Apr 7 2012, 15:19

[Tanya]

Спасибо. Если кому есть что добавить по методам оценки алгоритмов, буду только рад. Я разберусь с методами оценки и анализа интерполяторов, протестирую свой алгоритм и не буду скрывать результаты

Вот думаю-думаю, а понять не могу. Если есть N точек, то через них можно "провести" единственную кривую, описываемую полиномом N-1 порядка. В чем Ваша альтернативность?

[_Ivana]

Я не знаю как прокомментировать последние реплики. Во-первых, я нигде не утверждал что моя интерполяция полиномиальна. Хотя и обратного я тоже не утверждал И таки да, по объему вычислений она проще Фарроу.
К тому же, в первом посте я привел графики. Пусть и грубые, но там видно, что моя интерполяция не совпадает с Фарроу и не является ей, а значит и Лагранжем 3-го порядка как минимум. Если интересует проверка на других исходных данных - выложите ваши точки, покажу новые графики.

Лучше подскажите нормальный спектроанализатор вавок (из файла на не со входа звуковухи), и чтобы весил не как СаундФордж. А если ещё и будет выдавать численные коэффициенты - вообще было бы хорошо. Сразу выложу вам цифры и спектры ресемплинга простых и не очень простых сигналов. Мой GoldWave не видит разницы при ресемплинге на +500Гц (в оцифровку 44600) как на синусе 500Гц (оцифровка 44100) так и на сумме синусов с биением частоты в 1%.

Сообщение отредактировал _Ivana - Apr 7 2012, 17:56

[sup-sup]

Я не знаю как прокомментировать последние реплики. Во-первых, я нигде не утверждал что моя интерполяция полиномиальна. Хотя и обратного я тоже не утверждал И таки да, по объему вычислений она проще Фарроу.
К тому же, в первом посте я привел графики. Пусть и грубые, но там видно, что моя интерполяция не совпадает с Фарроу и не является ей, а значит и Лагранжем 3-го порядка как минимум. Если интересует проверка на других исходных данных - выложите ваши точки, покажу новые графики.

Лучше подскажите нормальный спектроанализатор вавок (из файла на не со входа звуковухи), и чтобы весил не как СаундФордж. А если ещё и будет выдавать численные коэффициенты - вообще было бы хорошо. Сразу выложу вам цифры и спектры ресемплинга простых и не очень простых сигналов. Мой GoldWave не видит разницы при ресемплинге на +500Гц (в оцифровку 44600) как на синусе 500Гц (оцифровка 44100) так и на сумме синусов с биением частоты в 1%.

SpectraPlus.
Действительно, спектр рассудит.
Можно найти в 'радиосканнере'

Сообщение отредактировал sup-sup - Apr 7 2012, 18:59

[thermit]

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

2. По вашим рисункам лично я вообще не сказал бы, что это результат какой-либо интерполяции.

3. Гениеф дофига. Нормальных инженеров не хватает.

[_Ivana]

Спасибо, скачал, поставил. Сгенерировать ключ не удалось, использую 30-дневную версию. Действительно много возможностей. Навскидку, что увидел:
Исходный синус 500Гц (44100) THD = 0.00063, IMD = 0.1798
ресемплится в 44600.
Мой алгоритм: THD = 0.00124, IMD = 0.2349
Фарроу: THD = 0.00109, IMD = 0.2349
Если ничего не напутал

Биения 500Гц + 1%
Исходный сигнал THD = 0.00054, IMD = 0.6593
Мой алгоритм: THD = 0.00180, IMD = 0.6442
Фарроу: THD = 0.00125, IMD = 0.6443

Термит, если нечего сказать конструктивно (кроме очевидной методики проверки отклонения от целевого сигнала). то можно просто не писать в этот топик?

2 картинки спектра ресемплинга биений 500Гц+1% (ресемплинг на 500Гц). Первая - мой алгоритм, вторая - Фарроу:

Принимаются пожелания какие ещё тестовые сигналы ресемплить и посмотреть, спектр, какие может показательные коэффициенты (сигнал/шум там или что ещё из имеющегося в СпектраПлюсе)

Сообщение отредактировал _Ivana - Apr 7 2012, 19:29

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

 

[sup-sup]

Спасибо, скачал, поставил. Сгенерировать ключ не удалось, использую 30-дневную версию. Действительно много возможностей. Навскидку, что увидел:
Исходный синус 500Гц (44100) THD = 0.00063, IMD = 0.1798
ресемплится в 44600.
Мой алгоритм: THD = 0.00124, IMD = 0.2349
Фарроу: THD = 0.00109, IMD = 0.2349
Если ничего не напутал

Биения 500Гц + 1%
Исходный сигнал THD = 0.00054, IMD = 0.6593
Мой алгоритм: THD = 0.00180, IMD = 0.6442
Фарроу: THD = 0.00125, IMD = 0.6443

Термит, если нечего сказать конструктивно (кроме очевидной методики проверки отклонения от целевого сигнала). то можно просто не писать в этот топик?

2 картинки спектра ресемплинга биений 500Гц+1% (ресемплинг на 500Гц). Первая - мой алгоритм, вторая - Фарроу:

Принимаются пожелания какие ещё тестовые сигналы ресемплить и посмотреть, спектр, какие может показательные коэффициенты (сигнал/шум там или что ещё из имеющегося в СпектраПлюсе)

Еще можно по спектрограмме оценить для сложного сигнала на разные артефакты. В SpectraPlus спектрограммы так себе. Самые лучшие, что я видел, делает SA (в том же 'радиосканнере' есть. Демо-версии достаточно.
А ресемплить сигналы со сложным спектром, меняющимся во времени. Еще, близким к частоте Найквиста. Если алгоритм позволяет это.

Сообщение отредактировал sup-sup - Apr 7 2012, 19:51

[_Ivana]

А ресемплить сигналы со сложным спектром, меняющимся во времени. Еще, близким к частоте Найквиста. Если алгоритм позволяет это.

Например какие? ГолдВэвйв много разных может генерить, СпектраПлюс поменьше вариантов, но тоже есть кое-что. Есть синус умноженный на экспоненту и много чего ещё диковинного, плюс самому можно формулы написать

Частотная модуляция . 500+-Гц. Ресемплинг тот же - с 44100 до 44600.
Мой THD = 0.00091, IMD = 0.6729, слева
Фарроу THD = 0.00112, IMD = 0.6729, справа

Сообщение отредактировал _Ivana - Apr 7 2012, 20:06

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

 

[thermit]

_Ivana:
Термит, если нечего сказать конструктивно (кроме очевидной методики проверки отклонения от целевого сигнала). то можно просто не писать в этот топик?

С начала этого топика кроме меня только 1 (один) человек ответил вам конструктивно. Пальцем тыкать не хочу.
Это что касается конструктивизма в частности.

Что касается конструктивизма вообще:
1 Да. Я хочу увидеть не непонятные картинки, а пиковую ошибку интерполяции гармонического сигнала в зависимости от частоты при частоте дискретизации = 1. Все остальное не несет никакой полезной информации.

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

в этой связи

3 Обнародуйте свой алгоритм. едва-ли этим вы уроните кого-либо под стол.

[sup-sup]

Например какие? ГолдВэвйв много разных может генерить, СпектраПлюс поменьше вариантов, но тоже есть кое-что. Есть синус умноженный на экспоненту и много чего ещё диковинного, плюс самому можно формулы написать

Частотная модуляция . 500+-Гц. Ресемплинг тот же - с 44100 до 44600.
Мой THD = 0.00091, IMD = 0.6729, слева
Фарроу THD = 0.00112, IMD = 0.6729, справа

Я предложил (по методике тестирования) для визуальной оценки смотреть спектрограммы в хорошем разрешении. Поиграться при этом размером FFT для нахождения лучшего разрешения по частоте и по времени. Это даст не усредненный замазанный результат, а более пиковый, как сказал thermit. Из сигналов можно взять ЛЧМ (chirp) от нулевой частоты до fs/2. Еще, сумму нескольких синусов в верху диапазона с разными фазами. И вообще, искать именно тот сигнал, который плохо интерполируется. В этом случае вопросов и не будет. Ищите неисправности в алгоритме и, если не находите, все хорошо. И смотреть в динамике, то есть спектрограммой. Это для качественной оценки одним взглядом. Более интересные места смотреть подробнее во временной и в частотной области.

Сообщение отредактировал sup-sup - Apr 7 2012, 20:30

[_Ivana]

Хорошо. Термит прав, надо отойти от цифровых применений и просто чисто математически проинтерполировать синус с разной частотой, при этом посчитать максимальный модуль отклонения от чистого синуса на интервале и средний модуль отклонения на нем же. Я и сам думал так сделать. Жаль в 1С нет функции синуса, придется как-то выкручиваться. Но я сделаю это, для нескольких частот точно. И сравню с Фарроу. И выложу сюда. Как найду где взять точный чистый синус любого аргумента, простите за банальные трудности

ЗЫ а насчет гениев и инженеров у меня своя версия - мало и тех и других. Критикофф с претензиями зато дофига

[thermit]

Да, вопрос от меня, как от критика с претензиями вам как к гению : При чем тут 1с?

[_Ivana]

ОФФ: потому что я её хорошо знаю, и использую как тестовый инструмент для отладки разной математики. Она отлично выводит и таблицы и графики и псевдографики и я могу удобно просматривать любые промежуточные результаты в любом виде и формате. И всем рекомендую
Например, для первого поста потребовалась псевдографика с хитрыми непечатными символами а-ля пробел, ибо пробелы с краев отрезаются. В 1С я это сделал за минуту.

ЗЫ я понимаю что критики с претензиями предпочитают наверное матлаб или маткад, и я с ними согласен, но его надо ещё найти, поставить и освоить.

[thermit]

_Ivana:
ЗЫ я понимаю что критики с претензиями предпочитают наверное матлаб или маткад, и я с ними согласен, но его надо ещё найти, поставить и освоить.

Угу. Предпочитают. Дык для гениев еще ексель есть... Там даже синус имеецца. Не говоря уж про тангенсы всякие...

Для гениев чуть меньшего размера осмелюсь порекомендовать это изделие www.jsoftware.com (v6.x).

Сообщение отредактировал thermit - Apr 7 2012, 21:26

[_Ivana]

За рекомендации всегда спасибо. Но от специалистов они всегда ценнее и эффективнее. Я уже нашел интегрируемый в 1С вариант:

CODE

// синус
function sin(value,sc=0)
try
if sc=0 then
sc=createObject("MSScriptControl.ScriptControl");
endif;
sc.language="VBscript";
except
return getEmptyValue();
endtry;
return sc.eval("sin("+value+")");
endFunction

Вроде работает, график на глаз гладкий. Буду пробовать с этой функцией.

[Pavel_SSS]

ОФФ: потому что я её хорошо знаю, и использую как тестовый инструмент для отладки разной математики. Она отлично выводит и таблицы и графики и псевдографики и я могу удобно просматривать любые промежуточные результаты в любом виде и формате. И всем рекомендую
Например, для первого поста потребовалась псевдографика с хитрыми непечатными символами а-ля пробел, ибо пробелы с краев отрезаются. В 1С я это сделал за минуту.

ЗЫ я понимаю что критики с претензиями предпочитают наверное матлаб или маткад, и я с ними согласен, но его надо ещё найти, поставить и освоить.

Интересно, попробовать вести бухгалтерию в Матлабе... Что-то в этом есть...
Насчет интерполяции - как тут уже отмечали, интерполяция - это нахождение непрерывной функции, желательно выраженной с помощью простых действий - умножения, деления, сложения вычитания, может быть еще извлечения корня, проходящей через заданные точки. Имея такую функцию, мы можем найти ее значения, а, следовательно, значения входного сигнала в промежуточных точках. Насколько значение функции будет соответствовать входному сигналу? Неизвестно. Теоретически возможно, что входной сигнал между отсчетами улетает до небес. Поэтому мы должны сделать предположение, что сигнал - физическая величина, следовательно обладает конечной энергией и конечным - и весьма нешироким, спектром. Иначе неуместно говорить об интерполяции сигнала, а надо говорить об математической абстракции - провести кривую через столько-то точек. Удобным методом интерполяции является полиномиальная интерполяция - нахождение полинома степени N-1, выражающего функцию, проходящую через N точек. Лагранж предложил метод нахождения таких полиномов, а Фарроу - способ построения интерполятора оптимально использующего вычислительные мощности. Но по сути - это полиномиальная интерполяция, дающая одинаковый результат независимо от способа вычисления полинома, так как такой полином существует один. Такая интерполяция хорошо подходит, в частности, к музыке, так как музыкальный сигнал - это некий набор гармонических составляющих - синусоид, параметры которых меняются относительно медленно. Эти синусоиды можно разложить в ряд Тейлора, сложить разложения, и получим, как раз, многочлен, причем невысокой степени.
Есть случаи, когда полиномиальная интерполяция малопригодна - например, при попытке интерполировать изображение, оно может потерять четкость переходов...

[ViKo]

Я думаю, что максимально верной интерполяцией нужно считать интерполяцию синусную, когда для каждой существующей выборки создают функцию sin(x)/x, с вершиной в этой выборке (а для остальных выборок эта функция будет в нуле), и затем складывают все эти функции. Дальше по полученной кривой вычисляют значения в любых других точках. Но это слишком сложный алгоритм.

Между выборками сигнал не должен меняться, как угодно, иначе частота дискретизации не будет удовлетворять теореме Котельникова.

[Tanya]

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

любую.
Даже аппроксимацию при бесконечном числе точек нельзя выбрать. Пока не выбран критерий близости функций.
Вот попробуйте - (0, 0), (1, 1), (100, 0), (101, -1). Можете продлить...
И.... главное - докажите, что Ваша интерполяция "максимально верная".

[ViKo]

Вот попробуйте - (0, 0), (1, 1), (100, 0), (101, -1). Можете продлить...
И.... главное - докажите, что Ваша интерполяция "максимально верная".

Имею в виду равномерную дискретизацию. Доказательство элементарное. Идеальный ФНЧ имеет именно такую импульсную характеристику.

А вот, собственно, доказательство. Первая же формула.
http://ru.wikipedia.org/wiki/%D0%9F%D0%B5%...%86%D0%B8%D1%8F

[Tanya]

А вот, собственно, доказательство. Первая же формула.

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

[_Ivana]

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

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

 

[petrov]

Я думаю, что максимально верной интерполяцией нужно считать интерполяцию синусную, когда для каждой существующей выборки создают функцию sin(x)/x, с вершиной в этой выборке (а для остальных выборок эта функция будет в нуле), и затем складывают все эти функции. Дальше по полученной кривой вычисляют значения в любых других точках. Но это слишком сложный алгоритм.

Он вообще не реализуемый, а вот обрезанный синк является ли самым лучшим ФНЧ из возможных?

[thermit]

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

Ну, это уже лучше. Теперь огласите число операций на выходной отсчет вашего интерполятора.

[_Ivana]

Понятно что лучше
Число операций на отсчет оглашу чуть позже. Картинки с графиками максимального и среднего отклонений тоже уже получил - вышло немного забавно Сейчас формирую полные графики отклонений с учетом любой возможной фазы интерполируемого синуса на любой частоте (сдвиг по фазе от 0 до 2пи внутри любой частоты). Графики опять же чуть позже.
А пока ТИПИЧНЫЕ (НЕ СЛУЧАЙНЫЕ) картинки сравнения интерполяций (мой зеленый, Фарроу синий):
Видна некая тенденция. Которая, кстати говоря, и позволяет мне предполагать, что по некоторым критериям мой алгоритм может оказаться явно лучше

Собственно, искомые графики Красный - мой, зеленый - Фароу
По оси абсцисс - отношение частоты дискретизации к частоте исходного синуса. Ссылаясь на Котельникова, я эту величину не брал меньше 2 Начинается с 2 и возрастает до.... можете на графике посмотреть.
Первый график - максимальный модуль отклонения (в абсолютных единицах, при синусе от -1 до +1). Второй - средний модуль отклонения.

ЗЫ собственно, недавно кто-то говорил, что "Интерполяция - большой раздел вычислительной математики" (С). Вот видимо эта умная математика и привела в итоге к таким результатам сравнения алгоритмов

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

 

[thermit]

_Ivana:
собственно, недавно кто-то говорил, что "Интерполяция - большой раздел вычислительной математики" (С). Вот видимо эта умная математика и привела в итоге к таким результатам сравнения алгоритмов

Пока что вы продолжаете плодить малоинформативные картинки. В тех условиях, что вы задали (4 точки на период) любой алгоритм интерполяции будет работать плохо. Что и видно на ваших картинках.
Задайте хотя бы 6 т/период, а лучше 8. А еще лучше выполнить измерение максимальной погрешности интерполяции для ряда частот 1/64 1/32 1/16 1/8 и выложить тут график.

[_Ivana]

thermit, странная ситуация
Кроме Вас, в этой теме похоже никого это не интересует и никто не хочет хоть какие-то графики от меня получить. Получается, что только Вы (при всем Вашем неоправданном снисходительно-пренебрежительном отношении) и пытаетесь разобраться. Так разбирайтесь, читайте внимательнее.
ПОСЛЕДНИЕ ГРАФИКИ - это И ЕСТЬ максимальная и средняя погрешности ОТ ОТНОШЕНИЯ ЧАСТОТЫ ДИСКРЕТИЗАЦИИ К САМОЙ ЧАСТОТЕ СИНУСА - Вашими словами, количества точек на период синуса!

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

Сообщение отредактировал _Ivana - Apr 8 2012, 12:23

[GetSmart]

_Ivana, для начала не нужно выбирать предельную дискретизацию. Мне видится на графиках 2 точки на период, то бишь почти по пределу Котельникова. Фарроу на такой частоте реально плохо работает. Синк, кстати, неплохо, обрезанный естественно.

Хотите чем-то возгордиться, покажите графики с 3-мя или 4-мя точками на период синуса.

[_Ivana]

Картинки при 2 точках на период (предел по Котельникову) я показал ТОЛЬКО с целью наглядной демонстрации различия алгоритмов. И сам график отклонений логично строить от значения 2 - предела по Котельникову,я уже об этом писал.
Вот картинки при повышении частоты дискретизации, по 2 - с разным фазовым сдвигом. Они не такие наглядные, поэтому я их и не показывал. При дальнейшем повышении частоты дискретизации будет еще ненагляднее

Сообщение отредактировал _Ivana - Apr 8 2012, 12:51

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

 

[thermit]

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

Ошибка при 2-4 точках на период маскирует ошибки для меньших частот. Поэтому ваши картинки по-прежнему малоинформативны.

Пора огласить число операций.

[GetSmart]

Картинки при 2 точках на период (предел по Котельникову) я показал ТОЛЬКО с целью наглядной демонстрации различия алгоритмов. И сам график отклонений логично строить от значения 2 - предела по Котельникову,я уже об этом писал.

НЕ ЛОГИЧНО, поверьте. Это крайний случай, который информативно ущербен. Иногда, очень редко, он может иметь значение. Обычно разумный проектировщик системы к нему (пределу) не приближается.

Поэтому демонстрируйте свой алгоритм во всей красе в широком диапазоне дискретизаций. Цифры в ранних постах пока никакой заметной относительно фэрроу точности не показали. По поводу кол-ва вычислений до сих пор загадка.

[_Ivana]

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

Хотите, построю этот же график но не с 2 а с 5 или 6 точек на период? Это будет тот же график, только отрезанный и отмасштибированный. Со скольки точек построить?

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

[thermit]

Хотите, построю этот же график но не с 2 а с 5 или 6 точек на период? Это будет тот же график, только отрезанный и отмасштибированный. Со скольки точек построить?

C 8-ми точек на период. Естественно, растянутый по оси y.

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

Не надо туман напускать. Умножение, сложение, деление, извлечение корня - операции.
Сколько и каких именно подобных операций на 1 выходной отсчет требует ваш интерполятор?

Сообщение отредактировал thermit - Apr 8 2012, 13:02

[_Ivana]

C 8-ми точек на период. Естественно, растянутый по оси y.

Наконец-то разумные пожелания. Сейчас построю, прикреплю.

Не надо туман напускать. Умножение, сложение, деление, извлечение корня - операции.
Сколько подобных операций на 1 выходной отсчет требует ваш интерполятор?

А почему бы мне не напустить тумана в этом вопросе? Я ведь пока не хочу обнародовать алгоритм Но повторяю в 100500-й раз - гарантированно меньше Фарроу при любых аппаратных возможностях При усеченных возможностях - ощутимо меньше

[ViKo]

Какое это доказательство? И чего? Ведь мы имеем некоторую функцию, заданную отсчетами в некоторых точках.
Потом мы придумываем непрерывную функцию, которая проходит через эти точки. Вариантов придумки бесконечно много. И ни один не лучше другого.

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

[thermit]

А почему бы мне не напустить тумана в этом вопросе? Я ведь пока не хочу обнародовать алгоритм sm.gif Но повторяю в 100500-й раз - гарантированно меньше Фарроу при любых аппаратных возможностях sm.gif При усеченных возможностях - ощутимо меньше

У нас есть такие приборы... Но мы вам про них не расскажем (цэ)
Зачем было тогда сюда писать?

[_Ivana]

Термит, Вы были правы - маскирует. Мой- красный, Фарроу - зеленый. И по моему (если я нигде ничего не напутал), это действительные показатели алгоритмов, а не накопленная ошибка округления или что-то ещё. Ну чтож - результат таков каков он есть Хотя я оставляю ещё шанс что я ошибся при построении графиков, но он честно говоря невелик
Спасибо вам за помощь в анализе (несмотря на остальные моменты)

У нас есть такие приборы... Но мы вам про них не расскажем (цэ)
Зачем было тогда сюда писать?

А вы не обижайтесь. В наш век Кали Юги на дворе не все всем все рассказывают
К тому же, в самом заголовке темы я предложил протестировать РЕЗУЛЬТАТЫ, а не говорил что буду обнародовать алгоритм

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

 

[thermit]

_Ivana:
В наш век Кали Юги на дворе не все всем все рассказывают

Гы. Я буду первым, кто скомуниздит ваш супералгоритм и сказочно обогатится...

[_Ivana]

Шутки шутками, а вот если бы моя кривая прошла ниже Фарроу, тогда уже был бы совсем другой разговор

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

Средняя ошибка. Тоже уступаю Фарроу....

Сообщение отредактировал _Ivana - Apr 8 2012, 13:48

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

 

[GetSmart]

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

Вот же ж идеалисты. Учебники прямо так и учат. Никаких нуансов.
Физику процесса надо бы понимать чтоб так категорично заявлять.
Если даже БПФ "затрудняется" точно посчитать амплитуду дробной частоты по _исходным сэмплам_.
То уж дубовый интерполятор и подавно не сможет (идеально точно и однозначно).
Не стоит забывать, что всегда кол-во сэмплов ограничено на практике.

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

Даже в общих чертах не проясните во сколько/насколько быстрее?

Сообщение отредактировал GetSmart - Apr 8 2012, 13:54

[_Ivana]

Даже в общих чертах не проясните во сколько/насколько быстрее?

В общих чертах... Ну например, при расчете Фарроу используется деление на 6 (хорошо, уговорили - на 3 ). Даже если есть аппаратное деление, то мой алгоритм быстрее ..... а вот не знаю насколько - но будет как раз протестировать! А уж если нет аппаратного деления - тогда за счет необходимости его эмуляции на асме получается ещё бОльшая разница. Кстати сказать, даже если бы Фарроу и не требовал этого пресловутого деления на 6, то мой алгоритм все равно бы был быстрее
Но больше не выпытывайте пожалуйста Если получится - протестирую на железе и скажу реальную разницу в скорости.

[des00]

В общих чертах... Ну например, при расчете Фарроу используется деление на 6 (хорошо, уговорили - на 3 )

может открою вам глаза, но у меня сплошь и рядом в модемах на ПЛИС используется интерполяторы лагранжа, по структуре Фарроу. Так вот там нет ни одного делителя совсем. И все прекрасно работает так что математика у вас не той системы

[thermit]

_Ivana:
Ну например, при расчете Фарроу используется деление на 6 (хорошо, уговорили - на 3 sm.gif ).

Деление на константу это умножение на константу. Так что в фарроу никаких деление нет в принципе.
Есть умножение на константу и перемножение 2-х операндов. Сумма есть.

[_Ivana]

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

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

[des00]

Но это не отменяет принципиального факта меньшего количества операций в моем алгоритме

т.е. вы делаете его быстрее чем 3 умножения и 4 сложения ?

[_Ivana]

т.е. вы делаете его быстрее чем 3 умножения и 4 сложения ?

Давайте разберемся. Вот Фарроу:

CODE

а3 = (1/6)*(у3 - у0) + (1/2)*(у1 - у2);
а1 = (1/2)*(у3 - у1) - а3;
а2 = у3 - у2 - а1 - а3;
а0 = у2;

х = т - 1;

у = а3*х;
у = (у + а2)*х;
у = (у + а1)*х;
у = у + а0;

Ваша фраза про "3 умножения и 4 сложения" относится к чему? Даже если оставить в стороне вычисление итогового значения полинома (что ещё 3 умножения и 4 сложения - одно на приведение параметра допустим), то получается 3 умножения и 8 сложений. Хотя, я не оптимизировал, может и можно ужать. Но если говорить о полном цикле операций Фарроу - то это 6 умножений и 12 сложений - если считать как я написал, не оптимизированно. Кто посчитает за меньшее? Прошу ваш вариант.

Сообщение отредактировал _Ivana - Apr 8 2012, 14:22

[des00]

Кто посчитает за меньшее? Прошу ваш вариант.

да весь интернет знает, только вы в неведении

На этом синтез фильтра Фарроу третьего порядка можно считать законченным и подведем некоторые итоги:
.....
2. Фильтр третьего порядка требует всего 3 операции умножения и может применяться в реальном времени.

[NiceParty]

да весь интернет знает, только вы в неведении

Прошу прощения, что вмешиваюсь в беседу уважаемых метров DSP.
На рисунке 7: "Модифицированный фильтр Фарроу третьего порядка" из приведенной Вами в пример статьи я насчитал 11 операций сложения и 6 операций умножения для итогового значения полинома в точке интерполяции.

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

P.S. Две операции умножения на 1/2, конечно, лучше заменить на арифметический сдвиг вправо (подразумевается целочисленная математика).
Но умножение на 1/6 я пока не знаю чем заменить. Буду признателен, если подскажете.

Сообщение отредактировал NiceParty - Apr 8 2012, 15:05

[_Ivana]

да весь интернет знает, только вы в неведении

des00 вы просто верите цитатам из интернета, даже не пытаясь проанализировать их смысл?
Если найдутся участники, понимающие о чем речь, а также понимающие разницу в количестве операций при умножении на 1/2 и 1/6 даже во флоатах, и при этом будут иметь желание продолжить дискуссию - я готов.

[Tanya]

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

Какая-такая -"Эта"?
Вы меня убиваете (без Котельникова).
Функция номер раз - кусочно-линейная.
Функция номер два - кусочно-квадратичная.
Дальше сами... Единственность опровергнута? Ваша.

[des00]

P.S. Две операции умножения на 1/2, конечно, лучше заменить на арифметический сдвиг вправо (подразумевается целочисленная математика).
Но умножение на 1/6 я пока не знаю чем заменить. Буду признателен, если подскажете.

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

в результате такой "моидификации", коэффициент усиления преобразования интерполятора будет 6. легко переводиться в 6/8. И совсем привередливым, кому хочеться получить 1цу, приводиться простым делением через умножение.

des00 вы просто верите цитатам из интернета, даже не пытаясь проанализировать их смысл?

а вы слишком узко смотрите на реализацию

[_Ivana]

а вы слишком узко смотрите на реализацию

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

[ViKo]

Какая-такая -"Эта"?
Вы меня убиваете (без Котельникова).
Функция номер раз - кусочно-линейная.
Функция номер два - кусочно-квадратичная.
Дальше сами... Единственность опровергнута? Ваша.

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

[thermit]

des00:
все переходит к набору сумматоров + сдвиги

Умножение на константу и без всяких извращений сведется к сдвигам/суммам. Дополнительная суета не требуется.

[NiceParty]

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

в результате такой "моидификации", коэффициент усиления преобразования интерполятора будет 6. легко переводиться в 6/8. И совсем привередливым, кому хочеться получить 1цу, приводиться простым делением через умножение.

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

P.S. Мне считать не лень, но я не смог увидеть преимуществ в кол-ве операций. Может, неправильно считал...

[Tanya]

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

1.Котельников никак и никого не учил оцифровывать.
2. Ни про какие частоты в разложении в ряд Фурье интерполяции нет никакого дела. Мухи отдельно. Ладно?
3. Изломы имеются только в точках пересечения прямых. Если новую функцию там не определять, то какие изломы?
Вы же сомневаетесь в их существовании для кусочно-квадратичной.
Так объясните, что есть та самая правильная интерполяция?
Лучше формулой. А то я Ваши слова плохо понимаю.

[_Ivana]

Подскажите кто знает, какие ещё показатели (кроме максимального и среднего модулей отклонений) и как я могу получить, зная отклонения моей интерполяции от целевой функции во множестве точек? Какая-нибудь там дисперсия и прочая вероятностно-статистические показатели? Что они каким-то образом характеризовали работу алгоритма. Я все ещё не теряю надежды выявить сильные/слабые стороны обоих методов

[Pavel_SSS]

Подскажите кто знает, какие ещё показатели (кроме максимального и среднего модулей отклонений) и как я могу получить, зная отклонения моей интерполяции от целевой функции во множестве точек? Какая-нибудь там дисперсия и прочая вероятностно-статистические показатели? Что они каким-то образом характеризовали работу алгоритма. Я все ещё не теряю надежды выявить сильные/слабые стороны обоих методов

А Вы знаете целевую функцию? Зачем Вам тогда интерполяция - просто вычислите значения функции в нужных точках...

[ViKo]

1. Котельников никак и никого не учил оцифровывать.
2. Ни про какие частоты в разложении в ряд Фурье интерполяции нет никакого дела. Мухи отдельно. Ладно?
3. Изломы имеются только в точках пересечения прямых. Если новую функцию там не определять, то какие изломы?
Вы же сомневаетесь в их существовании для кусочно-квадратичной.
Так объясните, что есть та самая правильная интерполяция?
Лучше формулой. А то я Ваши слова плохо понимаю.

1. Без комментариев
2. Есть очень тесная связь, для ЦОС.
3. Не понял... кусочно-линейная... состоит из кусков линейки()... в месте склейки кусков будут изломы... практически, в каждой точке исходной функции. В наличии изломов кусочно-квадратичной (из кусков квадратиков )... тоже не сомневаюсь. Говорю о том, это будет менее заметно (ближе к той, единственной), но дополнительные спектральные составляющие обязательно будут.
Вместо формул (я же уже дал ссылку, там русским по белому написано!) картинку покажу.

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

 

[sup-sup]

Подскажите кто знает, какие ещё показатели (кроме максимального и среднего модулей отклонений) и как я могу получить, зная отклонения моей интерполяции от целевой функции во множестве точек? Какая-нибудь там дисперсия и прочая вероятностно-статистические показатели? Что они каким-то образом характеризовали работу алгоритма. Я все ещё не теряю надежды выявить сильные/слабые стороны обоих методов

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

Сообщение отредактировал sup-sup - Apr 8 2012, 17:59

[ViKo]

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

[_Ivana]

Спасибо за дельные рекомендации, померю по возможности.

ЗЫ to all - вы меня только пожалуйста не бейте Я тут ещё почитал немного интернетов, и понял что я придумал так называемую "интерполяцию сплайнами Катмулла-Рома", точнее один из её вариантов И ещё и проанализировал её в сравнении с Лагранжем по максимальной и средней ошибке Я действительно не знал об этом раньше (я бы до таких фамилий никогда не додумался ), но раз дело такое - как обещал, выкладываю код, сравните по быстродействию с Фарроу:

CODE

SHORT InterpolateKatmullRom(int t, SHORT y0, SHORT y1, SHORT y2, SHORT y3)
{ int y, a0, a1, a2, a3;

y = y2 - y1;
a0 = y1;
a1 = (y2 - y0)>>1;
a3 = (y3 - y0);
a3 = (a3 - y)>>1;
a3 = (a3 - y);
a2 = y - a3 - a1;

y = a3*t;
y = y>>16;
y = (y + a2)*t;
y = y>>16;
y = (y + a1)*t;
y = y>>16;
y = y + a0;

if (y > 32767) y = 32767;
else if (y < -32768) y = -32768;

return (SHORT)y;
}

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

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

[ViKo]

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

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

[_Ivana]

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

Хорошо, пусть будет ограничение

ЗЫ есть у меня ещё одна идейка Надо бы её проверить теперь, благо, спасибо Термиту, я уже знаю как навскидку проверять разные алгоритмы интерполяции, и проверялка уже написана. Только хорошо бы выполнить 4 условия:
1) формализовать её сначала математически а потом в коде
2) победить Фарроу по отклонениям
3) оптимизировать по операциям чтобы было в разумных пределах
4) самое главное - чтобы это потом не оказалось методом какого-нибудь Фермана-Зингельшухера

Начинаю серию №2, надеюсь реализовать её тоже.

[Михаил_K]

_Ivana. Интересно было бы посмотреть на спектр сигнала полученного вашим интерполятором, при интреполяции в несколько раз, при условии на вход подается синус, частота которого близка к 1/2 частоты дискретизации. Например 0.49 частоты дискретизации. Оценить уровень аллиаса, который при этом получится.
А если серьезно...
Как то давно, на первом курсе института был у меня такой случай. Выполняли достаточно простенькое задание. Зову препода - типа сделал. Он подходит проверять. Вводит начальные значения. Получает результат. Вводит другие значения - программа отрабаывает не правильно. Ну я соответственно исправляю. Опять зову. Он опять вводит несколько вариантов исходных данных. Программа в очередной раз ломается. Я опять исправляю. Опять зову. Опять вводит несколько вариантов. Программа работает. Препод в задумчивости.. Пробует еще несколько вариантов. Программа работает. Препод задумчиво, со звуком "ХМ" лезет в код. Изучает внимательно. Вводит очередные исходные данные - программа ломается. На что он с криком "Ура" заставляет снова думать.
Мораль:
Вы бы описали свой алгоритм, а люди бы посмотрели на него, и сказали бы свое мнение. Не исключаю, что то, что вы придумали - уже известно, и имеет определенные границы применения, либо определенные "грабли".

[thermit]

http://electronix.ru/forum/index.php?showt...t&p=1046889

Собственно,мнение тут одно: Все уже давно придумано.

Сообщение отредактировал thermit - Apr 9 2012, 14:02

[_Ivana]

_Ivana. Интересно было бы посмотреть на спектр сигнала полученного вашим интерполятором, при интреполяции в несколько раз, при условии на вход подается синус, частота которого близка к 1/2 частоты дискретизации. Например 0.49 частоты дискретизации. Оценить уровень аллиаса, который при этом получится.

Запросто. Чуть позже, как будет время.

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

В данный момент размышления следующие. Фарроу (Лагранж 3-го порядка) - это имхо по сути (хотя многие могут со мной не согласиться) тот же кубический сплайн. Локальный, определяемый для 4-х точек, но действующий только на средний интервал из 3-х. Точно так же, как и сплайны Эрмита, например. Тот факт, что Лагранж при построении проходит через все 4 точки, не дает значения интерполирующей функции в левом и правом интервалах, поскольку при интерполяции собственно самих этих боковых интервалов (когда они являются центральными) полином проходит совершенно не так, как проходит в случае если интервал является правым или левым относительно центрального. Значит нельзя считать условие прохождение полинома через крайние точки самым оптимальным с точки зрения построения сплайна. Я придумал (а Катмулл и Ром до меня) другие условия: сплайн проходит через 2 точки и имеет определенные производные в них же, так чтобы полученная функция имела непрерывную первую производную. Сначала я придумал считать значение этой производной, как тангенс угла-биссектрисы исходных двух (словами трудно объяснить, надо на графике показывать). Только что проверил этот вариант - при малых дискретизациях ошибка выше Катмулла-Рома, при увеличении начиная где-то с 12 точек почти полностью сливается с Катмуллом-Ромом (график ниже). Но само значение производной рассчитывается через тангенсы-арктангенсы и достаточно сложно. Это и была моя вторая идея, которую я проверил только что. А потом я придумал считать производную как тангенс наклона соседних точек - собственно то, что потом узнал как Катмулла-Рома.

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

Собственно, на текущий момент мысли таковы. А снизу - график ошибки Фарроу, Катмулла-Рома и метода "биссектрисы" (как я его называю ) И обещанный ресемплинг синуса 3920Гц с исходных 8000 до 44100

Upd: добавлен график сравнения ошибок 4-х методов: зеленый - Катмулла-Рома, синий - Фарроу, желтый - полиномами Лагранжа 4 степени, красный - кубическими сплайнами при условии точного знания производной интерполируемой функции во точках отсчета. Получается, красный - теоретически возможный предел для сплайнов 3-го порядка с условиями по производным в точках отсчета. Но достаточно хороший предел, надо сказать, можно к нему и постремиться Осталось только придумать красивую и точную процедуру угадывания производных. Я придумал пока две - одна оказалась Катмуллом-Ромом, другая - сложнее неё и не лучше по результатам. Может придумаю ещё

Сообщение отредактировал _Ivana - Apr 9 2012, 18:40

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

 

[_Ivana]

В продолжение эпического мыслевыражения Задача: имеем 3 точки при равномерной дискретизации. Надо придумать, какую мы хотим производную интерполирующего сплайна в средней точке. Катмулл и Ром предложили просто - линия через крайние точки, средняя точка игнорируется - получается относительно большая ошибка. Я сейчас подумал - мы можем интерполировать эти 3 точки какой-нибудь красивой гладкой функцией, и посчитать значение её производной в средней точке. На ум сразу приходит парабола. Но, похоже, это и есть то, что имеется в виду под "сплайнами, заканчивающимися параболой" http://alglib.sources.ru/interpolation/spline3.php Куда не сунься - все уже украдено до нас! Реалистичный Термит похоже снова прав...

[Timmy]

гладкой функцией, и посчитать значение её производной в средней точке. На ум сразу приходит парабола. Но, похоже, это и есть то, что имеется в виду под "сплайнами, заканчивающимися параболой" http://alglib.sources.ru/interpolation/spline3.php Куда не сунься - все уже украдено до нас! Реалистичный Термит похоже снова прав...

В общем случае кубический сплайн можно построить сразу по большому числу точек, обеспечив непрерывность первой и второй производной во всех точках. Для этого требуется решение относительно сложной системы линейных уравнений. В "сплайнах, заканчивающихся параболой" для крайних точек дополнительно задаётся условие равенства нулю третьей производной, для определённости. Таким образом, крайние отрезки интерполируются полиномами второй степени, а не третьей.
Ещё замечу, что тест надо делать до 100 точек на период синусоиды, так как ошибка в 0.001 при 20 точках - это просто порнография. Для максимальной ошибки 10^-6 требуется около 50-60 точек при использовании полинома Лагранжа. И лучше использовать логарифмический масштаб по вертикали.

[_Ivana]

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

Про параболу и Катмулла-Рома: я с удивлением обнаружил, что если 3 точки интерполировать параболой, то производная в средней точке ВСЕГДА будет равна прямой через крайние Ну вот так Лагранж 2-го порядка проводит свои полиномы! Но я найду способ учесть положение средней точки, у меня уже есть идеи

А про 100 точек и логарифмический масштаб... Вы правы конечно, но пока я просто качественно проверяю разные алгоритмы, для оценки предпочтения одного перед другими, без особого упора на абсолютные значения, и надеюсь что кривые отклонений не пересекутся где-то там на 50-100 точках При более детальном количественном анализе конечно придется строить логарифмический масштаб и от 8 до 100 точек на период.

[_Ivana]

Пока все мои попытки придумать производную точнее чем Катмулл-Ром потерпели неудачу.
Но хотя бы выложу красивый график - как и просили, до 100 точек и в логарифмическом масштабе. Лагранж-4 обогнал-таки кубический сплайн с точными производными (если я нигде не ошибся). Да и сами графики какие-то подозрительно линейные, и весьма разнятся с представленными в недавней ссылке по сплайновой интерполяции... Хотя у меня и по оси абсцисс масштаб тоже логарифмический, а в тех графиках линейный, это может как-то объяснить красивую прямизну графиков Но сами значения у меня повыше, например для сплайна с точными производными.

Сообщение отредактировал _Ivana - Apr 9 2012, 22:44

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

 

[_Ivana]

С интересом разглядывая последний график, я могу сказать следующее:
1) можно ввести параметр так называемой "мощности" алгоритма интерполяции - угол наклона прямой линии (или почти прямой) на графике.
2) сравнивать алгоритмы можно по погрешности при малых дискретизациях (стартовой погрешности) и мощности. Например, Лагранж-4 проигрывает сплайнам с производными при малом количестве точек, но за счет большей мощности обгоняет последний при где-то 28 точках и дальше только увеличивает свое преимущество. А Катмулл-Ром при малом количество точек близок к Фарроу (хотя и все равно хуже), но за счет меньшей мощности при увеличении количества точек все больше и больше отстает от него. Интересно поведение Фарроу и сплайна по производным - они имеют одинаковую (насколько можно судить по графикам) мощность, но сплайн по производным всегда на порядок (в 10 раз) точнее
3) из всего перечисленного можно предположить, что максимально достижимая мощность сплайнов с угадыванием производных (типа Катмулла-Рома и т.п.) равна мощности сплайна с известными производными, и следовательно мощности Фарроу (Лагранжа 3 порядка). А значит, как бы я не исхитрялся с алгоритмом угадывания производных, он уступит по мощности точным производным, и следовательно Фарроу Я могу получить выигрыш только на начальном этапе - при малом (относительно) количестве точек, но потом Фарроу обгонит (как Лагранж-4 обогнал точные производные). Значит, теоретические изыскания могут рассчитывать только на то, что этот обгон состоится как можно позже, и математика алгоритма будет при этом проще Фарроу. У меня есть некоторые сомнения в целесообразности и осуществимости таких поисков. Фарроу оказался сильным (мощным ). И он победил меня на поле кубических сплайнов как таковых, хотя сам им в какой-то степени и является. Остается только признать свое поражение...

[Михаил_K]

В данный момент размышления следующие. Фарроу (Лагранж 3-го порядка) - это имхо по сути (хотя многие могут со мной не согласиться) тот же кубический сплайн. Локальный, определяемый для 4-х точек, но действующий только на средний интервал из 3-х. Точно так же, как и сплайны Эрмита, например. Тот факт, что Лагранж при построении проходит через все 4 точки, не дает значения интерполирующей функции в левом и правом интервалах, поскольку при интерполяции собственно самих этих боковых интервалов (когда они являются центральными) полином проходит совершенно не так, как проходит в случае если интервал является правым или левым относительно центрального. Значит нельзя считать условие прохождение полинома через крайние точки самым оптимальным с точки зрения построения сплайна. Я придумал (а Катмулл и Ром до меня) другие условия: сплайн проходит через 2 точки и имеет определенные производные в них же, так чтобы полученная функция имела непрерывную первую производную. Сначала я придумал считать значение этой производной, как тангенс угла-биссектрисы исходных двух (словами трудно объяснить, надо на графике показывать). Только что проверил этот вариант - при малых дискретизациях ошибка выше Катмулла-Рома, при увеличении начиная где-то с 12 точек почти полностью сливается с Катмуллом-Ромом (график ниже). Но само значение производной рассчитывается через тангенсы-арктангенсы и достаточно сложно. Это и была моя вторая идея, которую я проверил только что. А потом я придумал считать производную как тангенс наклона соседних точек - собственно то, что потом узнал как Катмулла-Рома.

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

[thermit]

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

[_Ivana]

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

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

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

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

[ViKo]

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

Предлагаю. Синус с частотой 1kHz оцифруйте с частотой 5kHz. А потом интерполируйте до частоты дискретизации 50kHz. И покажите спектр с полосе 0-25kHz.

[Михаил_K]

и спектр синуса 0.49 от частоты дискретизации ресемпленный в несколько раз выложил. А больше вы вроде ничего не предлагали.

А где? Че-то не увидел

[_Ivana]

Предлагаю. Синус с частотой 1kHz оцифруйте с частотой 5kHz. А потом интерполируйте до частоты дискретизации 50kHz. И покажите спектр с полосе 0-25kHz.

Сделаю. Вечером, как дома буду. По Катмуллу-Рому. Хотя могу для сравнения ещё по Фарроу, Лагранжу 4 порядка, линейному...

А где? Че-то не увидел

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

[Михаил_K]

Понятно. Просто там такой ужОс изображен, что я сразу и не понял. Я так понимаю что наш синус - это пик в районе 4 кГц, и рядышком же стоит его алиас, практически того же уровня. А вот теперь уменьшая частоту синуса, найдите такую, при которой алиас будет ниже дБ на 40 - 45. Ну и всякий мусор тоже должен быть не выше этого уровня. И напишите ее значение относительно частоты дискретизации.

[NiceParty]

Понятно. Просто там такой ужОс изображен, что я сразу и не понял. Я так понимаю что наш синус - это пик в районе 4 кГц, и рядышком же стоит его алиас, практически того же уровня...

Этот "ужОс" вполне предсказуем, т.к. при таком соотношении исходных частот сигнала и дискретизации ни одна полиномиальная интерполяция не даст нормального результата.
Лично по моим наблюдениям для хоть какого-то приемлемого результата при полиномиальной интерполяции максимальная частота исходного сигнала не должна быть больше 1/4 частоты семплирования. А лучше 1/8. Именно поэтому прежде чем применять алгоритмы полиномиальной интерполяции (в частности Лагранж 3-й степени с оптимизацией расчета по Фарроу) рекомендуют апсемплинг в 2, 4 и т.д. раз.
Если Вы знаете вариант полиномиальной интерполяции, который хорошо справляется с сигналами 0,48 частоты дискретизации, то расскажите, я буду весьма признателен.
Заранее спасибо.
С уважением, Александр.

[rudy_b]

...Остается только признать свое поражение...

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

[_Ivana]

Предлагаю. Синус с частотой 1kHz оцифруйте с частотой 5kHz. А потом интерполируйте до частоты дискретизации 50kHz. И покажите спектр с полосе 0-25kHz.

Первый - Катмулл-Ром, второй - Фарроу

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

 

[ViKo]

А линейную интерполяцию можете забомбить, для сравнения?
Что-то много лишнего. Откуда взялись частоты 2 kHz и 3 kHz? Если частота дискретизации была 5 kHz, то будут только 1 kHz и 4 kHz, ну и т.д. 14, 16, 19, 21... После интерполяции в идеале должны были остаться только 1 kHz и 49 kHz.
А много ли точек в вашем кадре, по которому спектр вычисляете? А не увеличить ли количество раз в 250? И окна еще разные можно попробовать.

P.S. А линейный масштаб по частоте приятнее будет смотреться.

[_Ivana]

По порядку - спектр исходного сигнала, после линейной интерполяции ( ) и Лагранж 4-го порядка.

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

 

[ViKo]

Где какой? А красный и синий что за спектры?

[ViKo]

P.S. А линейный масштаб по частоте приятнее будет смотреться.

[_Ivana]

Лагранж-4 из последних 3-х картинок самый зашумленный. Синий - мгновенный, красный - пиковый.

[ViKo]

Ну, и как это понимать? Что линейная интерполяция более качественная???

[_Ivana]

Подозреваю что на вашем примере - да. При ресемплировании того же синуса 1кГц из 44100 в те же 50000 будет совершенно другая картина, спектр линейного преобразования будет содержать больше ВЧ составляющих.

[ViKo]

Подозреваю что на вашем примере - да.

Странно...
То есть, имеем 5 точек на синусоиде. Забили внутри каждого промежутка между точками еще по 9 (увеличили частоту дискретизации в 10 раз). В первый раз тупо ровненько по прямой, во второй - по хитрому полиному Лагранжа. И оказалось, что первый вариант лучше?
Можно узнать ответы на все вопросы в сообщении 91?

[_Ivana]

Количество точек FFT = 16384. А на другие вопросы или ответил, или они риторические.
Впрочем, могу прислать вам вавки - сами проанализируете.

[ViKo]

Количество точек FFT = 16384. А на другие вопросы или ответил, или они риторические.
Впрочем, могу прислать вам вавки - сами проанализируете.

Тогда, если интерес у Вас не пропал, попробуйте интерполировать цифровыми фильтрами. Например, каскадным интегратором-гребенчатым фильтром, CIC. Похоже, интерполяция полиномами не даст качественного звука.
Та же ссылка, что уже давал, а от нее - дальше...
http://ru.wikipedia.org/wiki/%D0%9F%D0%B5%...%86%D0%B8%D1%8F

Ну, и сами синусоиды покажите, что ли? По паре периодов. Может, на глаз что-то неправильное попадется?
А если разрядность увеличить? А если там ограничение возникает?