45 vec
ac2rc(
const vec &ac);
49 vec
is2rc(
const vec &is);
51 vec
rc2ac(
const vec &rc);
53 vec
rc2is(
const vec &rc);
55 vec
autocorr(
const vec &x,
int order)
57 if (order < 0) order = x.size();
63 for (i = 0;i < order + 1;i++) {
65 for (j = 0;j < x.size() - i;j++) {
73 vec
levinson(
const vec &R2,
int order)
76 R[0] = R[0] * (1. + 1.e-9);
78 if (order < 0) order = R.length() - 1;
80 double *
any =
new double[order+1];
81 double *a =
new double[order+1];
92 for (m = 1;m <= order;m++) {
94 for (j = 1;j < m;j++) {
95 s = s + a[j] * R[m-j];
98 k = -(R[m] + s) / alfa;
100 cout <<
"levinson : panic! abs(k)>=1, order " << m <<
". Aborting..." << endl ;
101 for (j = m;j <= order;j++) {
106 for (j = 1;j < m;j++) {
107 any[j] = a[j] + k * a[m-j];
109 for (j = 1;j < m;j++) {
113 alfa = alfa * (1 - k * k);
115 for (j = 0;j < out.length();j++) {
123 vec
lpc(
const vec &x,
int order)
131 int order = a.length() - 1;
132 double alfa, s, *
any =
new double[order+1];
137 it_error_if(a[0] != 1,
"poly2ac : not an lpc filter");
140 for (m = 1;m <= order;m++) {
142 for (j = 1;j < m;j++) {
143 s = s + a[j] * r[m-j];
145 r[m] = -s - alfa * k[m-1];
146 for (j = 1;j < m;j++) {
147 any[j] = a[j] + k[m-1] * a[m-j];
149 for (j = 1;j < m;j++) {
153 alfa = alfa * (1 -
sqr(k[m-1]));
163 int order = a.size() - 1;
165 vec
any(order + 1), aold(a);
167 for (m = order - 1;m > 0;m--) {
169 if (fabs(k[m]) > 1) k[m] = 1.0 / k[m];
170 for (i = 0;i < m;i++) {
171 any[i+1] = (aold[i+1] - aold[m-i] * k[m]) / (1 - k[m] * k[m]);
176 if (fabs(k[0]) > 1) k[0] = 1.0 / k[0];
183 vec a(k.length() + 1),
any(k.length() + 1);
188 for (m = 1;m < k.size();m++) {
190 for (i = 0;i < m;i++) {
191 any[i+1] = a[i+1] + a[m-i] * k[m];
203 for (m = 0;m < k.size();m++) {
204 LAR[m] =
std::log((1 + k[m]) / (1 - k[m]));
209 vec
lar2rc(
const vec &LAR)
214 for (m = 0;m < LAR.size();m++) {
220 double FNevChebP_double(
double x,
const double c[],
int n)
223 double b0 = 0.0, b1 = 0.0, b2 = 0.0;
225 for (i = n - 1; i >= 0; --i) {
228 b0 = 2.0 * x * b1 - b2 + c[i];
230 return (0.5 * (b0 - b2 + c[0]));
233 double FNevChebP(
double x,
const double c[],
int n)
236 double b0 = 0.0, b1 = 0.0, b2 = 0.0;
238 for (i = n - 1; i >= 0; --i) {
241 b0 = 2.0 * x * b1 - b2 + c[i];
243 return (0.5 * (b0 - b2 + c[0]));
248 int np = pc.length() - 1;
251 vec fa((np + 1) / 2 + 1), fb((np + 1) / 2 + 1);
252 vec ta((np + 1) / 2 + 1), tb((np + 1) / 2 + 1);
254 double xlow, xmid, xhigh;
255 double ylow, ymid, yhigh;
262 double DW = (0.02 *
pi);
276 for (i = 1, j = np; i < na; ++i, --j)
277 fa[i] = pc[i] + pc[j];
280 for (i = 1, j = np; i < nb; ++i, --j)
281 fb[i] = pc[i] - pc[j];
284 for (i = 2; i < nb; ++i)
285 fb[i] = fb[i] + fb[i-2];
288 for (i = 1; i < na; ++i) {
289 fa[i] = fa[i] - fa[i-1];
290 fb[i] = fb[i] + fb[i-1];
295 for (i = 1, j = na - 2; i < na; ++i, --j)
299 for (i = 1, j = nb - 2; i < nb; ++i, --j)
307 ylow = FNevChebP_double(xlow, t, n);
312 while (xlow > -1.0 && nf < np) {
315 dx = aa * xhigh * xhigh + ss;
319 ylow = FNevChebP_double(xlow, t, n);
320 if (ylow * yhigh <= 0.0) {
322 for (i = 1; i <= NBIS; ++i) {
325 ymid = FNevChebP_double(xmid, t, n);
326 if (ylow * ymid <= 0.0) {
336 xmid = xlow + dx * ylow / (ylow - yhigh);
345 if (t == ta._data()) {
354 ylow = FNevChebP_double(xlow, t, n);
358 cout <<
"poly2lsf: WARNING: failed to find all lsfs" << endl ;
368 vec p(m + 1), q(m + 1);
371 it_error_if(m % 2 != 0,
"lsf2poly: THIS ROUTINE WORKS ONLY FOR EVEN m");
375 for (i = 0 ; i <= m ; i++) {
380 for (n = 1; n <= mq; n++) {
384 for (i = nor; i >= 2; i--) {
385 q[i] += q[i-2] - c1 * q[i-1];
386 p[i] += p[i-2] - c2 * p[i-1];
391 a[0] = 0.5 * (p[1] + q[1]);
392 for (i = 1, n = 2; i < m ; i++, n++)
393 a[i] = 0.5 * (p[i] + p[n] + q[n] - q[i]);
400 vec c(a.length() - 1);
402 for (
int n = 1;n <= c.length();n++) {
404 for (
int k = 1;k < n;k++) {
405 c[n-1] -= double(k) / n * a[n-k] * c[k-1];
413 it_error_if(num < a.length(),
"a2cepstrum : not allowed cepstrum length");
417 for (n = 1;n < a.length();n++) {
419 for (
int k = 1;k < n;k++) {
420 c[n-1] -= double(k) / n * a[n-k] * c[k-1];
423 for (n = a.length();n <= c.length();n++) {
425 for (
int k = n - a.length() + 1;k < n;k++) {
426 c[n-1] -= double(k) / n * a[n-k] * c[k-1];
434 vec a(c.length() + 1);
437 for (
int n = 1;n <= c.length();n++) {
439 for (
int k = 1;k < n;k++) {
440 a[n] += double(k) / n * a[n-k] * c[k-1];
446 vec
chirp(
const vec &a,
double factor)
448 vec temp(a.length());
452 it_error_if(a[0] != 1,
"chirp : a[0] should be 1");
454 for (i = 1;i < a.length();i++) {
461 vec
schurrc(
const vec &R,
int order)
463 if (order == -1) order = R.length() - 1;
465 vec k(order), scratch(2*order + 2);
473 ep = scratch._data();
474 en = scratch._data() + order + 1;
482 if (en[1] < 1.0) en[1] = 1.0;
486 k[h] = -ep[h+1] / en[1];
487 en[1] = en[1] + k[h] * ep[h+1];
488 if (h == (order - 1)) {
492 ep[order] = ep[order] + k[h] * en[order-h];
494 while (m < (order - 1)) {
496 ex = ep[m] + k[h] * en[m-h];
497 en[m-h] = en[m-h] + k[h] * ep[m];
512 r =
new double[2*M+1];
513 rny =
new double[2*M+1];
515 for (j = 0;j <= M;j++) {
516 r[M-j] = r[M+j] = R[j];
518 for (m = 1;m <= M;m++) {
519 k[m-1] = -r[M+m] / r[M];
520 for (j = -M;j <= M;j++) {
521 rny[M+j] = r[M+j] + k[m-1] * r[M+m-j];
523 for (j = -M;j <= M;j++) {
532 double sd(
const vec &In1,
const vec &In2)
538 double sd(
const vec &In1,
const vec &In2,
double highestfreq)
542 for (
int i = 0;i <
round(highestfreq*129);i++) {
545 S = S * 100 /
round(highestfreq * 129);
#define it_error_if(t, s)
Abort if t is true.
Various functions on vectors and matrices - header file.
vec acos(const vec &x)
Inverse cosine function.
ITPP_EXPORT double sd(const vec &In1, const vec &In2)
Spectral distortion between two vectors, in dB.
ITPP_EXPORT vec rc2lar(const vec &rc)
rc2lar - Reflection coefficients to log area ratios conversion.
vec log10(const vec &x)
log-10 of the elements
ITPP_EXPORT vec cepstrum2poly(const vec &c)
cepstrum2poly - Cepstrum to prediction polynomial conversion.
ITPP_EXPORT vec is2rc(const vec &is)
is2rc - Inverse sine parameters to reflection coefficients conversion.
T sum(const Vec< T > &v)
Sum of all elements in the vector.
ITPP_EXPORT vec ac2rc(const vec &ac)
ac2rc - Autocorrelation sequence to reflection coefficients conversion.
ITPP_EXPORT vec lar2rc(const vec &lar)
lar2rc - Log area ratios to reflection coefficients conversion.
ITPP_EXPORT vec autocorr(const vec &x, int order)
Computes the autocorrelation function.
Definitions of signal processing functions.
ITPP_EXPORT double round(double x)
Round to nearest integer, return result in double.
vec sin(const vec &x)
Sine function.
vec log(const vec &x)
The natural logarithm of the elements.
const double pi
Constant Pi.
vec exp(const vec &x)
Exp of the elements of a vector x.
double energy(const Vec< T > &v)
Calculate the energy: squared 2-norm. energy(v)=sum(abs(v).^2)
ITPP_EXPORT vec rc2ac(const vec &rc)
rc2ac - Reflection coefficients to autocorrelation sequence conversion.
ITPP_EXPORT vec chirp(const vec &a, double factor)
Returns a chirped version of the input vector.
ITPP_EXPORT vec schurrc(const vec &R, int order)
schurrc - Schur algorithm.
ITPP_EXPORT vec rc2is(const vec &rc)
rc2is - Reflection coefficients to inverse sine parameters conversion.
ITPP_EXPORT vec poly2ac(const vec &poly)
poly2ac - Prediction polynomial to autocorrelation sequence conversion.
Miscellaneous statistics functions and classes - header file.
Implementation of linear prediction functions, and conversion between common representations of linea...
vec sqr(const cvec &data)
Absolute square of elements.
ITPP_EXPORT vec ac2poly(const vec &ac)
ac2poly - Autocorrelation sequence to prediction polynomial conversion.
ITPP_EXPORT vec lerouxguegenrc(const vec &R, int order)
Computes reflection coefficients from autocorrelation, using the Le-Roux-Guegen algorithm.
ITPP_EXPORT vec poly2rc(const vec &poly)
poly2rc - Prediction polynomial to reflection coefficients conversion.
ITPP_EXPORT vec poly2cepstrum(const vec &a)
poly2cepstrum - Prediction polynomial to cepstrum conversion.
vec sqrt(const vec &x)
Square root of the elements.
ITPP_EXPORT vec poly2lsf(const vec &poly)
poly2lsf - Prediction polynomial to line spectral frequencies conversion.
bin abs(const bin &inbin)
absolute value of bin
vec cos(const vec &x)
Cosine function.
ITPP_EXPORT vec lpc(const vec &x, int order)
lpc - Linear Predictive Coefficients using autocorrelation method.
ITPP_EXPORT bool any(const bvec &testvec)
Returns true if any element is one and false otherwise.
ITPP_EXPORT vec rc2poly(const vec &rc)
rc2poly - Reflection coefficients to prediction polynomial conversion.
void poly(const vec &r, vec &p)
Create a polynomial of the given rootsCreate a polynomial p with roots r.
ITPP_EXPORT vec levinson(const vec &R2, int order)
Levinson - Levinson-Durbin recursion.
vec filter_spectrum(const vec &a, int nfft)
Power spectrum calculation of a filter.
ITPP_EXPORT vec lsf2poly(const vec &lsf)
lsf2poly - Line spectral frequencies to prediction polynomial conversion.