GSL庫之三次樣條插值的使用

背景

因爲工做緣由，對數值計算，fft還有小波使用頻率比較高；相關可選的方案較多，好比Matlab，簡易而又強大，然而用久了，看不到其源碼就成爲個人心結，順着貪婪的本性，找到了GSL，以三次樣條插值爲例，記錄其用法。

環境

GSL庫的安裝

sudo apt-get install libgsl0-dev libgsl0ldbl

編輯環境：emacs

編譯環境：gcc + make 

結果顯示：graph + evince

graph工具須要安裝plotutils包：

sudo apt-get install plotutils

evince爲ubuntu自帶的文檔查看器

固然，你也能夠在windows或Mac上使用。

場景

對單片機採樣與實際電流值創建對應關係表，以理想狀況論，該表應爲線性關係（感應磁路+AD轉換電路），但實驗條件下，磁路飽和與元器件差別問題使該表並不能保持較好的線性關係，綜合多方面因素，以三次樣條插值來構建其對應關係比較合理。

實現

GSL主頁    GSL庫Interpolation模塊文檔    三次樣條插值原理    （維基又被黑了？MLGB）

圖1 實驗數據

源碼

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_spline.h>

#define DATA_COUNTS 19
#define MAX_LINES 1024

int main()
{
   //double curr[DATA_COUNTS] = {1,3,6,8,10,20,30,50,80,100,200,250,300,350,400,450,500,600,630,650};
   //double ad[DATA_COUNTS] = {1,12,18,21,25,42,60,97,154,192,384,463,544,589,644,700,736,820,832,856};
	double curr[DATA_COUNTS] = {0,3,6,8,10,20,30,50,80,100,200,250,300,400,450,500,600,630,650};
	double ad[DATA_COUNTS] = {1,2,7,11,15,30,43,72,121,155,334,419,500,650,716,778,893,924,946};
	double xi = 0;
	double yi = 0;

	gsl_interp_accel *acc = gsl_interp_accel_alloc ();
	gsl_spline *spline = gsl_spline_alloc (gsl_interp_cspline, DATA_COUNTS);
	gsl_spline_init (spline, ad, curr, DATA_COUNTS);	

	for (xi=1; xi<=ad[DATA_COUNTS-1]; xi++){
		yi = gsl_spline_eval (spline, xi, acc);
		printf("%g %g\n", xi, yi);
	}
	
	gsl_spline_free (spline);
	gsl_interp_accel_free (acc);
		
	return 0;
}

該源碼參照 gsl參考手冊——三次樣條插值實例程序 。

Makefile文件

my_interp : my_interp.c
	gcc -Wall -g $< -L/usr/lib -lgsl -lgslcblas -lm -o $@

clean:
	rm -rf my_interp *~ *ps *dat

run:
	./my_interp>interp.dat
	graph -T ps<interp.dat>interp.ps

show:
	evince interp.ps

其中，run命令下的 「./my_interp>interp.dat」將插值後的數據重定向到文件中，便於繪圖工具處理；
「graph -T ps<interp.dat >interp.ps」將intrp.dat文件的數據輸入給PostScript程序並造成.ps文件；
show命令下的"evince interp.ps"使用ubuntu默認的文檔察看器顯示圖形。

可視化

Makefile文件命令依次執行以下：

$ make
gcc -Wall -g my_interp.c -L/usr/lib -lgsl -lgslcblas -lm -o my_interp

$ make run
./my_interp>interp.dat
graph -T ps<interp.dat>interp.ps

$ make show
evince interp.ps

圖形效果以下：

圖2 插值數據可視化

固然也能夠使用gnuplot可視化數據：

gnuplot> set term x11
Terminal type set to 'x11'
Options are ' nopersist'

gnuplot> set datafile separator " "

gnuplot> plot 'interp.dat' every::1::946 using 1:2 with lines

圖示以下：

圖3 插值數據可視化——gnuplot

程序源碼中註釋的數據，是更換了另外一個磁路結構後測量到的，相應地其圖形以下：

圖4 另外一個磁路結構實驗數據可視化

從產品研發的角度來看，圖2和圖3對應的磁路結構比圖4對應的磁路結構好得多，其優劣高下立判。