File: gradsearch_datafit_example.c

package info (click to toggle)
liquid-dsp 1.7.0-1
  • links: PTS, VCS
  • area: main
  • in suites: forky, sid, trixie
  • size: 9,216 kB
  • sloc: ansic: 115,859; sh: 3,513; makefile: 1,350; python: 274; asm: 11
file content (149 lines) | stat: -rw-r--r-- 3,673 bytes parent folder | download | duplicates (2)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
//
// gradsearch_datafit_example.c
//
// Fit 3-parameter curve to sampled data set in the minimum
// mean-squared error sense.
//

#include <stdio.h>
#include <stdlib.h>
#include <math.h>

#include "liquid.h"

#define OUTPUT_FILENAME "gradsearch_datafit_example.m"

// gradient search data set
struct gsdataset {
    float * x;
    float * y;
    unsigned int n;
};

// gradient search curve-fit error
float gserror(void * _dataset,
              float * _v,
              unsigned int _n);

// parameterized function
float gsfunc(float _x, float * _v)
{
    float c0 = _v[0];
    float c1 = _v[1];
    float c2 = _v[2];

    return c0 + sincf(c1*(_x-c2));
}


int main() {
    // options
    unsigned int num_samples = 400;     // number of samples
    float sig = 0.1f;                   // noise variance
    unsigned int num_iterations = 1000; // number of iterations to run

    float v[3] = {1, 1, 1};
    unsigned int i;

    // range
    float xmin = 0.0f;
    float xmax = 6.0f;
    float dx = (xmax - xmin) / (num_samples-1);

    // generate data set
    float x[num_samples];
    float y[num_samples];
    for (i=0; i<num_samples; i++) {
        x[i] = xmin + i*dx;
        y[i] = sincf(x[i]) + randnf()*sig;
    }
    struct gsdataset q = {x, y, num_samples};

    // create gradsearch object
    gradsearch gs = gradsearch_create((void*)&q, v, 3, gserror, LIQUID_OPTIM_MINIMIZE);

    float rmse;

    // execute search
    //rmse = gradsearch_run(gs, num_iterations, -1e-6f);

     // open output file
    FILE*fid = fopen(OUTPUT_FILENAME,"w");
    fprintf(fid,"%% %s : auto-generated file\n", OUTPUT_FILENAME);
    fprintf(fid,"clear all;\n");
    fprintf(fid,"close all;\n");

    // execute search one iteration at a time
    fprintf(fid,"u = zeros(1,%u);\n", num_iterations);
    for (i=0; i<num_iterations; i++) {
        rmse = gserror((void*)&q,v,3);
        fprintf(fid,"u(%3u) = %12.4e;\n", i+1, rmse);

        gradsearch_step(gs);

        if (((i+1)%100)==0)
            gradsearch_print(gs);
    }

    // print results
    printf("\n");
    gradsearch_print(gs);
    printf("  c0 = %12.8f, opt = 1\n", v[0]);
    printf("  c1 = %12.8f, opt = 0\n", v[1]);
    printf("  c2 = %12.8f, opt = 1\n", v[2]);
    printf("  rmse = %12.4e\n", rmse);

    fprintf(fid,"figure;\n");
    fprintf(fid,"semilogy(u);\n");
    fprintf(fid,"xlabel('iteration');\n");
    fprintf(fid,"ylabel('error');\n");
    fprintf(fid,"title('gradient search results');\n");
    fprintf(fid,"grid on;\n");

    // save sampled data set
    for (i=0; i<num_samples; i++) {
        fprintf(fid,"  x(%4u) = %12.8f;\n", i+1, x[i]);
        fprintf(fid,"  y(%4u) = %12.8f;\n", i+1, y[i]);
        fprintf(fid,"  y_hat(%4u) = %12.8f;\n", i+1, gsfunc(x[i],v));
    }
    fprintf(fid,"figure;\n");
    fprintf(fid,"plot(x,y,'x', x,y_hat,'-');\n");
    fprintf(fid,"xlabel('x');\n");
    fprintf(fid,"ylabel('f(x)');\n");
    fprintf(fid,"grid on;\n");
    fprintf(fid,"legend('data','fit',1);\n");


    fclose(fid);
    printf("results written to %s.\n", OUTPUT_FILENAME);

    gradsearch_destroy(gs);

    return 0;
}

// gradient search fit
float gserror(void * _dataset,
              float * _v,
              unsigned int _n)
{
    struct gsdataset * p = (struct gsdataset *) _dataset;

    float rmse = 0.0f; // mean-squared error
    unsigned int i;
    for (i=0; i<p->n; i++) {
        // compute function estimate
        float y_hat = gsfunc(p->x[i], _v);

        // compute error
        float e = p->y[i] - y_hat;

        // accumulate RMS error
        rmse += e*e;
    }

    // normalize error and return
    rmse = sqrtf(rmse / (float)(p->n));
    return rmse;
}