123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323 |
- #include "LPC_and_LineSpectralFrequencies.h"
- #include "NUM2.h"
- #include "Polynomial.h"
- static void cos2x (double *g, integer order) {
- for (integer i = 2; i <= order; i ++) {
- for (integer j = order; j > i; j --) {
- g [j - 2] -= g [j];
- }
- g [i - 2] -= 2.0 * g [i];
- }
- }
- static void Polynomial_fromLPC_Frame_lspsum (Polynomial me, LPC_Frame lpc) {
-
-
- integer order = lpc -> nCoefficients, g_order = (order + 1) / 2;
- my coefficients [order + 2] = 1.0;
- for (integer i = 1; i <= order; i ++) {
- my coefficients [order + 2 - i] = lpc -> a [i] + lpc -> a [order + 1 - i];
- }
- my coefficients [1] = 1.0;
- my numberOfCoefficients = order + 2;
- if (order % 2 == 0) {
- Polynomial_divide_firstOrderFactor (me, -1.0, nullptr);
- }
-
-
- for (integer i = 1; i <= g_order + 1; i ++) {
- my coefficients [i] = my coefficients [g_order + i];
- }
- my numberOfCoefficients = g_order + 1;
-
-
- cos2x (& my coefficients [1], g_order);
- }
- static void Polynomial_fromLPC_Frame_lspdif (Polynomial me, LPC_Frame lpc) {
-
- integer order = lpc -> nCoefficients;
- my coefficients [order + 2] = -1.0;
- for (integer i = 1; i <= order; i ++) {
- my coefficients [order + 2 - i] = - lpc -> a [i] + lpc -> a [order + 1 - i];
- }
- my coefficients [1] = 1.0;
- my numberOfCoefficients = order + 2;
- if (order % 2 == 0) {
- Polynomial_divide_firstOrderFactor (me, 1.0, nullptr);
- } else {
- Polynomial_divide_secondOrderFactor (me, 1.0);
- }
-
-
- integer g_order = my numberOfCoefficients / 2;
- for (integer i = 1; i <= g_order + 1; i ++) {
- my coefficients [i] = my coefficients [g_order + i];
- }
- my numberOfCoefficients = g_order + 1;
-
- cos2x (& my coefficients [1], g_order);
- }
- #if 0
- static void Roots_fromPolynomial (Roots me, Polynomial g, integer numberOfDerivatives, double *work) {
- if (numberOfDerivatives < 3) {
- Melder_throw (U"Number of derivatives should be at least 3.");
- }
- double xmin = -1.0, xmax = 1.0;
- integer numberOfRootsFound = 0;
- integer g_order = g -> numberOfCoefficients - 1;
- double *gabs = work, *fact = gabs + g_order + 1, *p2 = fact + numberOfDerivatives + 1;
- double *derivatives = p2 + numberOfDerivatives + 1, *constraints = derivatives + numberOfDerivatives + 1;
- double *intervals = constraints + numberOfDerivatives + 1;
-
-
- fact [0] = p2 [0] = 1.0;
- for (integer j = 1; j <= numberOfDerivatives; j ++) {
- fact [j] = fact [j - 1] * j;
- p2 [j] = p2 [j - 1] * 2.0;
- }
-
-
- for (integer k = 0; k <= g_order; k ++) {
- gabs [k] = fabs (g -> coefficients [k + 1]);
- }
- evaluatePolynomialAndDerivatives (gabs, g_order, 1.0, constraints, numberOfDerivatives);
- intervals [0] = 1.0;
- while (numberOfRootsFound < g_order || xmin == xmax) {
- double dsum1 = 0.0, dsum2 = 0.0;
- double xmid = (xmin + xmax) / 2.0;
- evaluatePolynomialAndDerivatives (g, g_order, xmid, derivatives, numberOfDerivatives);
- double fxmid = derivatives[0], fdxmin = derivatives[1];
- integer j = 1;
- bool rootsOnIntervalPossible_f = true, rootsOnIntervalPossible_df = true;
- while (j <= numberOfDerivatives && (rootsOnIntervalPossible_f || rootsOnIntervalPossible_df)) {
- intervals [j] = intervals [j - 1] * (xmax - xmin);
- integer k = j - 1;
- if (j > 1) {
- dsum1 += fabs (derivatives [k]) * intervals [k] / (p2 [k] * fact [k]);
- }
- if (j > 2) {
- dsum2 += fabs (derivatives [k]) * intervals [k - 1] / (p2 [k - 1] * fact [k - 1]);
- if (rootsOnIntervalPossible_f) {
- double testValue1 = dsum1 + constraints [j] * intervals [j] / (p2 [j] * fact [j]);
- rootsOnIntervalPossible_f = ! (fxmid + testValue1 < 0.0 || fxmid - testValue1 > 0.0);
- }
- if (rootsOnIntervalPossible_df) {
- double testValue2 = dsum2 + constraints [j] * intervals [j - 1] / (p2 [j - 1] * fact [j - 1]);
- rootsOnIntervalPossible_df = ! (fdxmin + testValue2 < 0.0 || fdxmin - testValue2 > 0.0);
- }
- }
- j++;
- }
- if (rootsOnIntervalPossible_f) {
- if (rootsOnIntervalPossible_df) {
- xmax = xmid;
- } else {
- double fxmin = evaluatePolynomial (g, g_order, xmin);
- double fxmax = evaluatePolynomial (g, g_order, xmax);
- if (fxmin * fxmax <= 0.0) {
- double root;
- NUMnrbis (dpoly, xmin, xmax, &poly, &root);
- roots [++numberOfRootsFound] = root;
- } else {
- xmin = xmax; xmax = 1.0;
- }
- }
- } else {
- xmin = xmax; xmax = 1.0;
- }
- }
- }
- #endif
- static integer Roots_fromPolynomial_grid (Roots me, Polynomial thee, double gridSize) {
- Melder_assert (my max >= thy numberOfCoefficients - 1);
- double xmin = thy xmin;
- integer numberOfRootsFound = 0;
- while (xmin < thy xmax && numberOfRootsFound < thy numberOfCoefficients - 1) {
- double xmax = xmin + gridSize;
- xmax = xmax > thy xmax ? thy xmax : xmax;
-
- double root = Polynomial_findOneSimpleRealRoot_ridders (thee, xmin, xmax);
- if (isdefined (root) && (numberOfRootsFound == 0 || my v [numberOfRootsFound].re != root)) {
- my v [++numberOfRootsFound].re = root;
- my v [numberOfRootsFound].im = 0.0;
- }
- xmin = xmax;
- }
-
- return numberOfRootsFound;
- }
- static void LineSpectralFrequencies_Frame_initFromLPC_Frame_grid (LineSpectralFrequencies_Frame me, LPC_Frame thee, Polynomial g1, Polynomial g2, Roots roots, double gridSize, double maximumFrequency) {
-
- LineSpectralFrequencies_Frame_init (me, thy nCoefficients);
- Polynomial_fromLPC_Frame_lspsum (g1, thee);
- integer half_order_g1 = g1 -> numberOfCoefficients - 1;
- Polynomial_fromLPC_Frame_lspdif (g2, thee);
- integer half_order_g2 = g2 -> numberOfCoefficients - 1;
-
- integer numberOfBisections = 0, numberOfRootsFound = 0;
- while (numberOfRootsFound < half_order_g1 && numberOfBisections < 10) {
- numberOfRootsFound = Roots_fromPolynomial_grid (roots, g1, gridSize);
- gridSize *= 0.5; numberOfBisections++;
- }
-
- Melder_require (numberOfBisections < 10, U"Too many bisections.");
-
-
-
- for (integer i = 1; i <= half_order_g1; i ++) {
- my frequencies [2 * i - 1] = acos (roots -> v [half_order_g1 + 1 - i].re / 2.0) / NUMpi * maximumFrequency;
- }
-
-
-
- for (integer i = 1; i <= half_order_g2; i ++) {
- double xmax = roots -> v [half_order_g1 + 1 - i].re;
- double xmin = i == half_order_g1 ? g1 -> xmin : roots -> v [half_order_g1 - i].re;
- double root = Polynomial_findOneSimpleRealRoot_ridders (g2, xmin, xmax);
- if (isdefined (root)) {
- my frequencies [2 * i] = acos (root / 2.0) / NUMpi * maximumFrequency;
- } else {
- my numberOfFrequencies --;
- }
- }
- }
- autoLineSpectralFrequencies LPC_to_LineSpectralFrequencies (LPC me, double gridSize) {
- try {
- if (gridSize == 0.0) {
- gridSize = 0.02;
- }
- double nyquistFrequency = 0.5 / my samplingPeriod;
- autoLineSpectralFrequencies thee = LineSpectralFrequencies_create (my xmin, my xmax, my nx, my dx, my x1, my maxnCoefficients, nyquistFrequency);
- autoPolynomial g1 = Polynomial_create (-2.0, 2.0, my maxnCoefficients + 1);
- autoPolynomial g2 = Polynomial_create (-2.0, 2.0, my maxnCoefficients + 1);
- autoRoots roots = Roots_create ((my maxnCoefficients + 1) / 2);
-
- for (integer iframe = 1; iframe <= my nx; iframe ++) {
- LPC_Frame lpc_frame = & my d_frames [iframe];
- LineSpectralFrequencies_Frame lsf_frame = & thy d_frames [iframe];
-
- LineSpectralFrequencies_Frame_initFromLPC_Frame_grid (lsf_frame, lpc_frame, g1.get(), g2.get(), roots.get(), gridSize, nyquistFrequency);
- }
- return thee;
- } catch (MelderError) {
- Melder_throw (me, U": no LineSpectralFrequencies created.");
- }
- }
- static void LPC_Frame_initFromLineSpectralFrequencies_Frame (LPC_Frame me, LineSpectralFrequencies_Frame thee, Polynomial fs, Polynomial fa, double maximumFrequency) {
- LPC_Frame_init (me, thy numberOfFrequencies);
-
- integer numberOfOmegas = (thy numberOfFrequencies + 1) / 2;
- for (integer i = 1; i <= numberOfOmegas; i ++) {
- double omega = thy frequencies [2 * i -1] / maximumFrequency * NUMpi;
- my a [i] = -2.0 * cos (omega);
- }
- Polynomial_initFromProductOfSecondOrderTerms (fs, {my a.at, numberOfOmegas});
-
- numberOfOmegas = thy numberOfFrequencies / 2;
- for (integer i = 1; i <= numberOfOmegas; i ++) {
- double omega = thy frequencies [2 * i] / maximumFrequency * NUMpi;
- my a [i] = -2.0 * cos (omega);
- }
- Polynomial_initFromProductOfSecondOrderTerms (fa, {my a.at, numberOfOmegas});
-
- if (thy numberOfFrequencies % 2 == 0) {
- Polynomial_multiply_firstOrderFactor (fs, -1.0);
- Polynomial_multiply_firstOrderFactor (fa, 1.0);
- } else {
- Polynomial_multiply_secondOrderFactor (fa, 1.0);
- }
- Melder_assert (fs -> numberOfCoefficients == fa -> numberOfCoefficients);
-
- for (integer i = 1; i <= fs -> numberOfCoefficients - 2; i ++) {
- my a [thy numberOfFrequencies - i + 1] = 0.5 * (fs -> coefficients [i + 1] + fa -> coefficients [i + 1]);
- }
- }
- autoLPC LineSpectralFrequencies_to_LPC (LineSpectralFrequencies me) {
- try {
- autoLPC thee = LPC_create (my xmin, my xmax, my nx, my dx, my x1, my maximumNumberOfFrequencies, 0.5 / my maximumFrequency);
- autoPolynomial fs = Polynomial_create (-1.0, 1.0, my maximumNumberOfFrequencies + 2);
- autoPolynomial fa = Polynomial_create (-1.0, 1.0, my maximumNumberOfFrequencies + 2);
-
- for (integer iframe = 1; iframe <= my nx; iframe ++) {
- LineSpectralFrequencies_Frame lsf_frame = & my d_frames [iframe];
- LPC_Frame lpc_frame = & thy d_frames [iframe];
-
- LPC_Frame_initFromLineSpectralFrequencies_Frame (lpc_frame, lsf_frame, fs.get(), fa.get(), my maximumFrequency);
- }
- return thee;
- } catch (MelderError) {
- Melder_throw (me, U": no LPC created from LineSpectralFrequencies.");
- }
- }
|