Емельянов Эдуард Владимирович (eddy_em) wrote,
Емельянов Эдуард Владимирович
eddy_em

И что я сразу о GSL не вспомнил?

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

И наткнулся я на эту идею в этой статье:
@ARTICLE{2010ApOpt..49.6489M,
   author = {{Mahajan}, V.~N. and {Aftab}, M.},
    title = "{Systematic comparison of the use of annular and Zernike circle polynomials for annular wavefronts}",
  journal = {\ao},
     year = 2010,
    month = nov,
   volume = 49,
    pages = {6489},
      doi = {10.1364/AO.49.006489},
   adsurl = {http://adsabs.harvard.edu/abs/2010ApOpt..49.6489M},
  adsnote = {Provided by the SAO/NASA Astrophysics Data System}
}

Вот простенький тестовый код:
#include <stdio.h>
#include <gsl/gsl_linalg.h>
#include <math.h>

void fillmatrix(gsl_matrix *m, double *x){
	int i, j;
	size_t W = m->size2, H = m->size1, T = m->tda;
	double *xptr = x;
	for(i = 0; i < H; i++, xptr++){
		double *ptr = &(m->data[i*T]);
		for(j = 0; j < W; j++, ptr++)
			*ptr = sin((double)(j+1) * (*xptr));
	}
}

int main(int argc, char **argv){
	double b_data[] = {0.37, 1.07, 0.95, 0.43, -0.27};
	double xx[] = {0.1, 0.34, 0.76, 0.93, 1.12};
	gsl_matrix *M = gsl_matrix_calloc(5, 4);
	fillmatrix(M, xx);
	gsl_vector *x = gsl_vector_calloc(4);
	gsl_vector_view b = gsl_vector_view_array (b_data, 5);
	gsl_vector *tau = gsl_vector_calloc(4); // min(size(M))
	gsl_linalg_QR_decomp(M, tau);
	gsl_vector *residual = gsl_vector_calloc(5);
	gsl_linalg_QR_lssolve(M, tau, &b.vector, x, residual);
	printf ("x = \n");
	gsl_vector_fprintf(stdout, x, "%g");
	printf ("\nresidual = \n");
	gsl_vector_fprintf(stdout, residual, "%g");
	gsl_matrix_free(M);
	gsl_vector_free(x);
	gsl_vector_free(tau);
	gsl_vector_free(residual);
	return 0;
}

Здесь я просто раскладываю функцию 1.25·sin(3·x) (округленную до двух знаков после запятой) на неравномерной сетке аргумента x по синусному базису (с первой по четвертую гармоники).
Результат получается вполне приличным:
gcc -lgsl -lgslcblas -lm G.c -o vect && ./vect 
x = 
0.130783
-0.143849
1.32437
-0.0124236

residual = 
-0.00101798
0.000479872
-0.00021396
0.00014092
-2.66255e-05

(если результаты округлять до четвертого знака, то получается более точно: 1.27 вместо 1.32). Понятно, что при увеличении количества точек раз в 100 погрешность раз в 10 уменьшится.
Эксперименты с разложением волнового фронта по полиномам Цернике на непрямоугольной сетке (точнее — для реальных координат отверстий в диафрагме Гартманна) оставлю на следующую неделю, т.к. после обеда надо помыть машину и ехать в командировку в Ставрополь, слушать защиту наших магистрантов (завтра с утра). Так что в понедельник-вторник выясню, подойдет ли этот метод для разложения волнового фронта на кольце (да еще и на неравномерной сетке), а также — для разложения градиента этого самого фронта. Если все будет ОК, можно будет приступить к анализу реальных данных (во всяком случае, по моделированию статью можно будет сразу написать, а применительно к реальным данным — чуть позже; все равно пока нет у нас ни одной нормальной пары гартманнограмм, все время с погодой не везло).
Tags: c, gsl, аппроксимация, велосипедостроение
Subscribe

  • DHT22/DHT11 на STM32F103

    Добил шайтана! Сначала ожидал, что нужно будет полноценным захватом ШИМ пользоваться, но т.к. в протоколе неинформативная часть имеет постоянную…

  • Свеженькие железячки

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

  • BMP180 на STM32F103

    Добавил еще один сниппет — работа с BMP180 (датчик температуры и давления). Опять в даташите формула для вычисления "компенсированных" значений…

promo eddy_em август 17, 2019 12:33 3
Buy for 10 tokens
Юра намедни напечатал корпус для хронометра. Для первого блина получилось неплохо: И еще немного фотографий:
  • Post a new comment

    Error

    Anonymous comments are disabled in this journal

    default userpic

    Your reply will be screened

  • 4 comments