Ещё по поводу генерации синусов

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

И пришлось написать табличную реализацию sincosf. Под катом код и возможно описание.




#define KPI		3.14159265f

static void
ksincosf(float a, float x[2])
{
	static float	tab[][4] = {

		{-4.06662660e-05f, -1.58916156e-08f, 6.25000000e-02f, 0.00000000e+00f},
		{-4.05074651e-05f, -1.22046410e-04f, 6.23779694e-02f, 6.24593178e-02f},
		{-4.01904834e-05f, -2.43600341e-04f, 6.20123542e-02f, 1.24674733e-01f},
		{-3.97165587e-05f, -3.64203017e-04f, 6.14045820e-02f, 1.86403296e-01f},
		{-3.90875417e-05f, -4.83383489e-04f, 6.05570263e-02f, 2.47403959e-01f},
		{-3.83058888e-05f, -6.00676358e-04f, 5.94729967e-02f, 3.07438514e-01f},
		{-3.73746521e-05f, -7.15623599e-04f, 5.81567263e-02f, 3.66272529e-01f},
		{-3.62974682e-05f, -8.27776345e-04f, 5.66133552e-02f, 4.23676257e-01f},
		{-3.50785435e-05f, -9.36696642e-04f, 5.48489101e-02f, 4.79425538e-01f},
		{-3.37226378e-05f, -1.04195915e-03f, 5.28702812e-02f, 5.33302673e-01f},
		{-3.22350460e-05f, -1.14315284e-03f, 5.06851949e-02f, 5.85097272e-01f},
		{-3.06215769e-05f, -1.23988254e-03f, 4.83021841e-02f, 6.34607080e-01f},
		{-2.88885313e-05f, -1.33177053e-03f, 4.57305543e-02f, 6.81638760e-01f},
		{-2.70426766e-05f, -1.41845798e-03f, 4.29803476e-02f, 7.26008655e-01f},
		{-2.50912208e-05f, -1.49960638e-03f, 4.00623036e-02f, 7.67543502e-01f},
		{-2.30417843e-05f, -1.57489886e-03f, 3.69878171e-02f, 8.06081108e-01f},
		{-2.09023702e-05f, -1.64404138e-03f, 3.37688941e-02f, 8.41470984e-01f},
		{-1.86813327e-05f, -1.70676396e-03f, 3.04181042e-02f, 8.73574935e-01f},
		{-1.63873450e-05f, -1.76282166e-03f, 2.69485322e-02f, 9.02267594e-01f},
		{-1.40293651e-05f, -1.81199559e-03f, 2.33737269e-02f, 9.27436917e-01f},
		{-1.16166009e-05f, -1.85409370e-03f, 1.97076476e-02f, 9.48984619e-01f},
		{-9.15847408e-06f, -1.88895162e-03f, 1.59646104e-02f, 9.66826556e-01f},
		{-6.66458359e-06f, -1.91643323e-03f, 1.21592317e-02f, 9.80893057e-01f},
		{-4.14466805e-06f, -1.93643120e-03f, 8.30637152e-03f, 9.91129190e-01f},
		{-1.60856766e-06f, -1.94886746e-03f, 4.42107510e-03f, 9.97494986e-01f},
		{ 9.33814142e-07f, -1.95369342e-03f, 5.18514476e-04f, 9.99965585e-01f},
	};

	float		*l, r, u;
	int		j, m = 0;

	if (a < 0.0) {

		m += 1;
		a = -a;
	}

	if (a > (KPI / 2.0f)) {

		m += 2;
		a = KPI - a;
	}

	r = a * 16.0f;
	j = (int) r;
	r = r - (float) j;
	l = tab[j];

	u = l[1] + l[0] * r;
	u = l[2] + u * r;
	u = l[3] + u * r;
	x[0] = (m & 1) ? -u : u;

	a = (KPI / 2.0f) - a;

	r = a * 16.0f;
	j = (int) r;
	r = r - (float) j;
	l = tab[j];

	u = l[1] + l[0] * r;
	u = l[2] + u * r;
	u = l[3] + u * r;
	x[1] = (m & 2) ? -u : u;
}


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

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

Ошибки приближаются к точности float, на картинке график разности моего синуса и sin откуда-то из glibc (libm).



Таблица гененриуется вот таким скриптом.


#!/usr/bin/octave -q

T = 0.125;

t = 0:T:1.65
ts = sin(t);
tc = cos(t) * T;
tab = [];

for i=1:13

	bet = 6 * (ts(i + 1) - ts(i)) - 3 * (tc(i + 1) + tc(i));
	p(1) = - bet / 3;
	p(2) = (tc(i + 1) - tc(i) + bet) / 2;
	p(3) = tc(i);
	p(4) = ts(i);

	tab = [tab; p];
end


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

Можно лучше? быстрее? в комменты.
  • +4
  • 04 ноября 2013, 21:41
  • amaora

Комментарии (8)

RSS свернуть / развернуть
осциллятор второго порядка вполне себе устойчивый
web.archive.org/web/20120308235958/http://www.ied.com/~petr/Oscillators.html
0
  • avatar
  • _pv
  • 04 ноября 2013, 23:01
С усточивостью проблем нет, частоту нельзя задать точно, ошибка ~1% для метода второго порядка. А боольший порядок уже слишком тяжел получается, лучше интегрировать сам угол без заметных ошибок а когда надо делать sincos.
0
какой еще больший порядок?
Y[i] = K*Y[i-1] — Y[i-2], где К = 2*cos(w)
косинус считается только один раз и определят частоту осциллятора с какой надо точностью, откуда 1%?
0
Частота переменная, это же вращение двигателя.
0
ммм как сладко=)) какое интересное решение=))
+1
полиномом от 0 до pi/2, 7 умножений, без таблиц разница те же 1e-7
cos(x) = sin(pi/2-x)

float sinf(float x){
  return x * (1.0000011907980209f + x * (-0.0000173804054785742f + x * (-0.16657758167275724f + x * (-0.00021883889593866195f + x * (0.008618152688037924f + x * (-0.00019515619194458684f + x * (-0.00013938986552924126f)))))));
}
+1
  • avatar
  • _pv
  • 05 ноября 2013, 02:45
А откуда полином, МНК по 7-ми точкам?

У меня всегда было ощущение, что точность будет падать с ростом степени полинома, даже проверять не стал насколько.
0
точек несколько тысяч задал, хотя и небольшому количеству, наверное, тоже самое получится.
можно еще конечно по Тэйлору sin(x) = x^(2n+1) * (-1)^n / (2n+1)!, но он ошибку при том же небольшом количестве степеней больше даст.
0
Только зарегистрированные и авторизованные пользователи могут оставлять комментарии.