LPC_and_Formant.cpp 5.4 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186
  1. /* LPC_and_Formant.cpp
  2. *
  3. * Copyright (C) 1994-2018 David Weenink
  4. *
  5. * This code is free software; you can redistribute it and/or modify
  6. * it under the terms of the GNU General Public License as published by
  7. * the Free Software Foundation; either version 2 of the License, or (at
  8. * your option) any later version.
  9. *
  10. * This code is distributed in the hope that it will be useful, but
  11. * WITHOUT ANY WARRANTY; without even the implied warranty of
  12. * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
  13. * General Public License for more details.
  14. *
  15. * You should have received a copy of the GNU General Public License
  16. * along with this work. If not, see <http://www.gnu.org/licenses/>.
  17. */
  18. /*
  19. djmw 20030616 Formant_Frame_into_LPC_Frame: remove formant with f >= Nyquist +
  20. change lpc indexing from -1..m
  21. djmw 20080122 float -> double
  22. */
  23. #include "LPC_and_Formant.h"
  24. #include "LPC_and_Polynomial.h"
  25. #include "NUM2.h"
  26. void Formant_Frame_init (Formant_Frame me, integer nFormants) {
  27. my nFormants = nFormants;
  28. if (nFormants > 0) {
  29. my formant = NUMvector<structFormant_Formant> (1, my nFormants);
  30. }
  31. }
  32. void Formant_Frame_scale (Formant_Frame me, double scale) {
  33. for (integer i = 1; i <= my nFormants; i ++) {
  34. my formant [i].frequency *= scale;
  35. my formant [i].bandwidth *= scale;
  36. }
  37. }
  38. void Roots_into_Formant_Frame (Roots me, Formant_Frame thee, double samplingFrequency, double margin) {
  39. integer n = my max - my min + 1;
  40. autoVEC fc = VECzero (n);
  41. autoVEC bc = VECzero (n);
  42. // Determine the formants and bandwidths
  43. thy nFormants = 0;
  44. double fLow = margin, fHigh = samplingFrequency / 2 - margin;
  45. for (integer i = my min; i <= my max; i ++) {
  46. if (my v [i].im < 0) {
  47. continue;
  48. }
  49. double f = fabs (atan2 (my v [i].im, my v [i].re)) * samplingFrequency / 2.0 / NUMpi;
  50. if (f >= fLow && f <= fHigh) {
  51. /*b = - log (my v [i].re * my v [i].re + my v [i].im * my v [i].im) * samplingFrequency / 2 / NUMpi;*/
  52. double b = - log (dcomplex_abs (my v [i])) * samplingFrequency / NUMpi;
  53. thy nFormants ++;
  54. fc [thy nFormants] = f;
  55. bc [thy nFormants] = b;
  56. }
  57. }
  58. Formant_Frame_init (thee, thy nFormants);
  59. for (integer i = 1; i <= thy nFormants; i ++) {
  60. thy formant [i].frequency = fc [i];
  61. thy formant [i].bandwidth = bc [i];
  62. }
  63. }
  64. void LPC_Frame_into_Formant_Frame (LPC_Frame me, Formant_Frame thee, double samplingPeriod, double margin) {
  65. thy intensity = my gain;
  66. if (my nCoefficients == 0) {
  67. return;
  68. }
  69. autoPolynomial p = LPC_Frame_to_Polynomial (me);
  70. autoRoots r = Polynomial_to_Roots (p.get());
  71. Roots_fixIntoUnitCircle (r.get());
  72. Roots_into_Formant_Frame (r.get(), thee, 1 / samplingPeriod, margin);
  73. }
  74. autoFormant LPC_to_Formant (LPC me, double margin) {
  75. try {
  76. double samplingFrequency = 1.0 / my samplingPeriod;
  77. integer nmax = my maxnCoefficients, err = 0;
  78. integer interval = nmax > 20 ? 1 : 10;
  79. Melder_require (nmax < 100, U"We cannot find the roots of a polynomial of order > 99.");
  80. Melder_require (margin < samplingFrequency / 4, U"Margin should be smaller than ", samplingFrequency / 4, U".");
  81. autoFormant thee = Formant_create (my xmin, my xmax, my nx, my dx, my x1, (nmax + 1) / 2);
  82. autoMelderProgress progress (U"LPC to Formant");
  83. for (integer i = 1; i <= my nx; i ++) {
  84. Formant_Frame formant = & thy d_frames [i];
  85. LPC_Frame lpc = & my d_frames [i];
  86. // Initialisation of Formant_Frame is taken care of in Roots_into_Formant_Frame!
  87. try {
  88. LPC_Frame_into_Formant_Frame (lpc, formant, my samplingPeriod, margin);
  89. } catch (MelderError) {
  90. Melder_clearError();
  91. err ++;
  92. }
  93. if ((interval == 1 || (i % interval) == 1)) {
  94. Melder_progress ( (double) i / my nx, U"LPC to Formant: frame ", i, U" out of ", my nx, U".");
  95. }
  96. }
  97. Formant_sort (thee.get());
  98. if (err > 0) {
  99. Melder_warning (err, U" formant frames out of ", my nx, U" are suspect.");
  100. }
  101. return thee;
  102. } catch (MelderError) {
  103. Melder_throw (me, U": no Formant created.");
  104. }
  105. }
  106. void Formant_Frame_into_LPC_Frame (Formant_Frame me, LPC_Frame thee, double samplingPeriod) {
  107. double nyquist = 2.0 / samplingPeriod;
  108. integer n = 2 * my nFormants;
  109. if (my nFormants < 1) {
  110. return;
  111. }
  112. autoVEC lpc = VECraw (n + 2);
  113. lpc [1] = 0;
  114. lpc [2] = 1;
  115. integer m = 2;
  116. for (integer i = 1; i <= my nFormants; i ++) {
  117. double f = my formant [i].frequency;
  118. if (f > nyquist) continue;
  119. // D(z): 1 + p z^-1 + q z^-2
  120. double r = exp (- NUMpi * my formant [i].bandwidth * samplingPeriod);
  121. double p = - 2 * r * cos (2 * NUMpi * f * samplingPeriod);
  122. double q = r * r;
  123. // By defining the two extra elements (0,1) in the lpc vector we can avoid testing
  124. // lpc[3..n+2] will store the coefficients
  125. for (integer j = m + 2; j > 2; j --) {
  126. lpc [j] += p * lpc [j - 1] + q * lpc [j - 2];
  127. }
  128. m += 2;
  129. }
  130. n = thy nCoefficients < n ? thy nCoefficients : n;
  131. for (integer i = 1; i <= n ; i ++) {
  132. thy a [i] = lpc [i + 2];
  133. }
  134. thy gain = my intensity;
  135. }
  136. autoLPC Formant_to_LPC (Formant me, double samplingPeriod) {
  137. try {
  138. autoLPC thee = LPC_create (my xmin, my xmax, my nx, my dx, my x1, 2 * my maxnFormants, samplingPeriod);
  139. for (integer i = 1; i <= my nx; i ++) {
  140. Formant_Frame f = & my d_frames [i];
  141. LPC_Frame lpc = & thy d_frames [i];
  142. integer m = 2 * f -> nFormants;
  143. LPC_Frame_init (lpc, m);
  144. Formant_Frame_into_LPC_Frame (f, lpc, samplingPeriod);
  145. }
  146. return thee;
  147. } catch (MelderError) {
  148. Melder_throw (me, U": no LPC created.");
  149. }
  150. }
  151. /* End of file LPC_and_Formant.cpp */