Интерполяция + (линейная | логарифмическая) шкала + С++
В этой статье я распишу теорию (а также базовые виртуальные классы), в следующей возьмусь за конкретные реализации средствами Qt.
Осторожно: в тексте много графики!
Откуда растут ноги у задачи
В общем, надо мне сделать регулятор холостого хода — такая штука для автомобиля, которая при холостом ходе в зависимости от температуры двигателя должна поддерживать определенные обороты. Поддерживает она их, регулируя заслонку шаговым двигателем.
В общем, мне надо знать текущую температуру. Было решено измерять ее штатными средствами — с терморезистора. Измеряем падение напряжения на нем — получаем сопротивление. Далее из таблицы (поскольку это микроконтроллер) получаем требуемые обороты.
Таблицу эту надо задать (для этого программа пишется средствами Qt). У меня есть несколько точек «сопротивление => температура». Мне надо для каждого кода АЦП (для ряда значений сорпотивлений) получить соответствующую температуру. Поскольку у разных автомобилей эти значения могут быть разными, то надо на экране, сверившись с таблицей, задать несколько точек на кривой.
По ходу дела оказалось, что график этот будет явно в логарифмическом масштабе. Значит, надо его вывести на экран. Как это сделать — читаем дальше.
Постановка задачи
Давайте немного подробнее опишем что нам надо:
Вот такое ТЗ… Ну да ничего, я справился! Давайте и вам помогу.
Да, пока не нырнули — спасибо Equation Editor-у от CodeCogs! С их помощью я лихо построил все математические формулы без всяких Microsoft Equation Editor, которые потом надо еще экспортировать в графику со вставкой сюда. Кстати, там есть и русский редактор. В общем, рекомендую!
Ну и если вместо формул вы видите пустые квадратики — это тоже «спасибо» Equation Editor-у…
Прикрепленный Excel-файл
По ходу написания статьи я все расчеты строил и проверял в таблице Excel с формулами. Оказалось очень удобно. И я решил его выложить для общественного пользования. Там внизу перечислены страницы по разделам. На каждой странице параметры, которые можно менять, отмечены как ячейки с желтым фоном. Остальные клетки лучше не трогать. Впрочем, все формулы можно смело смотреть. Скачивайте файлик и пробуйте на здоровье! Если проблемы с файлом — пишите, вышлю.
Функциональная зависимость
Итак, у нас есть некоторая зависимость — обозначим ее как 
. Здесь у нас — горизонтальная ось графика, — вертикальная. В моем случае было значение сопротивления, — температура.
Почему не ? Ведь вроде бы должно быть так? Так-то оно так, но только в школе в простейшем случае.
— это координаты точки на плоскости. Для простоты определимся использовать Декартову систему координат: 
задает вертикальное смещение горизонтальной оси относительно нуля, задает горизонтальное смещение вертикальной оси относительно нуля.
Все хорошо тогда, когда мы рисуем на бумаге эту самую систему координат и в ней ставим точки. Там и вправду — выбрали центр, линейкой отложили сюда, потом туда. А вот при построении графика в какой-то программе уже тонкости начинаются — что считать нулем? Что считать за «+», а что за «-«? Я рисую для этой статьи графику в CorelDRAW — там центр считается снизу слева (его можно передвинуть куда надо).
Да и в каких единицах график-то? В сантиметрах? А почему? У меня следующий этап будет реализация на С++ средствами Qt, так там я сделаю окно QWidget, у которого по умолчанию ноль — это слева сверху; единицы измерения — экранные пиксели.
Ну и не забываем о том, что это все эти красивые рассуждения справедливы пока что для линейной шкалы, а у нас маячит за горизонтом логарифмическая. Там вообще черт знает что будет!
Но это только лишь точка. А у нас будет какая-то линия, точнее — много линий. Что там будут за преобразования?
Вот именно поэтому мы с самого начала должны четко разделить функциональную зависимость и преобразования координат.
Итак, давайте договоримся о следующем: у нас есть некоторый абстрактный процесс, который описывается функциональной зависимостью . При отображении на экран используется преобразование в координаты , где , . Следующие шаги — это прояснить эти самые и .
Но отложим пока в сторону координаты — нам надо как-то задать нашу функцию (помните ТЗ)? Причем задать в тех самых абстрактных координатах . Этим и займемся.
Интерполяция
В моем случае был известен ряд точек :
| , Ω | , ˚ |
|---|---|
| 180 | 100 |
| 6 000 | 0 |
| 30 000 | -30 |
Методов интерполяции много, все рассматривать я тут не буду. Лично мне приглянулся вначале интерполяционный многочлен Лагранжа. Он весьма прост в расчете и реализации, а также в настройке. Там предполагается, что задано множество из точек вида (тут мы на время таки вернемся к заданию точек в виде — так уж принято в математике).
Многочлен вычисляется как , где .
Математика испугала? Хм… Ладно, напишу на языке С++:
Как видите, все достаточно тривиально (насколько тривиальными могут быть полиномы).
Еще одно большое достоинство полиномов Лагранжа — их легко можно промоделировать в таблице Excel-я, что я и делал.
Потом, правда, все стало немного печально, т. к. у этих полиномов, как и у любых других, на графике видны вибрации. Т. е. они не могут дать прямые линии — постоянные значения. В моем случае я не смог их настроить дОлжным образом — они выгибались в явно недопустимые числа. Поэтому мне пришлось от них отказаться…
Работая в Corel, я был близко знаком с кривыми Безье — тоже достаточно удобное и простое представление табличных данных. Весьма легко реализуется в программировании. Однако это уже не интерполяция, а, скорее, аппроксимация, т. к. тут приходится подгонять кривую к нужному виду.
В итоге, внимательно присмотревшись к своей функции, я понял, что у меня вполне прокатит кусочно-линейная интерполяция — прямые отрезки между заданными линиями. Не то, чтобы совсем уж фен-шуйно, но зато легко реализуемо и удобно настраиваемо.
Говоря языком математики, мы между точками и проводим прямые линии вида .
Опять же, на языке С++ это будет выглядеть так:
Тоже ничего революционного, не так ли?
Есть одно существенное различие между полиномом Лагранжа и линейной интерполяцией: у первого нельзя явно задать значения за пределами точек — они вычисляются, у второго можно это дело контролировать. Также и поэтому я в конечном итоге остановился на линейном варианте. Более того — в логарифмическом масштабе, к которому я стремился, линейные отрезки дают более подходящий мне вариант.
Впрочем, не будем заморачиваться сейчас на методах интерполяции. Давайте лучше мы сделаем базовый класс, от которого будем наследовать реализации различных методов $#*@!поляции.
Базовый класс для задания/расчета функции
Что этот класс должен уметь делать? Мне кажется, что такой класс должен:
Еще есть мысли? Если будут — пишите в комментариях, добавим!
Получается такой вот класс:
(Тем, кто недоволен моим стилем и структурой — предложите объективно лучше!)
(Тем, кто найдет ошибки в коде — спасибо!)
Думаю, тут все очевидно.
Для координат используется представление точки в виде QPointF (пара чисел в виде qreal, qreal. «На всех платформах, кроме ARM, используется double» — так написано для Qt 4.8).
В следующей статье мы распишем пару вариантов реализации этого класса.
Функция преобразования для вертикальной/горизонтальной шкалы
Есть линейные и логарифмические шкалы. Учитывая, что вертикальная шкала может быть сделана в одном формате, а горизонтальная — в другом, мы получаем четыре варианта графика:
Вариант первый — обе шкалы линейные. Вариант второй — обе логарифмические. Варианты третий и четвертый — смешанные графики. Кстати, в моем случае именно смешанный случай в итоге и подошел, т. к. по горизонтали у меня потребовался логарифмический масштаб, по вертикали — линейный.
Следовательно, задачу отображения нужно решать отдельно для обеих осей.
Напомним, что при отображении на экран используется преобразование в координаты , где , . Наша дальнейшая задача — построить эти функции для линейного и логарифмического случаев.
Что это за функции такие? На вход они получают координату в абстрактных (для компьютерной подпрограммы отображения на экран) координатах, на выход дают в экранных («экранные» координаты будут для разных операционных систем разными»). Для расчета им нужно знать следующее:
Базовый класс для преобразований шкал
Давайте сформулируем желаемую функциональность виртуального класса преобразований для шкалы, от которого будут унаследованы реализации шкал:
Реализация может выглядеть так:
Линейное преобразование
Давайте отдельно рассмотрим линейное преобразование для горизонтальной и вертикальной оси.
Расчет этих констант достаточно прост — это решение системы двух уравнений:
Еще важно уметь делать обратное преобразование — скажем, координаты указателя мыши перевести в абстрактные координаты. Также ничего сложного:
Шаг в данном случае для расчета не используется, но он потом нам пригодится в реализации на С++ для расчета смещения.
Как это будет использоваться на практике? Да все просто! Горизонтальное преобразование: — граница картинки графика, соответствующая (как правило, слева), — (как правило, справа), — шаг вывода картинки по горизонтали. Вертикальное преобразование — аналогично, но по вертикали (у нас в Qt будет нижней границей картинки, — верхней, причем scr_
Логарифмическое преобразование
А теперь окунемся туда, ради чего все это закрутилось:
(на графике не логарифм нарисован, а что-то похожее на него. Сделано это специально, т. к. логарифм тут будет не очень нагляден)
Вроде бы базовую математику рассмотрели. Нашли ошибки или неточности — пишите в комментариях, буду благодарен!
Со временем напишу следующую статью — реализацию этой математики средствами Qt языка C++.
Интерполяция данных: соединяем точки так, чтобы было красиво
Как построить график по n точкам? Самое простое — отметить их маркерами на координатной сетке. Однако для наглядности их хочется соединить, чтобы получить легко читаемую линию. Соединять точки проще всего отрезками прямых. Но график-ломаная читается довольно тяжело: взгляд цепляется за углы, а не скользит вдоль линии. Да и выглядят изломы не очень красиво. Получается, что кроме ломаных нужно уметь строить и кривые. Однако тут нужно быть осторожным, чтобы не получилось вот такого:
Немного матчасти
Восстановление промежуточных значений функции, которая в данном случае задана таблично в виде точек P1 .  Pn, называется интерполяцией. Есть множество способов интерполяции, но все они могут быть сведены к тому, что надо найти n – 1 функцию для расчёта промежуточных точек на соответствующих сегментах. При этом заданные точки обязательно должны быть вычислимы через соответствующие функции. На основе этого и может быть построен график:
Функции fi могут быть самыми разными, но чаще всего используют полиномы некоторой степени. В этом случае итоговая интерполирующая функция (кусочно заданная на промежутках, ограниченных точками Pi) называется сплайном.
В разных инструментах для построения графиков — редакторах и библиотеках — задача «красивой интерполяции» решена по-разному. В конце статьи будет небольшой обзор существующих вариантов. Почему в конце? Чтобы после ряда приведённых выкладок и размышлений можно было поугадывать, кто из «серьёзных ребят» какие методы использует.
Ставим опыты
Самый простой пример — линейная интерполяция, в которой используются полиномы первой степени, а в итоге получается ломаная, соединяющая заданные точки.
Давайте добавим немного конкретики. Вот набор точек (взяты почти с потолка):
Результат линейной интерполяции этих точек выглядит так:
Однако, как отмечалось выше, иногда хочется получить в итоге гладкую кривую.
Что есть гладкость? Бытовой ответ: отсутствие острых углов. Математический: непрерывность производных. При этом в математике гладкость имеет порядок, равный номеру последней непрерывной производной, и область, на которой эта непрерывность сохраняется. То есть, если функция имеет гладкость порядка 1 на отрезке [a; b], это означает, что на [a; b] она имеет непрерывную первую производную, а вот вторая производная уже терпит разрыв в каких-то точках.
У сплайна в контексте гладкости есть понятие дефекта. Дефект сплайна — это разность между его степенью и его гладкостью. Степень сплайна — это максимальная степень использованных в нём полиномов.
Важно отметить, что «опасными» точками у сплайна (в которых может нарушиться гладкость) являются как раз Pi, то есть точки сочленения сегментов, в которых происходит переход от одного полинома к другому. Все остальные точки «безопасны», ведь у полинома на области его определения нет проблем с непрерывностью производных.
Чтобы добиться гладкой интерполяции, нужно повысить степень полиномов и подобрать их коэффициенты так, чтобы в граничных точках сохранялась непрерывность производных.
Традиционно для решения такой задачи используют полиномы третьей степени и добиваются непрерывности первой и второй производной. То, что получается, называют кубическим сплайном дефекта 1. Вот как он выглядит для наших данных:
Кривая, действительно, гладкая. Но если предположить, что это график некоторого процесса или явления, который нужно показать заинтересованному лицу, то такой метод, скорее всего, не подходит. Проблема в ложных экстремумах. Появились они из-за слишком сильного искривления, которое было призвано обеспечить гладкость интерполяционной функции. Но зрителю такое поведение совсем не кстати, ведь он оказывается обманут относительно пиковых значений функции. А ради наглядной визуализации этих значений, собственно, всё и затевалось.
Так что надо искать другие решения.
Другое традиционное решение, кроме кубических сплайнов дефекта 1 — полиномы Лагранжа. Это полиномы степени n – 1, принимающие заданные значения в заданных точках. То есть членения на сегменты здесь не происходит, вся последовательность описывается одним полиномом.
Но вот что получается:
Гладкость, конечно, присутствует, но наглядность пострадала так сильно, что… пожалуй, стоит поискать другие методы. На некоторых наборах данных результат выходит нормальный, но в общем случае ошибка относительно линейной интерполяции (и, соответственно, ложные экстремумы) может получаться слишком большой — из-за того, что тут всего один полином на все сегменты.
В компьютерной графике очень широко применяются кривые Безье, представленные полиномами k-й степени.
Они не являются интерполирующими, так как из k + 1 точек, участвующих в построении, итоговая кривая проходит лишь через первую и последнюю. Остальные k – 1 точек играют роль своего рода «гравитационных центров», притягивающих к себе кривую.
Вот пример кубической кривой Безье:
Как это можно использовать для интерполяции? На основе этих кривых тоже можно построить сплайн. То есть на каждом сегменте сплайна будет своя кривая Безье k-й степени (кстати, k = 1 даёт линейную интерполяцию). И вопрос только в том, какое k взять и как найти k – 1 промежуточную точку.
Здесь бесконечно много вариантов (поскольку k ничем не ограничено), однако мы рассмотрим классический: k = 3.
Чтобы итоговая кривая была гладкой, нужно добиться дефекта 1 для составляемого сплайна, то есть сохранения непрерывности первой и второй производных в точках сочленения сегментов (Pi), как это делается в классическом варианте кубического сплайна.
Решение этой задачи подробно (с исходным кодом) рассмотрено здесь.
Вот что получится на нашем тестовом наборе:
Стало лучше: ложные экстремумы всё ещё есть, но хотя бы не так сильно отличаются от реальных.
Думаем и экспериментируем
Можно попробовать ослабить условие гладкости: потребовать дефект 2, а не 1, то есть сохранить непрерывность одной только первой производной.
Достаточное условие достижения дефекта 2 в том, что промежуточные контрольные точки кубической кривой Безье, смежные с заданной точкой интерполируемой последовательности, лежат с этой точкой на одной прямой и на одинаковом расстоянии:
Методом проб и ошибок эвристика для расчёта расстояния от точки интерполируемой последовательности до промежуточной контрольной получилась такой:
Первая и последняя промежуточные контрольные точки равны первой и последней точке графика соответственно (точки C1 (1) и Cn – 1 (2) совпадают с точками P1 и Pn соответственно).
В этом случае получается вот такая кривая:
Как видно, ложных экстремумов уже нет. Однако если сравнивать с линейной интерполяцией, местами ошибка очень большая. Можно сделать её ещё меньше, но тут в ход пойдут ещё более хитрые эвристики.
Эвристика для вычисления расстояний будет такой:
Результат получается такой:
В результате на шестом сегменте ошибка уменьшилась, а на седьмом — увеличилась: кривизна у Безье на нём оказалась больше, чем хотелось бы. Исправить ситуацию можно, принудительно уменьшив кривизну и тем самым «прижав» Безье ближе к отрезку прямой, которая соединяет граничные точки сегмента. Для этого используется следующая эвристика:
Результат следующий:
На этом было принято решение признать цель достигнутой.
Может быть, кому-то пригодится код.
А как люди-то делают?
Обещанный обзор. Конечно, перед решением задачи мы посмотрели, кто чем может похвастаться, а уже потом начали разбираться, как сделать самим и по возможности лучше. Но вот как только сделали, не без удовольствия ещё раз прошлись по доступным инструментам и сравнили их результаты с плодами наших экспериментов. Итак, поехали.
MS Excel
Это очень похоже на рассмотренный выше сплайн дефекта 1, основанный на кривых Безье. Правда, в отличие от него в чистом виде, тут всего два ложных экстремума — первый и второй сегменты (у нас было четыре). Видимо, к классическому поиску промежуточных контрольных точек тут добавляются ещё какие-то эвристики. Но ото всех ложных экстремумов они не спасли.
LibreOffice Calc
В настройках это названо кубическим сплайном. Очевидно, он тоже основан на Безье, и вот тут уже точная копия нашего результата: все четыре ложных экстремума на месте.
Есть там ещё один тип интерполяции, который мы тут не рассматривали: B-сплайн. Но для нашей задачи он явно не подходит, так как даёт вот такой результат 
Highcharts, одна из самых популярных JS-библиотек для построения диаграмм
Тут налицо «метод касательных» в варианте равенства расстояний от точки интерполируемой последовательности до промежуточных контрольных. Ложных экстремумов нет, зато есть сравнительно большая ошибка относительно линейной интерполяции (седьмой сегмент).
amCharts, ещё одна популярная JS-библиотека
Картина очень похожа на экселевскую, те же два ложных экстремума в тех же местах.
Coreplot, самая популярная библиотека построения графиков для iOS и OS X
Есть ложные экстремумы и видно, что используется сплайн дефекта 1 на основе Безье.
Библиотека открытая, так что можно посмотреть в код и убедиться в этом.
aChartEngine, вроде как самая популярная библиотека построения графиков для Android
Больше всего похоже на кривую Безье степени n – 1, хотя в самой библиотеке график называется «cubic line». Странно! Как бы то ни было, тут не только присутствуют ложные экстремумы, но и в принципе не выполняются условия интерполяции.