Edinburgh Speech Tools 2.4-release
 
Loading...
Searching...
No Matches
filter.cc
1/*************************************************************************/
2/* */
3/* Centre for Speech Technology Research */
4/* University of Edinburgh, UK */
5/* Copyright (c) 1994,1995,1996 */
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/* THE UNIVERSITY OF EDINBURGH 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 THE UNIVERSITY OF EDINBURGH 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 : Simon King */
34/* Date : October 1996 */
35/*-----------------------------------------------------------------------*/
36/* Filter functions */
37/* */
38/*=======================================================================*/
39#include <cstdlib>
40#include "EST_math.h"
41#include "sigpr/EST_filter.h"
42#include "sigpr/EST_fft.h"
43#include "EST_wave_aux.h"
44#include "EST_TBuffer.h"
45#include "sigpr/EST_Window.h"
46#include "EST_error.h"
47
48// there are multiple possible styles of this because of different
49// possibilities of where to change the filter coefficients.
50
51void lpc_filter(EST_Wave &sig, EST_FVector &a, EST_Wave &res)
52{
53 int i, j;
54 double s;
55
56 for (i = 0; i < sig.num_samples(); ++i)
57 {
58 s = 0;
59 for (j = 1; j < a.n(); ++j)
60 s += a(j) * (float)sig.a_safe(i - j);
61
62 sig.a(i) = (short) s + res.a(i);
63 }
64}
65
66void inv_lpc_filter(EST_Wave &sig, EST_FVector &a, EST_Wave &res)
67{
68 int i, j;
69 double r;
70
71 for (i = 0; i < a.n(); ++i)
72 {
73 r = sig.a_no_check(i);
74 for (j = 1; j < a.n(); ++j)
75 r -= a.a_no_check(j) * (float)sig.a_safe(i - j);
76 res.a(i) = (short) r;
77 }
78 for (i = a.n(); i < sig.num_samples(); ++i)
79 {
80 r = sig.a_no_check(i);
81 for (j = 1; j < a.n(); ++j)
82 r -= a.a_no_check(j) * (float)sig.a_no_check(i - j);
83 res.a(i) = (short) r;
84 }
85}
86
87/*void filter(EST_FVector &in, EST_FVector &out, EST_FVector &filter)
88{
89 double r;
90 for (int i = 0; i < in.length(); ++i)
91 {
92 r = in.a(i);
93 for (int j = 1; j < filter.length(); ++j)
94 r -= filter(j) * in.a(i - j);
95 out[i] = r;
96 }
97}
98*/
99
100void inv_lpc_filter_ola(EST_Wave &in_sig, EST_Track &lpc, EST_Wave &out_sig)
101{
102
103 int i, j, k, start, end, size;
104 EST_FVector filter;
105 EST_FVector window_vals;
106 EST_Wave in_sub, out_sub;
107
108 // copy attributes and size to waveform, fill with zeros.
109 out_sig.resize(in_sig.num_samples(), 1);
110 out_sig.set_sample_rate(in_sig.sample_rate());
111 out_sig.fill(0);
112
113 for(i = 1; i < lpc.num_frames() - 1; ++i)
114 {
115 start = (int)((float)lpc.t(i - 1) * in_sig.sample_rate());
116 end = (int)((float)lpc.t(i + 1) * in_sig.sample_rate());
117 if (end > out_sig.num_samples())
118 end = out_sig.num_samples();
119 size = end - start;
120
121 lpc.frame(filter, i);
122
123 if (size < filter.n())
124 break; // basically a mismatch between lpc and waveform
125
126 in_sig.sub_wave(in_sub, start, size);
127 out_sub.resize(size);
128
129 inv_lpc_filter(in_sub, filter, out_sub);
130
131
132 int centreIndex = (int)(lpc.t(i) * (float)in_sig.sample_rate());
133
134 EST_Window::make_window(window_vals, size, "hanning", centreIndex-start);
135
136 // printf( "%d %d %d (start centre end)\n", start, centreIndex, end );
137
138 // overlap and add using hanning window on original
139 for (k = 0, j = start; j < end; ++j, ++k){
140 out_sig.a_no_check(j) +=
141 (int)((float)out_sub.a_no_check(k) *
142 window_vals.a_no_check(k));
143 }
144 }
145}
146
147void lpc_filter_1(EST_Track &lpc, EST_Wave & res, EST_Wave &sig)
148{
149 int i, j, k;
150 int start, end;
151 EST_FVector filter;
152 float s;
153// int order = lpc.num_channels() - 1;
154// EST_Wave sig_sub, res_sub;
155
156 sig.resize(res.num_samples());
157 sig.set_sample_rate(res.sample_rate());
158 sig.fill(0);
159
160 for (start = 0, i = 0; i < lpc.num_frames() - 1; ++i)
161 {
162 end = int((lpc.t(i) + lpc.t(i + 1)) * (float)res.sample_rate()) / 2;
163 if (end > res.num_samples())
164 end = res.num_samples();
165
166 lpc.frame(filter, i);
167// res.sub_wave(res_sub, start, (end - start));
168// sig.sub_wave(sig_sub, start, (end - start));
169
170 // this should really be done by the lpc_frame filter function
171 // but it needs access to samples off the start of the frame
172 // in question.
173
174 if (start < filter.n())
175 for (k = start; k < end; ++k)
176 {
177 for (s = 0,j = 1; j < filter.n(); ++j)
178 s += filter.a_no_check(j) * (float)sig.a_safe(k - j);
179 sig.a_no_check(k) = (short) s + res.a_no_check(k);
180 }
181 else
182 for (k = start; k < end; ++k)
183 {
184 s = 0;
185 for (j = 1; j < filter.n(); ++j)
186 s += filter.a_no_check(j) * (float)sig.a_no_check(k - j);
187 sig.a_no_check(k) = (short) s + res.a_no_check(k);
188 }
189 start = end;
190 }
191}
192
193
194
195void lpc_filter_fast(EST_Track &lpc, EST_Wave & res, EST_Wave &sig)
196{
197 // An (unfortunately) faster version of the above. This is about
198 // three time faster than the above. Improvements would need to be
199 // done of signal access to make the above compete.
200 // Note the rescaling of the residual is packaged within here too
201 int i, j, k, m, n;
202 int start, end;
203 float s;
204 int order = lpc.num_channels() - 1;
205 if (order < 0) order = 0; // when lpc as no frames
206 float *buff = walloc(float,res.num_samples()+order);
207 float *filt = walloc(float,order+1);
208 short *residual = res.values().memory();
209
210 sig.resize(res.num_samples(),1,0); // no reseting of values
211 sig.set_sample_rate(res.sample_rate());
212 for (k=0; k<order; k++)
213 buff[k] = 0;
214 for (start = k, m = 0, i = 0; i < lpc.num_frames() - 1; ++i)
215 {
216 end = int((lpc.t(i) + lpc.t(i + 1)) * (float)res.sample_rate()) / 2;
217 if (end > res.num_samples())
218 end = res.num_samples();
219 for (j=1; j<lpc.num_channels(); j++)
220 filt[j]=lpc.a_no_check(i,j);
221 n = j;
222 for (k = start; k < end; ++k,++m)
223 {
224 s = 0;
225 for (j = 1; j < n; ++j)
226 s += filt[j] * buff[k-j];
227 // The 0.5 should be a parameter
228 // buff[k] = s + (residual[m]*0.5);
229 buff[k] = s + residual[m];
230 }
231 start = end;
232 }
233 short *signal = sig.values().memory();
234 for (j=0,i=order; i < k; j++,i++)
235 signal[j] = (short)buff[i];
236 wfree(buff);
237 wfree(filt);
238}
239
240void post_emphasis(EST_Wave &sig, float a)
241{
242 double last=0;
243
244 for (int j = 0; j < sig.num_channels(); ++j)
245 for (int i = 0; i < sig.num_samples(); i++)
246 {
247 sig.a(i, j) = (int)((float)sig.a(i, j) + a * last);
248 last = (float)sig(i, j);
249 // if (absval(sig(i) > maxval)
250 // maxval = absval(fddata[i]);
251 }
252}
253
254void pre_emphasis(EST_Wave &sig, float a)
255{
256 float x = 0.0;
257 float x_1 = 0.0;
258
259 for (int j = 0; j < sig.num_channels(); ++j)
260 for (int i = 0; i < sig.num_samples(); i++)
261 {
262 x = sig.a_no_check(i, j);
263 sig.a_no_check(i, j) =
264 sig.a_no_check(i, j) - int(a * x_1);
265 x_1 = x;
266 }
267}
268
269void pre_emphasis(EST_Wave &sig, EST_Wave &out, float a)
270{
271 out.resize(sig.num_samples(), sig.num_channels());
272
273 for (int j = 0; j < sig.num_channels(); ++j)
274 {
275 out.a_no_check(0, j) = sig.a_no_check(0, j);
276 for (int i = 1; i < sig.num_samples(); i++)
277 out.a_no_check(i, j) =
278 sig.a_no_check(i, j) - int(a * (float)sig.a_no_check(i-1, j));
279 }
280}
281
282void post_emphasis(EST_Wave &sig, EST_Wave &out, float a)
283{
284 out.resize(sig.num_samples(), sig.num_channels());
285
286 for (int j = 0; j < sig.num_channels(); ++j)
287 {
288 out.a_no_check(0, j) = sig.a_no_check(0, j);
289 for (int i = 1; i < sig.num_samples(); i++)
290 out.a_no_check(i, j) =
291 sig.a_no_check(i, j) + int(a * (float)sig.a_no_check(i-1, j));
292 }
293}
294
295void simple_mean_smooth(EST_Wave &c, int n)
296{ // simple mean smoother of order n
297 int i, j, h, k=1;
298 float *a = new float[c.num_samples()];
299 float sum;
300 h = n/2;
301
302 for (i = 0; (i < h); ++i)
303 {
304 k = (i * 2) + 1;
305 sum = 0.0;
306 for (j = 0; (j < k) && (k < c.num_samples()); ++j)
307 sum += c.a_no_check(j);
308 a[i] = sum /(float) k;
309 }
310
311 for (i = h; i < c.num_samples() - h; ++i)
312 {
313 sum = 0.0;
314 for (j = 0; j < n; ++j)
315 sum += c.a_no_check(i - h + j);
316 a[i] = sum /(float) k;
317 }
318
319 for (; i < c.num_samples(); ++i)
320 {
321 k = ((c.num_samples() - i)* 2) -1;
322 sum = 0.0;
323 for (j = 0; j < k; ++j)
324 sum += c.a_no_check(i - (k/2) + j);
325 a[i] = sum /(float) k;
326 }
327
328 for (i = 0; i < c.num_samples(); ++i)
329 c.a_no_check(i) = (short)(a[i] + 0.5);
330
331 delete [] a;
332}
333
334void FIRfilter(EST_Wave &in_sig, const EST_FVector &numerator,
335 int delay_correction)
336{
337 EST_Wave out_sig;
338
339 out_sig.resize(in_sig.num_samples());
340 out_sig.set_sample_rate(in_sig.sample_rate());
341 out_sig.set_file_type(in_sig.file_type());
342
343 FIRfilter(in_sig, out_sig, numerator, delay_correction);
344 in_sig = out_sig;
345}
346
347void FIRfilter(const EST_Wave &in_sig, EST_Wave &out_sig,
348 const EST_FVector &numerator, int delay_correction)
349{
350 if (delay_correction < 0)
351 EST_error("Can't have negative delay !\n");
352
353 if (numerator.n() <= 0)
354 EST_error("Can't filter EST_Wave with given filter");
355
356 int i,j,n=in_sig.num_samples();
357 out_sig.resize(n);
358
359 // could speed up here with three loops :
360 // 1 for first part (filter overlaps start of wave, one 'if')
361 // 2 for middle part (filter always within wave, no 'if's)
362 // 3 for last part (filter overlaps end of wave, one 'if')
363
364 // To make this faster do the float conversion once and hold it
365 // in a conventional array. Note this has been checked to be
366 // safe but if you change the code below you'll have to confirm it
367 // remains safe. If access through FVector was fast I'd use them
368 // but this is about twice as fast.
369 float *in = walloc(float,n);
370 for (i=0; i < n; ++i)
371 in[i] = (float)in_sig.a_no_check(i);
372 float *numer = walloc(float,numerator.n());
373 for (i=0; i < numerator.n(); ++i)
374 numer[i] = numerator.a_no_check(i);
375 float *out = walloc(float,n);
376
377 for (i = 0; i < n; ++i)
378 {
379 out[i] = 0;
380
381 int jlow=0;
382 int jhigh=numerator.n();
383
384 if(i+delay_correction >= n)
385 jlow = i + delay_correction - n + 1;
386
387 if(i+delay_correction - jhigh < 0)
388 jhigh = i + delay_correction;
389
390 for(j=jlow; j<jhigh; j++)
391 if (((i+delay_correction - j) >= 0) &&
392 ((i+delay_correction - j) < n))
393 out[i] += in[i+delay_correction - j] * numer[j];
394 }
395
396 for (i=0; i<n; i++)
397 out_sig.a_no_check(i) = (short)out[i];
398 out_sig.set_sample_rate(in_sig.sample_rate());
399 out_sig.set_file_type(in_sig.file_type());
400
401 wfree(in);
402 wfree(numer);
403 wfree(out);
404}
405
406void FIR_double_filter(EST_Wave &in_sig, EST_Wave &out_sig,
407 const EST_FVector &numerator)
408{
409 out_sig = in_sig;
410 FIRfilter(out_sig, numerator, 0);
411 reverse(out_sig);
412 FIRfilter(out_sig, numerator, 0);
413 reverse(out_sig);
414}
415
416EST_FVector design_FIR_filter(const EST_FVector &frequency_response,
417 int filter_order)
418{
419 // frequency_response contains the desired filter reponse,
420 // on a scale 0...sampling frequency
421
422 // check filter_order is odd
423 if((filter_order & 1) == 0){
424 cerr << "Requested filter order must be odd" << endl;
425 return EST_FVector(0);
426 }
427
428 // check frequency_response has dimension 2^N
429 int N = fastlog2(frequency_response.n());
430 if(frequency_response.n() != (int)pow(float(2.0),(float)N)){
431 cerr << "Desired frequency response must have dimension 2^N" << endl;
432 return EST_FVector(0);
433 }
434
435 int i;
436 EST_FVector filt(frequency_response);
437 EST_FVector dummy(frequency_response.n());
438 for(i=0;i<dummy.n();i++)
439 dummy[i] = 0.0;
440
441 int e=slowIFFT(filt,dummy);
442 if (e != 0)
443 {
444 cerr << "Failed to design filter because FFT failed" << endl;
445 return EST_FVector(0);
446 }
447
448 EST_FVector reduced_filt(filter_order);
449
450 int mid = filter_order/2;
451
452 reduced_filt[mid] = filt(0);
453 for(i=1; i<=mid ;i++)
454 {
455 // Hann window for zero ripple
456 float window = 0.5 + 0.5 * cos(PI*(float)i / (float)mid);
457 reduced_filt[mid+i] = filt(i) * window;
458 reduced_filt[mid-i] = filt(i) * window;
459 }
460
461 return reduced_filt;
462}
463
464
465
466EST_FVector design_high_or_low_pass_FIR_filter(int sample_rate,
467 int cutoff_freq, int order,
468 float gain1, float gain2)
469{
470 // change to bandpass filter .....
471
472 if (sample_rate <= 0){
473 cerr << "Can't design a FIR filter for a sampling rate of "
474 << sample_rate << endl;
475 return EST_FVector(0);
476 }
477
478 int i;
479 int N=10; // good minimum size
480
481 int fft_size = (int)pow(float(2.0), float(N));
482 while(fft_size < order*4){ // rule of thumb !?
483 N++;
484 fft_size = (int)pow(float(2.0), float(N));
485 }
486
487 // freq response is from 0 to sampling freq and therefore
488 // must be symmetrical about 1/2 sampling freq
489
490 EST_FVector freq_resp(fft_size);
491 int normalised_cutoff = (fft_size * cutoff_freq)/sample_rate;
492 for(i=0;i<normalised_cutoff;i++){
493 freq_resp[i] = gain1;
494 freq_resp[fft_size-i-1] = gain1;
495 }
496 for(;i<fft_size/2;i++){
497 freq_resp[i] = gain2;
498 freq_resp[fft_size-i-1] = gain2;
499 }
500
501 return design_FIR_filter(freq_resp, order);
502}
503
504EST_FVector design_lowpass_FIR_filter(int sample_rate, int freq, int order)
505{
506 return design_high_or_low_pass_FIR_filter(sample_rate,
507 freq, order, 1.0, 0.0);
508}
509
510EST_FVector design_highpass_FIR_filter(int sample_rate, int freq, int order)
511{
512 return design_high_or_low_pass_FIR_filter(sample_rate,
513 freq, order, 0.0, 1.0);
514}
515
516void FIRlowpass_filter(const EST_Wave &in_sig, EST_Wave &out_sig,
517 int freq, int order)
518{
519 EST_FVector filt = design_lowpass_FIR_filter(in_sig.sample_rate(),
520 freq, order);
521 FIRfilter(in_sig, out_sig, filt, filt.n()/2);
522}
523
524void FIRlowpass_filter(EST_Wave &in_sig, int freq, int order)
525{
526 EST_FVector filt = design_lowpass_FIR_filter(in_sig.sample_rate(),
527 freq, order);
528 FIRfilter(in_sig, filt, filt.n()/2);
529}
530
531
532void FIRhighpass_filter(EST_Wave &in_sig, int freq, int order)
533{
534 EST_FVector filt = design_highpass_FIR_filter(in_sig.sample_rate(),
535 freq,order);
536 FIRfilter(in_sig, filt, filt.n()/2);
537}
538
539void FIRhighpass_filter(const EST_Wave &in_sig, EST_Wave &out_sig,
540 int freq, int order)
541{
542 EST_FVector filt = design_highpass_FIR_filter(in_sig.sample_rate(),
543 freq,order);
544 FIRfilter(in_sig, out_sig, filt, filt.n()/2);
545}
546
547void FIRlowpass_double_filter(EST_Wave &in_sig, int freq, int order)
548{
549 EST_FVector filt = design_lowpass_FIR_filter(in_sig.sample_rate(),
550 freq, order);
551
552 FIRfilter(in_sig, filt, filt.n()/2);
553 reverse(in_sig);
554 FIRfilter(in_sig, filt, filt.n()/2);
555 reverse(in_sig);
556}
557
558void FIRlowpass_double_filter(const EST_Wave &in_sig, EST_Wave &out_sig,
559 int freq, int order)
560{
561 EST_FVector filt = design_lowpass_FIR_filter(in_sig.sample_rate(),
562 freq, order);
563
564 FIRfilter(in_sig, out_sig, filt, filt.n()/2);
565 reverse(out_sig);
566 FIRfilter(out_sig, filt, filt.n()/2);
567 reverse(out_sig);
568}
569
570void FIRhighpass_double_filter(const EST_Wave &in_sig, EST_Wave &out_sig,
571 int freq, int order)
572{
573 EST_FVector filt = design_highpass_FIR_filter(in_sig.sample_rate(),
574 freq,order);
575
576 FIRfilter(in_sig, out_sig, filt, filt.n()/2);
577 reverse(out_sig);
578 FIRfilter(out_sig, filt, filt.n()/2);
579 reverse(out_sig);
580}
581
582void FIRhighpass_double_filter(EST_Wave &in_sig, int freq, int order)
583{
584 EST_FVector filt = design_highpass_FIR_filter(in_sig.sample_rate(),
585 freq,order);
586
587 FIRfilter(in_sig, filt, filt.n()/2);
588 reverse(in_sig);
589 FIRfilter(in_sig, filt, filt.n()/2);
590 reverse(in_sig);
591}
const T * memory() const
INLINE int n() const
number of items in vector.
INLINE const T & a_no_check(int n) const
read-only const access operator: without bounds checking
float & t(int i=0)
return time position of frame i
Definition EST_Track.h:477
int num_channels() const
return number of channels in track
Definition EST_Track.h:656
int num_frames() const
return number of frames in track
Definition EST_Track.h:650
float & a_no_check(int i, int c=0)
Definition EST_Track.h:419
void frame(EST_FVector &fv, int n, int startf=0, int nf=EST_ALL)
Definition EST_Track.h:209
short & a(int i, int channel=0)
Definition EST_Wave.cc:128
int num_channels() const
return the number of channels in the waveform
Definition EST_Wave.h:145
short & a_safe(int i, int channel=0)
Definition EST_Wave.cc:152
int sample_rate() const
return the sampling rate (frequency)
Definition EST_Wave.h:147
void resize(int num_samples, int num_channels=EST_ALL, int set=1)
resize the waveform
Definition EST_Wave.h:184
int num_samples() const
return the number of samples in the waveform
Definition EST_Wave.h:143
void set_sample_rate(const int n)
Set sampling rate to n
Definition EST_Wave.h:149
static void make_window(EST_TBuffer< float > &window_vals, int size, const char *name, int window_centre)