History

Move FIR/IIR/Polynomial stuff to its own branch.

Erik de Castro Lopo committed 22 years ago

Commit 72257ee1 parent d981d3111 files · +1 −1344Browse files

Unified · Split

11 files
  1. FIR-POLY/Makefile +0 −40
  2. FIR-POLY/demin.cc +0 −270
  3. FIR-POLY/fir_interp.cc +0 −457
  4. FIR-POLY/mag_spectrum.cc +0 −84
  5. FIR-POLY/mag_spectrum.hh +0 −27
  6. FIR-POLY/timeline.m +0 −50
  7. src/Makefile.am +1 −2
  8. src/common.h +0 −6
  9. src/fip_best.h +0 −43
  10. src/samplerate.c +0 −9
  11. src/src_fir_iir_poly.c +0 −356
FIR-POLY/Makefile+0−40
@@ -1,40 +0,0 @@
1TARGETS = fir_interp
2
3CXX = g++
4CXXFLAGS = -g -O3 -Wall -W -Werror -Wstrict-prototypes -Wmissing-prototypes
5
6GMATH_INC = $(shell pkg-config --cflags GMath)
7GMATH_LIB = $(shell pkg-config --libs GMath)
8
9FFTW_INC = $(shell pkg-config --cflags fftw3)
10FFTW_LIB = $(shell pkg-config --libs fftw3)
11
12
13
14all: $(TARGETS)
15
16clean:
17 rm -f $(TARGETS) *.o
18
19check: $(TARGETS)
20 ./fir_interp
21
22fir_interp: fir_interp.o mag_spectrum.o demin.o
23 $(CXX) $(CXXFLAGS) $+ $(GMATH_LIB) $(FFTW_LIB) -o $@
24
25#---------------------------------------------------------------------
26
27fir_interp.o : fir_interp.cc mag_spectrum.hh
28 $(CXX) $(CXXFLAGS) $(GMATH_INC) -c fir_interp.cc -o $@
29
30mag_spectrum.o : mag_spectrum.cc mag_spectrum.hh
31 $(CXX) $(CXXFLAGS) $(FFTW_INC) -c mag_spectrum.cc -o $@
32
33demin.o : demin.cc
34 $(CXX) $(CXXFLAGS) $(GMATH_INC) -c $+ -o $@
35
36# Do not edit or modify anything in this comment block.
37# The arch-tag line is a file identity tag for the GNU Arch
38# revision control system.
39#
40# arch-tag: 34b6b183-dca1-4085-9089-c8b7d5b09099
FIR-POLY/demin.cc+0−270
@@ -1,270 +0,0 @@
1/*
2** Copyright (C) 1998-2004 Erik de Castro Lopo <erikd@mega-nerd.com>
3**
4** The copyright above and this notice must be preserved in all
5** copies of this source code. The copyright above does not
6** evidence any actual or intended publication of this source code.
7**
8** This is unpublished proprietary trade secret of Erik de Castro
9** Lopo. This source code may not be copied, disclosed,
10** distributed, demonstrated or licensed except as authorized by
11** Erik de Castro Lopo.
12**
13** No part of this program or publication may be reproduced,
14** transmitted, transcribed, stored in a retrieval system,
15** or translated into any language or computer language in any
16** form or by any means, electronic, mechanical, magnetic,
17** optical, chemical, manual, or otherwise, without the prior
18** written permission of Erik de Castro Lopo.
19*/
20
21#include <stdio.h>
22#include <unistd.h>
23
24#include <Minimize.hh>
25#include <PsuedoRand.hh>
26
27class DEmin
28{ private :
29 u_int PopSize, FuncCalcs, Generations ;
30
31 PsuedoRand Rand ;
32
33 GMatrix *MainVector, *Vector2 ;
34 double *MainCost, *Cost2 ;
35 double *ScaleFactor, *ScaleFactor2 ;
36 double *Crossover, *Crossover2 ;
37
38 MinFunc Function ;
39
40 public :
41 DEmin (void) ;
42 ~DEmin (void) ;
43
44 void Setup (MinFunc func, const GMatrix& start, u_int popsize, double range, u_long seed) ;
45
46 GMatrix Minimize (double tol, u_int maxgen) ;
47
48 void PrintStats (void) ;
49
50} ; // class DEmin
51
52//======================================================================================
53
54inline
55void Swap (GMatrix **a, GMatrix **b)
56{ GMatrix *temp ;
57 temp = *a ; *a = *b ; *b = temp ;
58} ; //
59
60inline void
61Swap (double **a, double **b)
62{ double *temp ;
63 temp = *a ; *a = *b ; *b = temp ;
64} ; // Swap
65
66inline double
67Max (double a, double b)
68{ return (a > b) ? a : b ;
69} ; // Max
70
71inline double
72Min (double a, double b)
73{ return (a < b) ? a : b ;
74} ; // Max
75
76//======================================================================================
77
78double
79MinDiffEvol (MinFunc Function, GMatrix& start, double Ftol, double range, u_long seedval)
80{ DEmin *pfmin ;
81 u_int dimensions ;
82
83 dimensions = start.GetElements () ;
84 if (dimensions < 20)
85 dimensions = 20 ;
86
87 pfmin = new DEmin ;
88
89 pfmin -> Setup (Function, start, dimensions * 5, range, seedval) ;
90
91 Ftol = (Ftol < 1e-8) ? 1e-8 : Ftol ;
92
93 start = pfmin -> Minimize (Ftol, dimensions * 5000) ;
94
95 delete pfmin ;
96
97 return Function (start) ;
98} ; // MinDiffEvol
99
100
101//======================================================================================
102
103DEmin::DEmin (void)
104{ PopSize = 0 ;
105 FuncCalcs = 0 ;
106 Generations = 0 ;
107 MainVector = NULL ;
108 MainCost = NULL ;
109 ScaleFactor = NULL ;
110 Crossover = NULL ;
111 Vector2 = NULL ;
112 Cost2 = NULL ;
113 ScaleFactor2 = NULL ;
114 Crossover2 = NULL ;
115} ; // DEmin constructor
116
117DEmin::~DEmin (void)
118{ if (PopSize)
119 { delete [] MainVector ;
120 delete MainCost ;
121 delete ScaleFactor ;
122 delete Crossover ;
123 delete [] Vector2 ;
124 delete Cost2 ;
125 delete ScaleFactor2 ;
126 delete Crossover2 ;
127 } ;
128 PopSize = 0 ;
129} ; // DEmin destructor
130
131void
132DEmin::Setup (MinFunc func, const GMatrix& start, u_int popsize, double range, u_long seed)
133{ u_int k ;
134 u_short rows, cols ;
135
136 Generations = 0 ;
137 Function = func ;
138
139 PopSize = (popsize < 2) ? popsize + 2 : popsize ;
140
141 MainVector = new GMatrix [PopSize] ;
142 MainCost = new double [PopSize] ;
143 ScaleFactor = new double [PopSize] ;
144 Crossover = new double [PopSize] ;
145 Vector2 = new GMatrix [PopSize] ;
146 Cost2 = new double [PopSize] ;
147 ScaleFactor2 = new double [PopSize] ;
148 Crossover2 = new double [PopSize] ;
149
150 Rand.Seed (seed + 1234567) ;
151
152 rows = start.GetRows () ;
153 cols = start.GetCols () ;
154
155 GMatrix::PsuedoRandSeed (123123) ;
156
157 if (IsComplex (start))
158 for (k = 0 ; k < PopSize ; k++)
159 { MainVector [k].SetSizeComplex (rows, cols) ;
160 MainVector [k].PsuedoRandomize () ;
161 MainVector [k] = (2.0 * range) * (MainVector [k] - Complex (0.5, 0.5)) ;
162 }
163 else
164 for (k = 0 ; k < PopSize ; k++)
165 { MainVector [k].SetSizeDouble (rows, cols) ;
166 MainVector [k].PsuedoRandomize () ;
167 MainVector [k] = (2.0 * range) * (MainVector [k] - 0.5) ;
168 } ;
169
170 MainVector [0] = start ; // Must include given vector.
171
172 for (k = 0 ; k < PopSize ; k++)
173 { MainCost [k] = Function (MainVector [k]) ;
174 FuncCalcs ++ ;
175 ScaleFactor [k] = 1.2 * Rand.UDDouble () ;
176 Crossover [k] = Rand.UDDouble () ;
177 } ;
178
179 return ;
180} ; // DEmin::Setup
181
182GMatrix
183DEmin::Minimize (double Ftol, u_int maxgen)
184{ u_int k, a, b, c ;
185 GMatrix temp ;
186 double scalefactor, crossover, cost, improvement, best ;
187
188 improvement = 100 * Ftol ;
189 best = 10e10 ;
190
191 while (Generations < maxgen && best > Ftol)
192 { Swap (&MainVector, &Vector2) ;
193 Swap (&MainCost, &Cost2) ;
194 Swap (&ScaleFactor, &ScaleFactor2) ;
195 Swap (&Crossover, &Crossover2) ;
196
197 improvement = 0.0 ;
198
199 for (k = 0 ; k < PopSize ; k++)
200 { do
201 a = u_int (double (PopSize) * Rand.UDDouble ()) ;
202 while (a == k || a >= PopSize) ;
203 do
204 b = u_int (double (PopSize) * Rand.UDDouble ()) ;
205 while (b == k || b == a || b >= PopSize) ;
206 do
207 c = u_int (double (PopSize) * Rand.UDDouble ()) ;
208 while (c == k || c == a || c == b || c >= PopSize) ;
209
210if (a >= PopSize) { printf ("a - very bad (%u, %u).\n", a, PopSize) ; exit (0) ; } ;
211if (b >= PopSize) { printf ("b - very bad (%u, %u).\n", b, PopSize) ; exit (0) ; } ;
212if (c >= PopSize) { printf ("c - very bad (%u, %u).\n", c, PopSize) ; exit (0) ; } ;
213
214 scalefactor = 0.5 * (ScaleFactor2 [b] + ScaleFactor2 [c]) ;
215
216 temp = Vector2 [a] + scalefactor * (Vector2 [b] - Vector2 [c]) ;
217
218 crossover = 0.5 * (Crossover2 [k] + Crossover2 [a]) ;
219
220 temp.CrossOver (Vector2 [k], crossover) ;
221
222 cost = Function (temp) ;
223 FuncCalcs ++ ;
224
225 if (cost < Cost2 [k])
226 { MainVector [k] = temp ;
227 MainCost [k] = cost ;
228 ScaleFactor [k] = scalefactor ;
229 Crossover [k] = crossover ;
230 improvement = 0.5 * (improvement + Max (improvement, Cost2 [k] - cost)) ;
231 best = Min (best, cost) ;
232 }
233 else
234 { MainVector [k] = Vector2 [k] ;
235 MainCost [k] = Cost2 [k] ;
236 ScaleFactor [k] = ScaleFactor2 [k] ;
237 Crossover [k] = Crossover2 [k] ;
238 } ; // if
239 } ; // for k
240
241 Generations ++ ;
242 } ; // while
243
244 // Now find best vector.
245
246 for (a = 0, k = 1 ; k < PopSize ; k++)
247 if (MainCost [k] < MainCost [a])
248 a = k ;
249 temp = MainVector [a] ;
250
251 return temp ;
252} ; // DEmin::Minimize
253
254void
255DEmin::PrintStats (void)
256{ printf ("Generations : %8d\n", Generations) ;
257 printf ("Function evaluations : %8d\n", FuncCalcs) ;
258 putchar ('\n') ;
259} ; // DEmin::PrintStats
260
261//======================================================================================
262
263
264
265// Do not edit or modify anything in this comment block.
266// The arch-tag line is a file identity tag for the GNU Arch
267// revision control system.
268//
269// arch-tag: e874ad11-7294-4757-9370-fd1e6e1261d0
270
FIR-POLY/fir_interp.cc+0−457
@@ -1,457 +0,0 @@
1/*
2** Copyright (C) 2004 Erik de Castro Lopo <erikd@mega-nerd.com>
3**
4** The copyright above and this notice must be preserved in all
5** copies of this source code. The copyright above does not
6** evidence any actual or intended publication of this source code.
7**
8** This is unpublished proprietary trade secret of Erik de Castro
9** Lopo. This source code may not be copied, disclosed,
10** distributed, demonstrated or licensed except as authorized by
11** Erik de Castro Lopo.
12**
13** No part of this program or publication may be reproduced,
14** transmitted, transcribed, stored in a retrieval system,
15** or translated into any language or computer language in any
16** form or by any means, electronic, mechanical, magnetic,
17** optical, chemical, manual, or otherwise, without the prior
18** written permission of Erik de Castro Lopo.
19*/
20
21
22#include <cstdio>
23#include <cstdlib>
24#include <unistd.h>
25#include <ctime>
26#include <cstring>
27#include <complex>
28
29#include <GMatrix.hh>
30#include <Minimize.hh>
31
32#include "mag_spectrum.hh"
33
34#define UPSAMPLE_RATIO 2
35#define ARRAY_LEN(x) (sizeof (x) / sizeof ((x) [0]))
36
37#define INTERP_LEN 1600
38/*
39** Number of half cycles of impulse response to the left and right
40** of the center.
41** left_half_offset <= right_len
42*/
43
44static int left_len = 150 ;
45static int right_len = 750 ;
46
47
48typedef struct
49{
50 /* Window coefficients (all in range [0, 1]). */
51 double y_left [10], y_right [10] ;
52
53 double sinc_time_fudge ;
54
55} FIR_PARAMS ;
56
57typedef struct
58{
59 union
60 { FIR_PARAMS fir_params ;
61 double data [30] ;
62 } ;
63
64 unsigned middle, total_len ;
65
66 double left_window_poly [20] ;
67 double right_window_poly [20] ;
68
69 double sinc [INTERP_LEN] ;
70 double window [INTERP_LEN] ;
71 double filter [INTERP_LEN] ;
72 double mag [INTERP_LEN / 2] ;
73 double logmag [INTERP_LEN / 2] ;
74 double phase [INTERP_LEN / 2] ;
75} FIR_INTERP ;
76
77static double fir_error (const GMatrix& gm) ;
78static int calc_window (FIR_INTERP *interp, double * error) ;
79static void calc_sinc (FIR_INTERP *interp) ;
80static void calc_filter (FIR_INTERP *interp) ;
81
82static void oct_save (const FIR_INTERP *interp) ;
83static void randomize_data (FIR_INTERP *interp) ;
84
85static double evaluate_filter (FIR_INTERP *interp) ;
86
87int
88main (void)
89{ FIR_INTERP interp ;
90 GMatrix gm_start ;
91 double error ;
92 int param_count ;
93
94 puts ("Running ........") ;
95 param_count = sizeof (FIR_PARAMS) / sizeof (double) ;
96
97 randomize_data (&interp) ;
98
99 gm_start = GMatrix (param_count, 1, interp.data) ;
100
101 fir_error (gm_start) ;
102
103 do
104 error = MinDiffEvol (fir_error, gm_start, 1e-15, 2.0, random ()) ;
105 while (error > 10.0) ;
106
107 printf ("error : %f\n", error) ;
108
109 return 0 ;
110} /* main */
111
112/*==============================================================================
113*/
114
115static inline double
116poly_evaluate (int order, const double *coeff, double x)
117{ double result ;
118 int k ;
119
120 /*
121 ** Use Horner's Rule for evaluating the value of
122 ** the polynomial at x.
123 */
124 result = coeff [order] ;
125 for (k = order - 1 ; k >= 0 ; k--)
126 result = result * x + coeff [k] ;
127
128 return result ;
129} /* poly_evaluate */
130
131static double
132evaluate_filter (FIR_INTERP *interp)
133{ double error = 0.0, temp, stop_max = 0.0, stop_sum = 0.0 ;
134 unsigned k ;
135
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 } ;
145
146 error = stop_max * stop_sum ;
147 return error ;
148} /* evaluate_filter */
149
150static double
151fir_error (const GMatrix& gm)
152{ static FIR_INTERP interp ;
153 static int test_count = 0 ;
154 static double best_error = 1e200 ;
155 double error = 0.0 ;
156 unsigned param_count ;
157
158 memset (&interp, 0, sizeof (interp)) ;
159
160 param_count = sizeof (FIR_PARAMS) / sizeof (double) ;
161
162 if (ARRAY_LEN (interp.data) < param_count)
163 { printf ("\n\nError : ARRAY_LEN (interp.data) < param_count.\n\n") ;
164 exit (1) ;
165 } ;
166
167 if (gm.GetData (param_count, interp.data) != param_count)
168 { printf ("\n\nError : GetData should return %d.\n\n", param_count) ;
169 exit (1) ;
170 } ;
171
172 /* Eval error in sinc_time_fudge. */
173 if (fabs (interp.fir_params.sinc_time_fudge - 1.0) > 0.5)
174 return 1e30 * (fabs (interp.fir_params.sinc_time_fudge - 1.0)) ;
175
176
177 interp.middle = left_len ;
178 interp.total_len = left_len + right_len ;
179
180 if (interp.total_len > ARRAY_LEN (interp.sinc))
181 { printf ("\n\nError : interp.total_len > ARRAY_LEN (interp.sinc).\n") ;
182 exit (1) ;
183 } ;
184
185 if (calc_window (&interp, &error))
186 return error ;
187
188 calc_sinc (&interp) ;
189 calc_filter (&interp) ;
190
191 error = evaluate_filter (&interp) ;
192
193 test_count ++ ;
194// printf ("%s %d : exit\n", __func__, __LINE__) ;
195// exit (1) ;
196
197 if (error < best_error)
198 { oct_save (&interp) ;
199 best_error = error ;
200 printf ("%12d best : %f\n", test_count, best_error) ;
201 } ;
202
203 return error ;
204} /* fir_error */
205
206static int
207calc_window (FIR_INTERP *interp, double * returned_error)
208{ unsigned k ;
209 double error, temp, x ;
210
211 /* Make sure left poly params are in range [0, 1]. */
212 error = 0.0 ;
213 for (k = 0 ; k < ARRAY_LEN (interp->fir_params.y_left) ; k++)
214 { temp = fabs (interp->fir_params.y_left [k] - 0.5) ;
215 if (temp > 0.5 && temp > error)
216 error = temp ;
217 } ;
218
219 if (error > 0.0)
220 { *returned_error = 1e20 * error ;
221 return 1 ;
222 } ;
223
224 /* Make sure right poly params are in range [0, 1]. */
225 error = 0.0 ;
226 for (k = 0 ; k < ARRAY_LEN (interp->fir_params.y_right) ; k++)
227 { temp = fabs (interp->fir_params.y_right [k] - 0.5) ;
228 if (temp > 0.5 && temp > error)
229 error = temp ;
230 } ;
231
232 if (error > 0.0)
233 { *returned_error = 1e18 * error ;
234 return 1 ;
235 } ;
236
237 /* Generate polynomial for left side of window. */
238 interp->left_window_poly [0] = interp->fir_params.y_left [0] ;
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] ;
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] ;
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] ;
245
246 /* Generate left side of window. */
247 error = 0.0 ;
248 for (k = 0 ; k < interp->middle ; k++)
249 { x = k / (interp->middle * 1.0) ;
250 interp->window [k] = poly_evaluate (ARRAY_LEN (interp->left_window_poly) - 1, interp->left_window_poly, x) ;
251 } ;
252
253 /* Ensure start value in range [0, 0,2]. */
254 if (fabs (interp->window [0] - 0.1) > 0.1)
255 { *returned_error = 1e16 * fabs (interp->window [0] - 0.1) ;
256 return 1 ;
257 } ;
258
259 /* Generate polynomial for right side of window. */
260 interp->right_window_poly [0] = 1 ;
261 interp->right_window_poly [1] = 0 ;
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] ;
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] ;
267
268 /* Generate right side of window. */
269 error = 0.0 ;
270 for (k = interp->middle ; k < interp->total_len ; k++)
271 { x = (k - interp->middle) / (1.0 * (interp->total_len - interp->middle)) ;
272 interp->window [k] = poly_evaluate (ARRAY_LEN (interp->right_window_poly) - 1, interp->right_window_poly, x) ;
273 } ;
274
275 /* Ensure end value in range [0, 0,4]. */
276 k = interp->total_len - 1 ;
277 if (fabs (interp->window [k] - 0.1) > 0.1)
278 { *returned_error = 1e16 * fabs (interp->window [k] - 0.1) ;
279 return 1 ;
280 } ;
281
282 /* Ensure whole window is in [0, 1] range. */
283 error = 0.0 ;
284 for (k = 0 ; k < interp->total_len ; k++)
285 { temp = fabs (interp->window [k] - 0.5) ;
286 if (temp > 0.5 && temp > error)
287 error = temp ;
288 } ;
289
290 if (error > 0.0)
291 { *returned_error = 1e16 * error ;
292 return 1 ;
293 } ;
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
321 return 0 ;
322} /* calc_window */
323
324static void
325calc_sinc (FIR_INTERP *interp)
326{ unsigned k ;
327 double x ;
328
329 for (k = 0 ; k < interp->middle ; k++)
330 { x = M_PI * (interp->middle - k) * interp->fir_params.sinc_time_fudge / UPSAMPLE_RATIO ;
331 interp->sinc [k] = sin (x) / x ;
332 } ;
333 interp->sinc [interp->middle] = 1.0 ;
334
335 for (k = interp->middle + 1 ; k < interp->total_len ; k++)
336 { x = M_PI * (k - interp->middle) * interp->fir_params.sinc_time_fudge / UPSAMPLE_RATIO ;
337 interp->sinc [k] = sin (x) / x ;
338 } ;
339
340} /* calc_sinc */
341
342static void
343calc_filter (FIR_INTERP *interp)
344{ unsigned k ;
345 double sum = 0.0 ;
346
347 for (k = 0 ; k < interp->total_len ; k++)
348 { interp->filter [k] = interp->sinc [k] * interp->window [k] ;
349 sum += interp->filter [k] ;
350 } ;
351
352 /* Now normalize. */
353 for (k = 0 ; k < interp->total_len ; k++)
354 interp->filter [k] /= sum ;
355
356} /* calc_filter */
357
358static void
359oct_save (const FIR_INTERP *interp)
360{ const char * filename = "a.mat" ;
361 FILE * file ;
362 unsigned k ;
363
364 unlink (filename) ;
365
366 if ((file = fopen (filename, "w")) == NULL)
367 { printf ("\nError : fopen failed.\n") ;
368 exit (1) ;
369 } ;
370
371 fprintf (file, "# Not created by Octave\n") ;
372
373 fprintf (file, "# name: sinc_time_fudge\n") ;
374 fprintf (file, "# type: scalar\n%16.14f\n", interp->fir_params.sinc_time_fudge) ;
375
376 fprintf (file, "# name: middle\n") ;
377 fprintf (file, "# type: scalar\n%d\n", interp->middle) ;
378
379 fprintf (file, "# name: total_len\n") ;
380 fprintf (file, "# type: scalar\n%d\n", interp->total_len) ;
381
382 fprintf (file, "# name: sinc\n") ;
383 fprintf (file, "# type: matrix\n") ;
384 fprintf (file, "# rows: %d\n", interp->total_len) ;
385 fprintf (file, "# columns: 1\n") ;
386
387 for (k = 0 ; k < interp->total_len ; k++)
388 fprintf (file, "% f\n", interp->sinc [k]) ;
389
390 fprintf (file, "# name: win\n") ;
391 fprintf (file, "# type: matrix\n") ;
392 fprintf (file, "# rows: %d\n", interp->total_len) ;
393 fprintf (file, "# columns: 1\n") ;
394
395 for (k = 0 ; k < interp->total_len ; k++)
396 fprintf (file, "% f\n", interp->window [k]) ;
397
398 fprintf (file, "# name: filt\n") ;
399 fprintf (file, "# type: matrix\n") ;
400 fprintf (file, "# rows: %d\n", interp->total_len) ;
401 fprintf (file, "# columns: 1\n") ;
402
403 for (k = 0 ; k < interp->total_len ; k++)
404 fprintf (file, "% f\n", interp->filter [k] * UPSAMPLE_RATIO) ;
405
406 fprintf (file, "# name: mag\n") ;
407 fprintf (file, "# type: matrix\n") ;
408 fprintf (file, "# rows: %d\n", ARRAY_LEN (interp->mag)) ;
409 fprintf (file, "# columns: 1\n") ;
410
411 for (k = 0 ; k < ARRAY_LEN (interp->mag) ; k++)
412 fprintf (file, "% f\n", interp->mag [k]) ;
413
414 fprintf (file, "# name: logmag\n") ;
415 fprintf (file, "# type: matrix\n") ;
416 fprintf (file, "# rows: %d\n", ARRAY_LEN (interp->logmag)) ;
417 fprintf (file, "# columns: 1\n") ;
418
419 for (k = 0 ; k < ARRAY_LEN (interp->logmag) ; k++)
420 fprintf (file, "% f\n", interp->logmag [k]) ;
421
422 fprintf (file, "# name: phase\n") ;
423 fprintf (file, "# type: matrix\n") ;
424 fprintf (file, "# rows: %d\n", ARRAY_LEN (interp->phase)) ;
425 fprintf (file, "# columns: 1\n") ;
426
427 for (k = 0 ; k < ARRAY_LEN (interp->phase) ; k++)
428 fprintf (file, "% f\n", interp->phase [k]) ;
429
430 fclose (file) ;
431} /* oct_save */
432
433static void
434randomize_data (FIR_INTERP *interp)
435{ FILE * file ;
436 unsigned k, param_count, seed ;
437
438 file = fopen ("/dev/urandom", "r") ;
439 fread (&seed, 1, sizeof (seed), file) ;
440 fclose (file) ;
441
442 srandom (seed) ;
443
444 param_count = sizeof (FIR_PARAMS) / sizeof (double) ;
445
446 for (k = 0 ; k < param_count ; k++)
447 interp->data [k] = 3.0 * ((1.0 * random ()) / INT_MAX - 0.5) ;
448
449} /* randomize_data */
450
451/*
452** Do not edit or modify anything in this comment block.
453** The following line is a file identity tag for the GNU Arch
454** revision control system.
455**
456** arch-tag: 3ce7ca6f-394e-432c-aeaa-f228afc05afd
457*/
FIR-POLY/mag_spectrum.cc+0−84
@@ -1,84 +0,0 @@
1/*
2** Copyright (C) 2004 Erik de Castro Lopo <erikd@mega-nerd.com>
3**
4** The copyright above and this notice must be preserved in all
5** copies of this source code. The copyright above does not
6** evidence any actual or intended publication of this source code.
7**
8** This is unpublished proprietary trade secret of Erik de Castro
9** Lopo. This source code may not be copied, disclosed,
10** distributed, demonstrated or licensed except as authorized by
11** Erik de Castro Lopo.
12**
13** No part of this program or publication may be reproduced,
14** transmitted, transcribed, stored in a retrieval system,
15** or translated into any language or computer language in any
16** form or by any means, electronic, mechanical, magnetic,
17** optical, chemical, manual, or otherwise, without the prior
18** written permission of Erik de Castro Lopo.
19*/
20
21#include <cstdio>
22#include <cstdlib>
23#include <cstring>
24#include <cmath>
25
26#include <fftw3.h>
27
28void
29mag_spectrum (double *input, int len, double *magnitude, double *logmag, double *phase)
30{ static fftw_plan plan = NULL ;
31 static int saved_len = -1 ;
32 static double *fft_data = NULL ;
33
34 double maxval ;
35 int k, pi_factor ;
36
37 if (input == NULL || magnitude == NULL || phase == NULL)
38 { printf ("%s %d : input == NULL || magnitude == NULL || phase == NULL\n", __func__, __LINE__) ;
39 exit (1) ;
40 } ;
41
42 if (saved_len > 0 && saved_len != len)
43 { printf ("%s %d : saved_len != len.\n", __func__, __LINE__) ;
44 exit (1) ;
45 } ;
46
47 if (fft_data == NULL)
48 fft_data = (double *) calloc (sizeof (double), len) ;
49
50 if (plan == NULL)
51 { plan = fftw_plan_r2r_1d (len, input, fft_data, FFTW_R2HC, FFTW_ESTIMATE | FFTW_PRESERVE_INPUT) ;
52 if (plan == NULL)
53 { printf ("%s %d : create plan failed.\n", __func__, __LINE__) ;
54 exit (1) ;
55 } ;
56 } ;
57
58 fftw_execute (plan) ;
59
60 /* (k < N/2 rounded up) */
61 maxval = 0.0 ;
62 for (k = 0 ; k < len / 2 ; k++)
63 { magnitude [k] = sqrt (fft_data [k] * fft_data [k] + fft_data [len - k - 1] * fft_data [len - k - 1]) ;
64 logmag [k] = (magnitude [k] < 1e-10) ? -200.0 : 20 * log10 (magnitude [k]) ;
65 phase [k] = atan2 (fft_data [len - k - 1], fft_data [k]) ;
66 } ;
67
68 pi_factor = 0 ;
69 for (k = 1 ; k < len / 2 ; k++)
70 { phase [k] -= pi_factor * M_PI ;
71 if (fabs (phase [k] - phase [k - 1]) > 3.0)
72 { phase [k] -= 2 * M_PI ;
73 pi_factor += 2 ;
74 }
75 } ;
76
77 return ;
78} /* mag_spectrum */
79
80// Do not edit or modify anything in this comment block.
81// The following line is a file identity tag for the GNU Arch
82// revision control system.
83//
84// arch-tag: 58aa43ee-b8e1-4178-8861-7621399fa049
FIR-POLY/mag_spectrum.hh+0−27
@@ -1,27 +0,0 @@
1/*
2** Copyright (C) 2004 Erik de Castro Lopo <erikd@mega-nerd.com>
3**
4** The copyright above and this notice must be preserved in all
5** copies of this source code. The copyright above does not
6** evidence any actual or intended publication of this source code.
7**
8** This is unpublished proprietary trade secret of Erik de Castro
9** Lopo. This source code may not be copied, disclosed,
10** distributed, demonstrated or licensed except as authorized by
11** Erik de Castro Lopo.
12**
13** No part of this program or publication may be reproduced,
14** transmitted, transcribed, stored in a retrieval system,
15** or translated into any language or computer language in any
16** form or by any means, electronic, mechanical, magnetic,
17** optical, chemical, manual, or otherwise, without the prior
18** written permission of Erik de Castro Lopo.
19*/
20
21void mag_spectrum (double *input, int len, double *magnitude, double *logmag, double *phase) ;
22
23// Do not edit or modify anything in this comment block.
24// The following line is a file identity tag for the GNU Arch
25// revision control system.
26//
27// arch-tag: cc5db1c8-4c67-49f7-9c2f-6c5021a69970
FIR-POLY/timeline.m+0−50
@@ -1,50 +0,0 @@
1function [t, mt, w] = timeline (len, stretch)
2
3
4if (stretch < 0 || stretch > 1)
5 error ("stretch out of range [0, 1].\n");
6 endif
7
8# Mofified time
9
10x_w = 0.2 ;
11x = [-0.5 -0.2 0.2 0.5];
12y = [-0.5 -0.4 -0.1 0.5];
13
14p = polyfit (x, y, 4);
15
16t = linspace (-0.5, 0.5, 100);
17y = polyval (p,t);
18plot (t, y)
19return
20
21mt = ((1 - stretch) * t + stretch * t.^2) - 0.5 ;
22
23
24# Basically a Nuttall window
25ak = [0.42323 0.49755 0.07922] ;
26
27w = zeros (1, len) ;
28
29k = 0 ;
30for a = ak,
31 w = w + a * cos (2 * pi * k * mt) ;
32 k += 1 ;
33 endfor
34
35w = fliplr (w / len) ;
36
37# Linear t
38
39t = linspace (-(len-2)/2, (len-2)/2, len) / len ;
40
41
42end
43
44# Do not edit or modify anything in this comment block.
45# The following line is a file identity tag for the GNU Arch
46# revision control system.
47#
48# arch-tag: 830f701e-5e9c-478a-a929-609513d4b1a6
49
50
src/Makefile.am+1−2
@@ -11,8 +11,7 @@ COEFF_HDRS = high_qual_coeffs.h mid_qual_coeffs.h fastest_coeffs.h · 10 unmodified lines
11 11
12noinst_HEADERS = common.h float_cast.h $(COEFF_HDRS)12noinst_HEADERS = common.h float_cast.h $(COEFF_HDRS)
13 13
14SRC_SOURCES = samplerate.c src_sinc.c $(COEFF_HDRS) src_zoh.c src_linear.c \14SRC_SOURCES = samplerate.c src_sinc.c $(COEFF_HDRS) src_zoh.c src_linear.c
15 src_fir_iir_poly.c
16 15
17# MinGW requires -no-undefined if a DLL is to be built.16# MinGW requires -no-undefined if a DLL is to be built.
18libsamplerate_la_LDFLAGS = -no-undefined version-info @SHARED_VERSION_INFO@ @SHLIB_VERSION_ARG@17libsamplerate_la_LDFLAGS = -no-undefined version-info @SHARED_VERSION_INFO@ @SHLIB_VERSION_ARG@
src/common.h+0−6
@@ -116,12 +116,6 @@ const char* zoh_get_description (int src_enum) ; · 115 unmodified lines
116 116
117int zoh_set_converter (SRC_PRIVATE *psrc, int src_enum) ;117int zoh_set_converter (SRC_PRIVATE *psrc, int src_enum) ;
118 118
119/* In src_fir_iir_poly.c */
120const char* fip_get_name (int src_enum) ;
121const char* fip_get_description (int src_enum) ;
122
123int fip_set_converter (SRC_PRIVATE *psrc, int src_enum) ;
124
125#endif /* COMMON_H_INCLUDED */119#endif /* COMMON_H_INCLUDED */
126 120
127/*121/*
src/fip_best.h+0−43
@@ -1,43 +0,0 @@
1/*
2** Copyright (C) 2004 Erik de Castro Lopo <erikd@mega-nerd.com>
3**
4** This program is free software; you can redistribute it and/or modify
5** it under the terms of the GNU General Public License as published by
6** the Free Software Foundation; either version 2 of the License, or
7** (at your option) any later version.
8**
9** This program is distributed in the hope that it will be useful,
10** but WITHOUT ANY WARRANTY; without even the implied warranty of
11** MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12** GNU General Public License for more details.
13**
14** You should have received a copy of the GNU General Public License
15** along with this program; if not, write to the Free Software
16** Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307, USA.
17*/
18
19
20static const float fip_best_coeffs [250] =
21{ 0.0, 1.2, 2.3
22} ;
23
24/* FIP_COEFF_SET needs to be defined before this header is included. */
25static const FIP_COEFF_SET fip_best =
26{
27 /* count */
28 ARRAY_LEN (fip_best_coeffs),
29
30 /* increment */
31 1,
32
33 /* c */
34 fip_best_coeffs
35} ;
36
37/*
38** Do not edit or modify anything in this comment block.
39** The following line is a file identity tag for the GNU Arch
40** revision control system.
41**
42** arch-tag: 86739f0c-c446-45f9-8ad0-35b6903c6135
43*/
src/samplerate.c+0−9
@@ -313,9 +313,6 @@ src_get_name (int converter_type) · 312 unmodified lines
313 if ((desc = linear_get_name (converter_type)) != NULL)313 if ((desc = linear_get_name (converter_type)) != NULL)
314 return desc ;314 return desc ;
315 315
316 if ((desc = fip_get_name (converter_type)) != NULL)
317 return desc ;
318
319 return NULL ;316 return NULL ;
320} /* src_get_name */317} /* src_get_name */
321 318
@@ -332,9 +329,6 @@ src_get_description (int converter_type) · 10 unmodified lines
332 if ((desc = linear_get_description (converter_type)) != NULL)329 if ((desc = linear_get_description (converter_type)) != NULL)
333 return desc ;330 return desc ;
334 331
335 if ((desc = fip_get_description (converter_type)) != NULL)
336 return desc ;
337
338 return NULL ;332 return NULL ;
339} /* src_get_description */333} /* src_get_description */
340 334
@@ -488,9 +482,6 @@ psrc_set_converter (SRC_PRIVATE *psrc, int converter_type) · 147 unmodified lines
488 if (linear_set_converter (psrc, converter_type) == SRC_ERR_NO_ERROR)482 if (linear_set_converter (psrc, converter_type) == SRC_ERR_NO_ERROR)
489 return SRC_ERR_NO_ERROR ;483 return SRC_ERR_NO_ERROR ;
490 484
491 if (fip_set_converter (psrc, converter_type) == SRC_ERR_NO_ERROR)
492 return SRC_ERR_NO_ERROR ;
493
494 return SRC_ERR_BAD_CONVERTER ;485 return SRC_ERR_BAD_CONVERTER ;
495} /* psrc_set_converter */486} /* psrc_set_converter */
496 487
src/src_fir_iir_poly.c+0−356
@@ -1,356 +0,0 @@
1/*
2** Copyright (C) 2004 Erik de Castro Lopo <erikd@mega-nerd.com>
3**
4** This program is free software; you can redistribute it and/or modify
5** it under the terms of the GNU General Public License as published by
6** the Free Software Foundation; either version 2 of the License, or
7** (at your option) any later version.
8**
9** This program is distributed in the hope that it will be useful,
10** but WITHOUT ANY WARRANTY; without even the implied warranty of
11** MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12** GNU General Public License for more details.
13**
14** You should have received a copy of the GNU General Public License
15** along with this program; if not, write to the Free Software
16** Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307, USA.
17*/
18
19#include <stdio.h>
20#include <stdlib.h>
21#include <string.h>
22
23#include "config.h"
24#include "float_cast.h"
25#include "common.h"
26
27#define FIP_MAGIC_MARKER MAKE_MAGIC ('F', '/', 'I', '/', 'P', ' ')
28
29#define BUFFER_LEN 1024
30
31typedef struct
32{ int count ;
33 int increment ;
34 const float *c ;
35} FIP_COEFF_SET ;
36
37#include "fip_best.h"
38
39typedef struct
40{ int fip_magic_marker ;
41
42 int channels ;
43 long in_count, in_used ;
44 long out_count, out_gen ;
45
46 double src_ratio, input_index ;
47
48 /* The raw coeff set for this converter. */
49 const FIP_COEFF_SET *coeff_set ;
50
51 /* The current FIR coefficents for the times 2 upsampler. */
52 float *coeffs0, * coeffs1 ;
53 int max_half_coeff_len, half_coeff_len ;
54
55 /* Buffer definitions. */
56 float *input_buffer ;
57 int in_buf_len, in_current ;
58
59 float *fir_out ;
60 int fir_out_len, fir_start, fir_end ;
61
62 float *iir_out ;
63 int iir_out_len, iir_start, iir_end ;
64
65 float dummy [1] ;
66} FIR_IIR_POLY ;
67
68static void fip_reset (SRC_PRIVATE *psrc) ;
69static int fip_process (SRC_PRIVATE *psrc, SRC_DATA *data) ;
70
71static int fip_process_var_down (FIR_IIR_POLY * fip, SRC_DATA * data) ;
72static int fip_process_const_down (FIR_IIR_POLY * fip, SRC_DATA * data) ;
73static int fip_process_var_up (FIR_IIR_POLY * fip, SRC_DATA * data) ;
74static int fip_process_const_up (FIR_IIR_POLY * fip, SRC_DATA * data) ;
75
76
77const char*
78fip_get_name (int src_enum)
79{
80 switch (src_enum)
81 { case SRC_FIR_IIR_POLY_BEST :
82 return "Best FIR/IIR/Polynomial Interpolator" ;
83#if 0
84 case SRC_FIR_IIR_POLY_MEDIUM :
85 return "Medium FIR/IIR/Polynomial Interpolator" ;
86
87 case SRC_FIR_IIR_POLY_FASTEST :
88 return "Fastest FIR/IIR/Polynomial Interpolator" ;
89#endif
90 } ;
91
92 return NULL ;
93} /* fip_get_descrition */
94
95const char*
96fip_get_description (int src_enum)
97{
98 switch (src_enum)
99 { case SRC_FIR_IIR_POLY_BEST :
100 return "Three stage FIR/IIR/Polynomial Interpolator, best quality, XXdb SNR, XX% BW." ;
101#if 0
102 case SRC_FIR_IIR_POLY_MEDIUM :
103 return "Three stage FIR/IIR/Polynomial Interpolator, medium quality, XXdb SNR, XX% BW." ;
104
105 case SRC_FIR_IIR_POLY_FASTEST :
106 return "Three stage FIR/IIR/Polynomial Interpolator, fastest, XXdb SNR, XX% BW." ;
107#endif
108 } ;
109
110 return NULL ;
111} /* fip_get_descrition */
112
113int
114fip_set_converter (SRC_PRIVATE *psrc, int src_enum)
115{ FIR_IIR_POLY *fip = NULL, temp_fip ;
116 int buffer_total ;
117
118 if (psrc->private_data != NULL)
119 { fip = (FIR_IIR_POLY *) psrc->private_data ;
120 if (fip->fip_magic_marker != FIP_MAGIC_MARKER)
121 { free (psrc->private_data) ;
122 psrc->private_data = NULL ;
123 } ;
124 } ;
125
126 memset (&temp_fip, 0, sizeof (temp_fip)) ;
127
128 if ((OFFSETOF (FIR_IIR_POLY, dummy) & 0xF) != 0)
129 { printf ("(OFFSETOF (FIR_IIR_POLY, dummy) & 0xF) != 0\n") ;
130 exit (1) ;
131 } ;
132
133 temp_fip.fip_magic_marker = FIP_MAGIC_MARKER ;
134 temp_fip.channels = psrc->channels ;
135
136 psrc->reset = fip_reset ;
137 psrc->process = fip_process ;
138
139 switch (src_enum)
140 { case SRC_FIR_IIR_POLY_BEST :
141 temp_fip.coeff_set = &fip_best ;
142 break ;
143#if 0
144 case SRC_FIR_IIR_POLY_MEDIUM :
145 temp_fip.coeff_set = &fip_medium ;
146 break ;
147
148 case SRC_FIR_IIR_POLY_FASTEST :
149 temp_fip.coeff_set = &fip_fastest ;
150 break ;
151#endif
152 default :
153 return SRC_ERR_BAD_CONVERTER ;
154 } ;
155
156 temp_fip.max_half_coeff_len = SRC_MAX_RATIO * temp_fip.coeff_set->count / temp_fip.coeff_set->increment ;
157 temp_fip.in_buf_len = temp_fip.channels * 2 * temp_fip.max_half_coeff_len ;
158
159 temp_fip.fir_out_len = BUFFER_LEN ;
160 temp_fip.iir_out_len = 2 * BUFFER_LEN ;
161
162 buffer_total = 4 + temp_fip.max_half_coeff_len + psrc->channels * (temp_fip.in_buf_len + temp_fip.fir_out_len + temp_fip.iir_out_len) ;
163
164 if ((fip = calloc (1, sizeof (FIR_IIR_POLY) + sizeof (fip->dummy [0]) * buffer_total)) == NULL)
165 return SRC_ERR_MALLOC_FAILED ;
166
167 *fip = temp_fip ;
168 memset (&temp_fip, 0xEE, sizeof (temp_fip)) ;
169
170 psrc->private_data = fip ;
171
172 /* Allocate buffers. */
173 fip->coeffs0 = fip->dummy ;
174 while ((((unsigned long) fip->coeffs0) & 0xF) > 0)
175 fip->coeffs0 ++ ;
176
177 fip->coeffs1 = fip->coeffs0 + fip->max_half_coeff_len ;
178 fip->input_buffer = fip->coeffs1 + fip->max_half_coeff_len ;
179
180 fip->fir_out = fip->input_buffer + fip->channels * fip->in_buf_len ;
181 fip->iir_out = fip->fir_out + fip->channels * fip->fir_out_len ;
182
183 printf ("total : %d / %d\n", (fip->iir_out + fip->iir_out_len) - fip->dummy, buffer_total) ;
184
185 fip_reset (psrc) ;
186
187 return SRC_ERR_NO_ERROR ;
188} /* fip_set_converter */
189
190static void
191fip_reset (SRC_PRIVATE *psrc)
192{ FIR_IIR_POLY *fip ;
193
194 fip = (FIR_IIR_POLY *) psrc->private_data ;
195 if (fip == NULL)
196 return ;
197
198 memset (fip->input_buffer, 0, fip->channels * fip->in_buf_len * sizeof (fip->dummy [0])) ;
199 memset (fip->fir_out, 0, fip->channels * fip->fir_out_len * sizeof (fip->dummy [0])) ;
200 memset (fip->iir_out, 0, fip->channels * fip->iir_out_len * sizeof (fip->dummy [0])) ;
201
202 fip->in_current = 0 ;
203 fip->fir_start = 0 ;
204 fip->iir_start = fip->iir_end = 0 ;
205
206 fip->src_ratio = fip->input_index = 0.0 ;
207
208 fip->half_coeff_len = 0 ;
209
210 return ;
211} /* fip_reset */
212
213/*========================================================================================
214*/
215
216static int
217fip_process (SRC_PRIVATE *psrc, SRC_DATA *data)
218{ FIR_IIR_POLY *fip ;
219 int error ;
220
221 if (psrc->private_data == NULL)
222 return SRC_ERR_NO_PRIVATE ;
223
224 fip = (FIR_IIR_POLY*) psrc->private_data ;
225
226 /* If there is not a problem, this will be optimised out. */
227 if (sizeof (fip->dummy [0]) != sizeof (data->data_in [0]))
228 return SRC_ERR_SIZE_INCOMPATIBILITY ;
229
230 /* Just after an fip_reset(), the src_ratio field is invalid. */
231 if (fip->src_ratio < 1.0 / SRC_MAX_RATIO)
232 fip->src_ratio = data->src_ratio ;
233
234 fip->in_count = data->input_frames * fip->channels ;
235 fip->out_count = data->output_frames * fip->channels ;
236 fip->in_used = fip->out_gen = 0 ;
237
238 /* Choose a more specialised process() function. */
239 if (fip->src_ratio < 1.0 || data->src_ratio < 1.0)
240 { if (fabs (fip->src_ratio - data->src_ratio) > 1e-20)
241 error = fip_process_var_down (fip, data) ;
242 else
243 error = fip_process_const_down (fip, data) ;
244 }
245 else if (fabs (fip->src_ratio - data->src_ratio) > 1e-20)
246 error = fip_process_var_up (fip, data) ;
247 else
248 error = fip_process_const_up (fip, data) ;
249
250 data->input_frames_used = fip->in_used / fip->channels ;
251 data->output_frames_gen = fip->out_gen / fip->channels ;
252
253 return error ;
254} /* fip_process */
255
256/*----------------------------------------------------------------------------------------
257*/
258
259static void
260fip_generate_fir_current_coeffs (FIR_IIR_POLY * fip)
261{ int k ;
262
263 if (fip->src_ratio >= 1.0)
264 { fip->half_coeff_len = fip->coeff_set->count / fip->coeff_set->increment / 2 ;
265
266 for (k = 0 ; k < fip->half_coeff_len ; k++)
267 { fip->coeffs0 [k] = fip->coeff_set->c [fip->coeff_set->increment * k] ;
268 fip->coeffs1 [k] = fip->coeff_set->c [fip->coeff_set->increment * k + fip->coeff_set->increment / 2] ;
269 } ;
270
271 /* Round it up to a multiple of 4. */
272 while (fip->half_coeff_len & 3)
273 { fip->coeffs0 [fip->half_coeff_len] = 0.0 ;
274 fip->coeffs1 [fip->half_coeff_len] = 0.0 ;
275 fip->half_coeff_len ++ ;
276 } ;
277
278 return ;
279 } ;
280
281 printf ("%s : not implemented yet.\n", __func__) ;
282 exit (1) ;
283} /* fip_generate_fir_current_coeffs */
284
285static int
286fip_process_var_down (FIR_IIR_POLY * fip, SRC_DATA * data)
287{ fip = NULL ;
288 data = NULL ;
289 printf ("%s : not implemented yet.\n", __func__) ;
290 return SRC_ERR_NO_ERROR ;
291} /* fip_process_var_down */
292
293static int
294fip_process_const_down (FIR_IIR_POLY * fip, SRC_DATA * data)
295{
296 if (fip->half_coeff_len == 0)
297 fip_generate_fir_current_coeffs (fip) ;
298
299 data = NULL ;
300 printf ("%s : not implemented yet.\n", __func__) ;
301 return SRC_ERR_NO_ERROR ;
302} /* fip_process_const_down */
303
304static int
305fip_process_var_up (FIR_IIR_POLY * fip, SRC_DATA * data)
306{
307 if (fip->half_coeff_len == 0)
308 fip_generate_fir_current_coeffs (fip) ;
309
310 data = NULL ;
311 printf ("%s : not implemented yet.\n", __func__) ;
312 return SRC_ERR_NO_ERROR ;
313} /* fip_process_var_up */
314
315static int
316fip_process_const_up (FIR_IIR_POLY * fip, SRC_DATA * data)
317{ long max_in_used ;
318 int ch, indx ;
319
320 /* If we don't yet have the coeffs, generate them. */
321 if (fip->half_coeff_len == 0)
322 fip_generate_fir_current_coeffs (fip) ;
323
324 max_in_used = lrintf (fip->in_count / fip->src_ratio - 0.5) ;
325 max_in_used -= max_in_used & 1 ;
326
327 if (fip->channels != 1)
328 { printf ("%s : fip->channels != 1\n", __func__) ;
329 exit (1) ;
330 }
331
332 /* Main processing loop. */
333 while (fip->out_gen < fip->out_count && fip->in_used < max_in_used)
334 {
335 for (ch = 0 ; ch < fip->channels ; ch ++)
336 { indx = (fip->in_current + fip->half_coeff_len + ch) % fip->in_buf_len ;
337 fip->input_buffer [indx] = data->data_in [fip->in_used ++] ;
338 } ;
339
340
341
342
343 fip->in_current = (fip->in_current + 1) % fip->in_buf_len ;
344 } ;
345
346 return SRC_ERR_NO_ERROR ;
347} /* fip_process_const_up */
348
349
350/*
351** Do not edit or modify anything in this comment block.
352** The arch-tag line is a file identity tag for the GNU Arch
353** revision control system.
354**
355** arch-tag: a6c8bad4-740c-4e4f-99f1-e274174ab540
356*/

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