1 //
2 // lpc_example.c
3 //
4 // This example demonstrates linear prediction in liquid. An input signal
5 // is generated which exhibits a strong temporal correlation. The linear
6 // predictor generates an approximating all-pole filter which minimizes
7 // the squared error between the prediction and the actual output.
8 //
9
10 #include <stdio.h>
11 #include <math.h>
12 #include <complex.h>
13
14 #include "liquid.h"
15
16 #define OUTPUT_FILENAME "lpc_example.m"
17
main()18 int main() {
19 // options
20 unsigned int n = 200; // input sequence length
21 unsigned int p = 4; // prediction filter order
22
23 // create low-pass filter object
24 iirfilt_rrrf f = iirfilt_rrrf_create_lowpass(2, 0.05f);
25 iirfilt_rrrf_print(f);
26
27 unsigned int i;
28
29 // allocate memory for data arrays
30 float y[n]; // input signal (filtered noise)
31 float a_hat[p+1]; // lpc output
32 float g_hat[p+1]; // lpc output
33
34 // generate input signal (filtered noise)
35 for (i=0; i<n; i++)
36 iirfilt_rrrf_execute(f, randnf(), &y[i]);
37
38 // destroy filter object
39 iirfilt_rrrf_destroy(f);
40
41 // run linear prediction algorithm
42 liquid_lpc(y,n,p,a_hat,g_hat);
43
44 // run prediction filter
45 float a_lpc[p+1];
46 float b_lpc[p+1];
47 for (i=0; i<p+1; i++) {
48 a_lpc[i] = (i==0) ? 1.0f : 0.0f;
49 b_lpc[i] = (i==0) ? 0.0f : -a_hat[i];
50 }
51 f = iirfilt_rrrf_create(b_lpc,p+1, a_lpc,p+1);
52 iirfilt_rrrf_print(f);
53 float y_hat[n];
54 for (i=0; i<n; i++)
55 iirfilt_rrrf_execute(f, y[i], &y_hat[i]);
56 iirfilt_rrrf_destroy(f);
57
58 // compute prediction error
59 float err[n];
60 for (i=0; i<n; i++)
61 err[i] = y[i] - y_hat[i];
62
63 // compute autocorrelation of prediction error
64 float lag[n];
65 float rxx[n];
66 for (i=0; i<n; i++) {
67 lag[i] = (float)i;
68 rxx[i] = 0.0f;
69 unsigned int j;
70 for (j=i; j<n; j++)
71 rxx[i] += err[j] * err[j-i];
72 }
73 float rxx0 = rxx[0];
74 for (i=0; i<n; i++)
75 rxx[i] /= rxx0;
76
77 // print results
78 for (i=0; i<p+1; i++)
79 printf(" a[%3u] = %12.8f, g[%3u] = %12.8f\n", i, a_hat[i], i, g_hat[i]);
80
81 printf(" prediction rmse = %12.8f\n", sqrtf(rxx0 / n));
82
83 //
84 // plot results to output file
85 //
86 FILE * fid = fopen(OUTPUT_FILENAME,"w");
87 fprintf(fid,"%% %s : auto-generated file\n", OUTPUT_FILENAME);
88 fprintf(fid,"clear all;\n");
89 fprintf(fid,"close all;\n");
90 fprintf(fid,"\n");
91 fprintf(fid,"p=%u;\n", p);
92 fprintf(fid,"n=%u;\n",n);
93
94 #if 0
95 fprintf(fid,"b = zeros(1,p+1);\n");
96 fprintf(fid,"a = zeros(1,p+1);\n");
97 for (i=0; i<p+1; i++) {
98 fprintf(fid,"b(%4u) = %12.4e;\n", i+1, b_lpc[i]);
99 fprintf(fid,"a(%4u) = %12.4e;\n", i+1, a_lpc[i]);
100 }
101 #endif
102
103 fprintf(fid,"y = zeros(1,n);\n");
104 fprintf(fid,"y_hat = zeros(1,n);\n");
105 fprintf(fid,"lag = zeros(1,n);\n");
106 fprintf(fid,"rxx = zeros(1,n);\n");
107 for (i=0; i<n; i++) {
108 fprintf(fid,"y(%4u) = %12.4e + j*%12.4e;\n", i+1, crealf(y[i]), cimagf(y[i]));
109 fprintf(fid,"y_hat(%4u) = %12.4e + j*%12.4e;\n", i+1, crealf(y_hat[i]), cimagf(y_hat[i]));
110 fprintf(fid,"lag(%4u) = %12.4e + j*%12.4e;\n", i+1, crealf(lag[i]), cimagf(lag[i]));
111 fprintf(fid,"rxx(%4u) = %12.4e + j*%12.4e;\n", i+1, crealf(rxx[i]), cimagf(rxx[i]));
112 }
113
114 // plot output
115 fprintf(fid,"t=0:(n-1);\n");
116 fprintf(fid,"figure;\n");
117 fprintf(fid,"subplot(5,1,1:3);\n");
118 fprintf(fid," hold on;\n");
119 fprintf(fid," plot(t,y, 'LineWidth',1.5,'Color',[0 0.2 0.5]);\n");
120 fprintf(fid," plot(t,y_hat,'LineWidth',1.5,'Color',[0 0.5 0.2]);\n");
121 fprintf(fid," hold off;\n");
122 fprintf(fid," ylabel('signal');\n");
123 fprintf(fid," legend('input','LPC estimate','location','northeast');\n");
124 fprintf(fid," grid on;\n");
125 fprintf(fid,"subplot(5,1,4);\n");
126 fprintf(fid," plot(t,y-y_hat,'LineWidth',2,'Color',[1 1 1]*0.6);\n");
127 fprintf(fid," ylabel('error');\n");
128 fprintf(fid," axis([0 n -0.1 0.1]);\n");
129 fprintf(fid," grid on;\n");
130 fprintf(fid,"subplot(5,1,5);\n");
131 fprintf(fid," plot(t,rxx,'LineWidth',2,'Color',[0.5 0 0]);\n");
132 fprintf(fid," xlabel('time');\n");
133 fprintf(fid," ylabel('r_{xx}(e)');\n");
134 fprintf(fid," axis([0 n -1.0 1]);\n");
135 fprintf(fid," grid on;\n");
136
137 fclose(fid);
138 printf("results written to %s.\n", OUTPUT_FILENAME);
139
140 printf("done.\n");
141 return 0;
142 }
143
144