diff --git a/apps/plugins/test_grey.c b/apps/plugins/test_grey.c
index 31f6310..1fe8759 100644
--- a/apps/plugins/test_grey.c
+++ b/apps/plugins/test_grey.c
@@ -100,6 +100,16 @@ PLUGIN_HEADER
 
 #define STEPS 16
 
+/* scale factor for fixed version of natural b-spline interpolation */
+#define FACTOR 2048
+
+/* b-spline interpolation helper macros 
+ * used in spline_z_param() function
+ */
+#define h(i) (x[i+1] - x[i])
+#define b(i) (((y[i+1] - y[i])*FACTOR)/(x[i+1] - x[i]))
+
+
 GREY_INFO_STRUCT
 
 static const unsigned char dither_matrix[16][16] = {
@@ -124,6 +134,11 @@ static const unsigned char dither_matrix[16][16] = {
 static unsigned char input_levels[STEPS+1];
 static unsigned char lcd_levels[STEPS+1];
 
+/* b-spline interpolation helper arrays */
+static int u[STEPS];
+static int v[STEPS];
+static int z[STEPS+1];
+
 static unsigned char *gbuf;
 static size_t gbuf_size = 0;
 
@@ -153,6 +168,38 @@ static void fill_rastered(int bx, int by, int bw, int bh, int step)
     }
 }
 
+/* b-spline interpolation function */
+static inline int spline(int ai,int bi,int ci,int di, int dx)
+{
+    return (ai*FACTOR + (dx)*(bi+(dx)*(ci+(dx)*di)))/FACTOR;
+}
+
+static void spline_z_param(unsigned char *x, unsigned char *y, 
+                            int *u, int *v, int *z, unsigned char n)
+{
+    /* Calculate z coefficient of natural cubic b-spline
+     * with n+1 joint points with coordinates (x,y)
+     */
+
+    int i;
+
+    /* Gaussian elimination */
+    u[1] = 2*FACTOR*(h(0)+h(1));
+    v[1] = 6*(b(1)-b(0));
+    for (i=2; i<=n-1;i++)
+    {
+        u[i] = 2*FACTOR*(h(i-1)+h(i))-(((h(i-1)*h(i-1))*FACTOR*FACTOR)/u[i-1]);
+        v[i] = 6*(b(i)-b(i-1))-(((h(i-1)*v[i-1])*FACTOR)/u[i-1]);
+    }
+
+    /* Back substitution */
+    z[n] = 0;
+    for (i=n-1;i>=1;i--)
+        z[i] = ((v[i]-h(i)*z[i+1])*FACTOR)/u[i];
+
+    z[0] = 0;
+}
+                          
 /* plugin entry point */
 enum plugin_status plugin_start(const void* parameter)
 {
@@ -161,6 +208,9 @@ enum plugin_status plugin_start(const void* parameter)
     int button, i, j, l, fd;
     unsigned char filename[MAX_PATH];
 
+    /* interpolation parameters */
+    int bi,ci,di;
+
     /* standard stuff */
     (void)parameter;
     
@@ -249,15 +299,24 @@ enum plugin_status plugin_start(const void* parameter)
                          rb->fdprintf(fd, "%3d: %3d\n", input_levels[i],
                                       lcd_levels[i]);
                     rb->fdprintf(fd,"\n\nInterpolated lcdlinear matrix\n");
+                    rb->fdprintf(fd,"with natural cubic b-spline\n");
+
+                    /* calculate spline z coefficient */
+                    spline_z_param(input_levels,lcd_levels,u,v,z,STEPS);
 
-                    for (i = 1; i <= STEPS; i++)
+                    for (i = 0; i < STEPS; i++)
                     {
+                        bi = -((input_levels[i+1]-input_levels[i])*z[i+1])/6 - 
+                             (input_levels[i+1]-input_levels[i])*z[i]/3 + 
+                             ((lcd_levels[i+1]-lcd_levels[i])*FACTOR)/
+                             (input_levels[i+1]-input_levels[i]);
+                        ci = z[i]/2;
+                        di = (z[i+1] - z[i])/(6*(input_levels[i+1]-input_levels[i]));
+
                         for (j=0; j < STEPS; j++)
-                        {
-                            rb->fdprintf(fd, "%3d, ", 
-                                lcd_levels[i-1] + 
-                                ((lcd_levels[i] - lcd_levels[i-1])*j)/STEPS);
-                        }
+                            rb->fdprintf(fd,"%3d, ", 
+                                         spline(lcd_levels[i],bi,ci,di,j));
+
                         rb->fdprintf(fd, "\n");
                     }
                     rb->close(fd);