History

Larhe update of fir_interp.cc.

Erik de Castro Lopo committed 22 years ago

Commit 1c4bee31 parent dc719a51 file · +90 −40Browse files

Unified · Split

FIR-POLY/fir_interp.cc+90−40
@@ -31,23 +31,24 @@ · 30 unmodified lines
31 31
32#include "mag_spectrum.hh"32#include "mag_spectrum.hh"
33 33
34#define UPSAMPLE_RATIO 2.034#define UPSAMPLE_RATIO 2
35#define ARRAY_LEN(x) (sizeof (x) / sizeof ((x) [0]))35#define ARRAY_LEN(x) (sizeof (x) / sizeof ((x) [0]))
36 36
37#define INTERP_LEN 1600
37/*38/*
38** Number of half cycles of impulse response to the left and right39** Number of half cycles of impulse response to the left and right
39** of the center.40** of the center.
40** left_half_offset <= right_half_cycles41** left_half_offset <= right_len
41*/42*/
42 43
43static int left_half_cycles = 6 ;44static int left_len = 150 ;
44static int right_half_cycles = 20 ;45static int right_len = 750 ;
45 46
46 47
47typedef struct48typedef struct
48{49{
49 /* Window coefficients (all in range [0, 1]). */50 /* Window coefficients (all in range [0, 1]). */
50 double y_left [3], y_right [3] ;51 double y_left [10], y_right [10] ;
51 52
52 double sinc_time_fudge ;53 double sinc_time_fudge ;
53 54
@@ -57,20 +58,20 @@ typedef struct · 3 unmodified lines
57{58{
58 union59 union
59 { FIR_PARAMS fir_params ;60 { FIR_PARAMS fir_params ;
60 double data [20] ;61 double data [30] ;
61 } ;62 } ;
62 63
63 unsigned middle, total_len ;64 unsigned middle, total_len ;
64 65
65 double left_window_poly [5] ;66 double left_window_poly [20] ;
66 double right_window_poly [5] ;67 double right_window_poly [20] ;
67 68
68 double sinc [2048] ;69 double sinc [INTERP_LEN] ;
69 double window [2048] ;70 double window [INTERP_LEN] ;
70 double filter [2048] ;71 double filter [INTERP_LEN] ;
71 double mag [1024] ;72 double mag [INTERP_LEN / 2] ;
72 double logmag [1024] ;73 double logmag [INTERP_LEN / 2] ;
73 double phase [1024] ;74 double phase [INTERP_LEN / 2] ;
74} FIR_INTERP ;75} FIR_INTERP ;
75 76
76static double fir_error (const GMatrix& gm) ;77static double fir_error (const GMatrix& gm) ;
@@ -90,6 +91,7 @@ main (void) · 13 unmodified lines
90 double error ;91 double error ;
91 int param_count ;92 int param_count ;
92 93
94 puts ("Running ........") ;
93 param_count = sizeof (FIR_PARAMS) / sizeof (double) ;95 param_count = sizeof (FIR_PARAMS) / sizeof (double) ;
94 96
95 randomize_data (&interp) ;97 randomize_data (&interp) ;
@@ -128,15 +130,27 @@ poly_evaluate (int order, const double *coeff, double x) · 32 unmodified lines
128 130
129static double131static double
130evaluate_filter (FIR_INTERP *interp)132evaluate_filter (FIR_INTERP *interp)
131{133{ double error = 0.0, temp, stop_max = 0.0, stop_sum = 0.0 ;
134 unsigned k ;
135
132 mag_spectrum (interp->filter, ARRAY_LEN (interp->filter), interp->mag, interp->logmag, interp->phase) ;136 mag_spectrum (interp->filter, ARRAY_LEN (interp->filter), interp->mag, interp->logmag, interp->phase) ;
137
138 for (k = ARRAY_LEN (interp->logmag) / UPSAMPLE_RATIO ; k < ARRAY_LEN (interp->logmag) ; k++)
139 { temp = interp->logmag [k] + 150.0 ;
140 if (temp > stop_max)
141 stop_max = temp ;
142 if (temp > 0.0)
143 stop_sum += temp ;
144 } ;
133 145
134 return 1.0 ;146 error = stop_max * stop_sum ;
147 return error ;
135} /* evaluate_filter */148} /* evaluate_filter */
136 149
137static double150static double
138fir_error (const GMatrix& gm)151fir_error (const GMatrix& gm)
139{ static FIR_INTERP interp ;152{ static FIR_INTERP interp ;
153 static int test_count = 0 ;
140 static double best_error = 1e200 ;154 static double best_error = 1e200 ;
141 double error = 0.0 ;155 double error = 0.0 ;
142 unsigned param_count ;156 unsigned param_count ;
@@ -146,24 +160,22 @@ fir_error (const GMatrix& gm) · 3 unmodified lines
146 param_count = sizeof (FIR_PARAMS) / sizeof (double) ;160 param_count = sizeof (FIR_PARAMS) / sizeof (double) ;
147 161
148 if (ARRAY_LEN (interp.data) < param_count)162 if (ARRAY_LEN (interp.data) < param_count)
149 { printf ("\n\nError : ARRAY_LEN (interp.data) < param_count.\n") ;163 { printf ("\n\nError : ARRAY_LEN (interp.data) < param_count.\n\n") ;
150 exit (1) ;164 exit (1) ;
151 } ;165 } ;
152 166
153 if (gm.GetData (param_count, interp.data) != param_count)167 if (gm.GetData (param_count, interp.data) != param_count)
154 { printf ("\n\nError : GetData should return %d.\n", param_count) ;168 { printf ("\n\nError : GetData should return %d.\n\n", param_count) ;
155 exit (1) ;169 exit (1) ;
156 } ;170 } ;
157 171
158 /* Eval error in sinc_time_fudge. */172 /* Eval error in sinc_time_fudge. */
159 if (interp.fir_params.sinc_time_fudge < 1.0)173 if (fabs (interp.fir_params.sinc_time_fudge - 1.0) > 0.5)
160 return 1e30 * (fabs (interp.fir_params.sinc_time_fudge - 1.0)) ;174 return 1e30 * (fabs (interp.fir_params.sinc_time_fudge - 1.0)) ;
161 175
162 if (interp.fir_params.sinc_time_fudge > 1.25)
163 return 1e30 * (fabs (interp.fir_params.sinc_time_fudge - 1.25)) ;
164 176
165 interp.middle = lrint (floor (left_half_cycles * UPSAMPLE_RATIO / interp.fir_params.sinc_time_fudge)) ;177 interp.middle = left_len ;
166 interp.total_len = interp.middle + lrint (floor (right_half_cycles * UPSAMPLE_RATIO / interp.fir_params.sinc_time_fudge)) ;178 interp.total_len = left_len + right_len ;
167 179
168 if (interp.total_len > ARRAY_LEN (interp.sinc))180 if (interp.total_len > ARRAY_LEN (interp.sinc))
169 { printf ("\n\nError : interp.total_len > ARRAY_LEN (interp.sinc).\n") ;181 { printf ("\n\nError : interp.total_len > ARRAY_LEN (interp.sinc).\n") ;
@@ -178,13 +190,14 @@ fir_error (const GMatrix& gm) · 8 unmodified lines
178 190
179 error = evaluate_filter (&interp) ;191 error = evaluate_filter (&interp) ;
180 192
181 oct_save (&interp) ;193 test_count ++ ;
182 printf ("%s %d : exit\n", __func__, __LINE__) ;194// printf ("%s %d : exit\n", __func__, __LINE__) ;
183 exit (1) ;195// exit (1) ;
184 196
185 if (error < best_error)197 if (error < best_error)
186 { oct_save (&interp) ;198 { oct_save (&interp) ;
187 best_error = error ;199 best_error = error ;
200 printf ("%12d best : %f\n", test_count, best_error) ;
188 } ;201 } ;
189 202
190 return error ;203 return error ;
@@ -195,7 +208,8 @@ calc_window (FIR_INTERP *interp, double * returned_error) · 4 unmodified lines
195{ unsigned k ;208{ unsigned k ;
196 double error, temp, x ;209 double error, temp, x ;
197 210
198 error = 0 ;211 /* Make sure left poly params are in range [0, 1]. */
212 error = 0.0 ;
199 for (k = 0 ; k < ARRAY_LEN (interp->fir_params.y_left) ; k++)213 for (k = 0 ; k < ARRAY_LEN (interp->fir_params.y_left) ; k++)
200 { temp = fabs (interp->fir_params.y_left [k] - 0.5) ;214 { temp = fabs (interp->fir_params.y_left [k] - 0.5) ;
201 if (temp > 0.5 && temp > error)215 if (temp > 0.5 && temp > error)
@@ -207,7 +221,8 @@ calc_window (FIR_INTERP *interp, double * returned_error) · 5 unmodified lines
207 return 1 ;221 return 1 ;
208 } ;222 } ;
209 223
210 error = 0 ;224 /* Make sure right poly params are in range [0, 1]. */
225 error = 0.0 ;
211 for (k = 0 ; k < ARRAY_LEN (interp->fir_params.y_right) ; k++)226 for (k = 0 ; k < ARRAY_LEN (interp->fir_params.y_right) ; k++)
212 { temp = fabs (interp->fir_params.y_right [k] - 0.5) ;227 { temp = fabs (interp->fir_params.y_right [k] - 0.5) ;
213 if (temp > 0.5 && temp > error)228 if (temp > 0.5 && temp > error)
@@ -219,11 +234,14 @@ calc_window (FIR_INTERP *interp, double * returned_error) · 5 unmodified lines
219 return 1 ;234 return 1 ;
220 } ;235 } ;
221 236
222 interp->left_window_poly [0] = interp->fir_params.y_left [0] ;237 /* Generate polynomial for left side of window. */
223 interp->left_window_poly [1] = 13.0/2 - 13.0/2 * interp->fir_params.y_left [0] + 27.0/2 * interp->fir_params.y_left [1] - 27.0/2 * interp->fir_params.y_left [2] ;238 interp->left_window_poly [0] = interp->fir_params.y_left [0] ;
224 interp->left_window_poly [2] = -139.0/4 + 29.0/2 * interp->fir_params.y_left [0] - 189.0/4 * interp->fir_params.y_left [1] + 135.0/2 * interp->fir_params.y_left [2] ;239 interp->left_window_poly [1] = 12.4166666667 - 12.4166666667 * interp->fir_params.y_left [0] + 31.25 * interp->fir_params.y_left [1] - 41.6666666667 * interp->fir_params.y_left [2] + 41.6666666667 * interp->fir_params.y_left [3] - 31.25 * interp->fir_params.y_left [4] ;
225 interp->left_window_poly [3] = 54 - 27.0/2 * interp->fir_params.y_left [0] + 54 * interp->fir_params.y_left [1] - 189.0/2 * interp->fir_params.y_left [2] ;240 interp->left_window_poly [2] = -140.756944444 + 58.2916666667 * interp->fir_params.y_left [0] - 231.770833333 * interp->fir_params.y_left [1] + 413.194444444 * interp->fir_params.y_left [2] - 447.916666667 * interp->fir_params.y_left [3] + 348.958333333 * interp->fir_params.y_left [4] ;
226 interp->left_window_poly [4] = -99.0/4 + 9.0/2 * interp->fir_params.y_left [0] - 81.0/4 * interp->fir_params.y_left [1] + 81.0/2 * interp->fir_params.y_left [2] ;241 interp->left_window_poly [3] = 571.614583333 - 135.416666667 * interp->fir_params.y_left [0] + 662.760416667 * interp->fir_params.y_left [1] - 1395.83333333 * interp->fir_params.y_left [2] + 1682.29166667 * interp->fir_params.y_left [3] - 1385.41666667 * interp->fir_params.y_left [4] ;
242 interp->left_window_poly [4] = -1062.93402778 + 166.666666667 * interp->fir_params.y_left [0] - 917.96875 * interp->fir_params.y_left [1] + 2152.77777778 * interp->fir_params.y_left [2] - 2838.54166667 * interp->fir_params.y_left [3] + 2500. * interp->fir_params.y_left [4] ;
243 interp->left_window_poly [5] = 917.96875 - 104.166666667 * interp->fir_params.y_left [0] + 618.489583333 * interp->fir_params.y_left [1] - 1562.5 * interp->fir_params.y_left [2] + 2213.54166667 * interp->fir_params.y_left [3] - 2083.33333333 * interp->fir_params.y_left [4] ;
244 interp->left_window_poly [6] = -297.309027778 + 26.0416666667 * interp->fir_params.y_left [0] - 162.760416667 * interp->fir_params.y_left [1] + 434.027777778 * interp->fir_params.y_left [2] - 651.041666667 * interp->fir_params.y_left [3] + 651.041666667 * interp->fir_params.y_left [4] ;
227 245
228 /* Generate left side of window. */246 /* Generate left side of window. */
229 error = 0.0 ;247 error = 0.0 ;
@@ -232,16 +250,20 @@ calc_window (FIR_INTERP *interp, double * returned_error) · 2 unmodified lines
232 interp->window [k] = poly_evaluate (ARRAY_LEN (interp->left_window_poly) - 1, interp->left_window_poly, x) ;250 interp->window [k] = poly_evaluate (ARRAY_LEN (interp->left_window_poly) - 1, interp->left_window_poly, x) ;
233 } ;251 } ;
234 252
235 if (fabs (interp->window [0] - 0.2) > 0.2)253 /* Ensure start value in range [0, 0,2]. */
254 if (fabs (interp->window [0] - 0.1) > 0.1)
236 { *returned_error = 1e16 * fabs (interp->window [0] - 0.1) ;255 { *returned_error = 1e16 * fabs (interp->window [0] - 0.1) ;
237 return 1 ;256 return 1 ;
238 } ;257 } ;
239 258
240 interp->right_window_poly [0] = 1 ;259 /* Generate polynomial for right side of window. */
241 interp->right_window_poly [1] = 0 ;260 interp->right_window_poly [0] = 1 ;
242 interp->right_window_poly [2] = -85.0/4 + 27 * interp->fir_params.y_right [0] - 27.0/4 * interp->fir_params.y_right [1] + interp->fir_params.y_right [2] ;261 interp->right_window_poly [1] = 0 ;
243 interp->right_window_poly [3] = 45 - 135.0/2 * interp->fir_params.y_right [0] + 27 * interp->fir_params.y_right [1] - 9.0/2 * interp->fir_params.y_right [2] ;262 interp->right_window_poly [2] = -83.4652777778 + 125. * interp->fir_params.y_right [0] - 62.5 * interp->fir_params.y_right [1] + 27.7777777778 * interp->fir_params.y_right [2] - 7.8125 * interp->fir_params.y_right [3] + interp->fir_params.y_right [4] ;
244 interp->right_window_poly [4] = -99.0/4 + 81.0/2 * interp->fir_params.y_right [0] - 81.0/4 * interp->fir_params.y_right [1] + 9.0/2 * interp->fir_params.y_right [2] ;263 interp->right_window_poly [3] = 446.614583333 - 802.083333333 * interp->fir_params.y_right [0] + 557.291666667 * interp->fir_params.y_right [1] - 270.833333333 * interp->fir_params.y_right [2] + 79.4270833333 * interp->fir_params.y_right [3] - 10.4166666667 * interp->fir_params.y_right [4] ;
264 interp->right_window_poly [4] = -932.725694444 + 1848.95833333 * interp->fir_params.y_right [0] - 1536.45833333 * interp->fir_params.y_right [1] + 850.694444444 * interp->fir_params.y_right [2] - 266.927083333 * interp->fir_params.y_right [3] + 36.4583333333 * interp->fir_params.y_right [4] ;
265 interp->right_window_poly [5] = 865.885416667 - 1822.91666667 * interp->fir_params.y_right [0] + 1692.70833333 * interp->fir_params.y_right [1] - 1041.66666667 * interp->fir_params.y_right [2] + 358.072916667 * interp->fir_params.y_right [3] - 52.0833333333 * interp->fir_params.y_right [4] ;
266 interp->right_window_poly [6] = -297.309027778 + 651.041666667 * interp->fir_params.y_right [0] - 651.041666667 * interp->fir_params.y_right [1] + 434.027777778 * interp->fir_params.y_right [2] - 162.760416667 * interp->fir_params.y_right [3] + 26.0416666667 * interp->fir_params.y_right [4] ;
245 267
246 /* Generate right side of window. */268 /* Generate right side of window. */
247 error = 0.0 ;269 error = 0.0 ;
@@ -250,12 +272,14 @@ calc_window (FIR_INTERP *interp, double * returned_error) · 2 unmodified lines
250 interp->window [k] = poly_evaluate (ARRAY_LEN (interp->right_window_poly) - 1, interp->right_window_poly, x) ;272 interp->window [k] = poly_evaluate (ARRAY_LEN (interp->right_window_poly) - 1, interp->right_window_poly, x) ;
251 } ;273 } ;
252 274
275 /* Ensure end value in range [0, 0,4]. */
253 k = interp->total_len - 1 ;276 k = interp->total_len - 1 ;
254 if (fabs (interp->window [k] - 0.2) > 0.2)277 if (fabs (interp->window [k] - 0.1) > 0.1)
255 { *returned_error = 1e16 * fabs (interp->window [k] - 0.1) ;278 { *returned_error = 1e16 * fabs (interp->window [k] - 0.1) ;
256 return 1 ;279 return 1 ;
257 } ;280 } ;
258 281
282 /* Ensure whole window is in [0, 1] range. */
259 error = 0.0 ;283 error = 0.0 ;
260 for (k = 0 ; k < interp->total_len ; k++)284 for (k = 0 ; k < interp->total_len ; k++)
261 { temp = fabs (interp->window [k] - 0.5) ;285 { temp = fabs (interp->window [k] - 0.5) ;
@@ -268,6 +292,32 @@ calc_window (FIR_INTERP *interp, double * returned_error) · 6 unmodified lines
268 return 1 ;292 return 1 ;
269 } ;293 } ;
270 294
295 /* Monotonicity of left side of window. */
296 error = 0.0 ;
297 for (k = 1 ; k < interp->middle ; k++)
298 { temp = (interp->window [k - 1] - interp->window [k]) ;
299 if (temp > 0.0)
300 error += temp ;
301 } ;
302
303 if (error > 0.0)
304 { *returned_error = 1e14 * error ;
305 return 1 ;
306 } ;
307
308 /* Monotonicity of right side of window. */
309 error = 0.0 ;
310 for (k = interp->middle + 1 ; k < interp->total_len ; k++)
311 { temp = (interp->window [k] - interp->window [k - 1]) ;
312 if (temp > 0.0)
313 error += temp ;
314 } ;
315
316 if (error > 0.0)
317 { *returned_error = 1e14 * error ;
318 return 1 ;
319 } ;
320
271 return 0 ;321 return 0 ;
272} /* calc_window */322} /* calc_window */
273 323

Keyboard shortcuts

General

Show keyboard shortcuts
?
Focus the current search
/

Issues and pull requests

Open filters
f

Issues, pull requests and Inbox

Select next item
j↓
Select previous item
k↑
Open selected item
Enter
Dismiss filters or clear selection
Esc

Repository code

Go to file
t