1 /*************************************************************************/
2 /* */
3 /* Language Technologies Institute */
4 /* Carnegie Mellon University */
5 /* Copyright (c) 2001 */
6 /* All Rights Reserved. */
7 /* */
8 /* Permission is hereby granted, free of charge, to use and distribute */
9 /* this software and its documentation without restriction, including */
10 /* without limitation the rights to use, copy, modify, merge, publish, */
11 /* distribute, sublicense, and/or sell copies of this work, and to */
12 /* permit persons to whom this work is furnished to do so, subject to */
13 /* the following conditions: */
14 /* 1. The code must retain the above copyright notice, this list of */
15 /* conditions and the following disclaimer. */
16 /* 2. Any modifications must be clearly marked as such. */
17 /* 3. Original authors' names are not deleted. */
18 /* 4. The authors' names are not used to endorse or promote products */
19 /* derived from this software without specific prior written */
20 /* permission. */
21 /* */
22 /* CARNEGIE MELLON UNIVERSITY AND THE CONTRIBUTORS TO THIS WORK */
23 /* DISCLAIM ALL WARRANTIES WITH REGARD TO THIS SOFTWARE, INCLUDING */
24 /* ALL IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS, IN NO EVENT */
25 /* SHALL CARNEGIE MELLON UNIVERSITY NOR THE CONTRIBUTORS BE LIABLE */
26 /* FOR ANY SPECIAL, INDIRECT OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES */
27 /* WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN */
28 /* AN ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, */
29 /* ARISING OUT OF OR IN CONNECTION WITH THE USE OR PERFORMANCE OF */
30 /* THIS SOFTWARE. */
31 /* */
32 /*************************************************************************/
33 /* Author: Alan W Black (awb@cs.cmu.edu) */
34 /* Date: June 2001 */
35 /*************************************************************************/
36 /* */
37 /* residuals etc */
38 /* */
39 /*************************************************************************/
40 #include <stdio.h>
41 #include <math.h>
42 #include "cst_wave.h"
43 #include "cst_track.h"
44 #include "cst_sigpr.h"
45
46
inv_lpc_filterd(short * sig,float * a,int order,double * res,int size)47 void inv_lpc_filterd(short *sig, float *a, int order, double *res, int size)
48 {
49 int i, j;
50 double r;
51 for (i = 0; i < order; i++)
52 {
53 r = sig[i];
54 for (j = 1; j < order; j++)
55 /* if (i-j >= 0) */
56 r -= a[j] * (double)sig[i - j];
57 res[i] = r;
58 }
59 for (i = order; i < size; i++)
60 {
61 r = sig[i];
62 for (j = 1; j < order; j++)
63 r -= a[j] * (double)sig[i - j];
64 res[i] = r;
65 }
66 }
67
find_residual(cst_wave * sig,cst_track * lpc)68 cst_wave *find_residual(cst_wave *sig, cst_track *lpc)
69 {
70 cst_wave *res;
71 int i;
72 double *resd;
73 int size,start,end;
74
75 res = new_wave();
76 cst_wave_resize(res,sig->num_samples,1);
77 res->sample_rate = sig->sample_rate;
78 resd = cst_alloc(double,sig->num_samples);
79
80 start = lpc->num_channels;
81 for (i=0; i<lpc->num_frames; i++)
82 {
83 end = (int)((float)sig->sample_rate * lpc->times[i]);
84 size = end - start;
85
86 inv_lpc_filterd(&sig->samples[start],lpc->frames[i],lpc->num_channels,
87 &resd[start],size);
88
89 start = end;
90 }
91
92 for (i=0; i<res->num_samples; i++)
93 res->samples[i] = resd[i];
94
95 cst_free(resd);
96
97 return res;
98 }
99
lpc_filterd(short * res,float * a,int order,double * sig,int size)100 void lpc_filterd(short *res, float *a, int order, double *sig, int size)
101 {
102 int i, j;
103 double r;
104 for (i = 0; i < order; i++)
105 {
106 r = res[i];
107 for (j = 1; j < order; j++)
108 /* if (i-j >= 0) */
109 r += a[j] * (double)sig[i - j];
110 sig[i] = r;
111 }
112 for (i = order; i < size; i++)
113 {
114 r = res[i];
115 for (j = 1; j < order; j++)
116 r += a[j] * (double)sig[i - j];
117 sig[i] = r;
118 }
119 }
120
reconstruct_wave(cst_wave * res,cst_track * lpc)121 cst_wave *reconstruct_wave(cst_wave *res, cst_track *lpc)
122 {
123 cst_wave *sig;
124 int i;
125 double *sigd;
126 int start, end, size;
127
128 sig = new_wave();
129 cst_wave_resize(sig,res->num_samples,1);
130 sig->sample_rate = res->sample_rate;
131 sigd = cst_alloc(double,sig->num_samples);
132
133 start = 0;
134 for (i=0; i<lpc->num_frames; i++)
135 {
136 end = (int)((float)sig->sample_rate * lpc->times[i]);
137 size = end - start;
138
139 lpc_filterd(&res->samples[start],
140 lpc->frames[i],lpc->num_channels,
141 &sigd[start],size);
142
143 start = end;
144 }
145
146 for (i=0; i<sig->num_samples; i++)
147 sig->samples[i] = sigd[i];
148
149 cst_free(sigd);
150
151 return sig;
152 }
153
154 float lpc_min = -2.709040;
155 float lpc_range = 2.328840+2.709040;
156
lpc_filterde(unsigned char * resulaw,unsigned short * as,int order,short * sig,int size)157 void lpc_filterde(unsigned char *resulaw,
158 unsigned short *as, int order,
159 short *sig,
160 int size)
161 {
162 int i, j;
163 float r;
164 float *a;
165 float *sigf;
166
167 a = cst_alloc(float,order);
168 for (i=1; i<order; i++)
169 a[i] = (((float)as[i]/65535.0)*lpc_range)+lpc_min;
170 sigf = cst_alloc(float,size);
171 for (i=0; i<size; i++)
172 sigf[i] = sig[i];
173
174 for (i = 0; i < order; i++)
175 {
176 r = cst_ulaw_to_short(resulaw[i]);
177 for (j = 1; j < order; j++)
178 /* if (i-j >= 0) */
179 r += a[j] * (float)sigf[i - j];
180 sigf[i] = r;
181 }
182 for (i = order; i < size; i++)
183 {
184 r = cst_ulaw_to_short(resulaw[i]);
185 for (j = 1; j < order; j++)
186 r += a[j] * (float)sigf[i - j];
187 sigf[i] = r;
188 }
189
190 for (i=0; i<size; i++)
191 sig[i] = sigf[i];
192
193 cst_free(sigf);
194 /* cst_free(a); */
195 }
196
reconstruct_wavee(cst_wave * res,cst_track * lpc)197 cst_wave *reconstruct_wavee(cst_wave *res, cst_track *lpc)
198 {
199 cst_wave *sig;
200 int i,j;
201 short *sigd;
202 int start, end, size;
203 unsigned short *as;
204 unsigned char *resulaw;
205
206 sig = new_wave();
207 cst_wave_resize(sig,res->num_samples,1);
208 sig->sample_rate = res->sample_rate;
209 sigd = cst_alloc(short,sig->num_samples);
210 as = cst_alloc(unsigned short,lpc->num_channels);
211
212 resulaw = cst_alloc(unsigned char, res->num_samples);
213 for (i=0; i<res->num_samples; i++)
214 resulaw[i] = cst_short_to_ulaw(res->samples[i]);
215
216 start = 0;
217 for (i=0; i<lpc->num_frames; i++)
218 {
219 end = (int)((float)sig->sample_rate * lpc->times[i]);
220 size = end - start;
221
222 for (j=1; j<lpc->num_channels; j++)
223 as[j] = (((lpc->frames[i][j]-lpc_min)/lpc_range)*65535.0);
224
225 lpc_filterde(&resulaw[start],
226 as,lpc->num_channels,
227 &sigd[start],size);
228
229 start = end;
230 }
231
232 for (i=0; i<sig->num_samples; i++)
233 sig->samples[i] = sigd[i];
234
235 cst_free(sigd);
236
237 return sig;
238 }
239
dumdum(cst_wave * res,cst_track * lpc)240 cst_wave *dumdum(cst_wave *res,cst_track *lpc)
241 {
242 cst_lpcres *lpcres;
243 int i,j,s;
244 unsigned short *bbb;
245 int start,end,size;
246 FILE *ofd;
247
248 lpcres = new_lpcres();
249 lpcres_resize_frames(lpcres,lpc->num_frames);
250 lpcres->num_channels = lpc->num_channels-1;
251 start = 0;
252 for (i=0; i<lpc->num_frames; i++)
253 {
254 bbb = cst_alloc(unsigned short,lpc->num_channels-1);
255 for (j=1; j<lpc->num_channels; j++)
256 bbb[j-1] = (unsigned short)
257 (((lpc->frames[i][j]-lpc_min)/lpc_range)*65535);
258 lpcres->frames[i] = bbb;
259 end = lpc->times[i]*res->sample_rate;
260 lpcres->times[i] = lpc->times[i];
261 size = end - start;
262 lpcres->sizes[i] = size;
263 start = end;
264 }
265 lpcres_resize_samples(lpcres,res->num_samples);
266 lpcres->lpc_min = lpc_min;
267 lpcres->lpc_range = lpc_range;
268 lpcres->sample_rate = res->sample_rate;
269 for (i=0; i<res->num_samples; i++)
270 lpcres->residual[i] = cst_short_to_ulaw(res->samples[i]);
271
272 ofd = fopen("lpc1.lpc","w");
273 for (s=0,i=0; i<lpcres->num_frames; i++)
274 {
275 fprintf(ofd,"%d %d %d\n",i,lpcres->times[i],lpcres->sizes[i]);
276 for (j=0; j < lpcres->num_channels; j++)
277 fprintf(ofd,"%d ",lpcres->frames[i][j]);
278 fprintf(ofd,"\n");
279 for (j=0; j < lpcres->sizes[i]; j++,s++)
280 fprintf(ofd,"%d ",lpcres->residual[s]);
281 fprintf(ofd,"\n");
282 }
283 fclose(ofd);
284
285 return lpc_resynth(lpcres);
286 }
287
main(int argc,char ** argv)288 int main(int argc, char **argv)
289 {
290 /* need that cool argument parser */
291 cst_wave *w,*res,*sig,*w2, *w3;
292 cst_track *lpc;
293 int i;
294 double r;
295
296 if (argc != 3)
297 {
298 fprintf(stderr,"usage: lpc_test2 LPC WAVEFILE\n");
299 return 1;
300 }
301
302 lpc = new_track();
303 cst_track_load_est(lpc,argv[1]);
304 sig = new_wave();
305 cst_wave_load_riff(sig,argv[2]);
306
307 res = find_residual(sig,lpc);
308
309 w = reconstruct_wave(res,lpc);
310
311 cst_wave_save_riff(res,"res.wav");
312 cst_wave_save_riff(w,"new.wav");
313
314 w2 = reconstruct_wavee(res,lpc);
315 cst_wave_save_riff(w2,"newe.wav");
316
317 w3 = dumdum(res,lpc);
318 cst_wave_save_riff(w3,"newr.wav");
319
320 for (r=0.0,i=0; i<sig->num_samples; i++)
321 r += ((float)sig->samples[i]-(float)w->samples[i]) *
322 ((float)sig->samples[i]-(float)w->samples[i]);
323 r /= sig->num_samples;
324 printf("orig/new %f\n",sqrt(r));
325 for (r=0.0,i=0; i<sig->num_samples; i++)
326 r += ((float)sig->samples[i]-(float)w2->samples[i]) *
327 ((float)sig->samples[i]-(float)w2->samples[i]);
328 r /= sig->num_samples;
329 printf("orig/newe %f\n",sqrt(r));
330 for (r=0.0,i=0; i<sig->num_samples; i++)
331 r += ((float)sig->samples[i]-(float)w3->samples[i]) *
332 ((float)sig->samples[i]-(float)w3->samples[i]);
333 r /= sig->num_samples;
334 printf("orig/newr %f\n",sqrt(r));
335 return 0;
336 }
337