Sound_and_Cepstrum.cpp 3.4 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131
  1. /* Sound_and_Cepstrum.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 20020516 GPL header
  20. djmw 20020529 changed NUMrealft to NUM...RealFastFourierTransform_f
  21. djmw 20041124 Changed call to Sound_to_Spectrum.
  22. */
  23. #include "Sound_and_Cepstrum.h"
  24. #include "Sound_and_Spectrum.h"
  25. #include "Cepstrum_and_Spectrum.h"
  26. #include "NUM2.h"
  27. /*
  28. Algorithm:
  29. J.B. Bednar & T.L. Watt (1985), Calculating the Complex Cepstrum
  30. without Phase Unwrapping or Integration, IEEE Trans. on ASSP 33,
  31. 1014-1017 (Does not work yet).
  32. */
  33. autoCepstrum Sound_to_Cepstrum_bw (Sound me) {
  34. try {
  35. integer nfft = 2;
  36. while (nfft < my nx) nfft *= 2;
  37. double qmax = (my xmax - my xmin) * nfft / my nx;
  38. autoCepstrum thee = Cepstrum_create (qmax, nfft);
  39. autoVEC x = VECraw (nfft);
  40. autoVEC nx = VECraw (nfft);
  41. for (integer i = 1; i <= my nx; i ++) {
  42. x [i] = my z [1] [i];
  43. nx [i] = (i - 1) * x [i];
  44. }
  45. // Step 1: Fourier transform x(n) -> X(f)
  46. // and n*x(n) -> NX(f)
  47. NUMforwardRealFastFourierTransform (x.get());
  48. NUMforwardRealFastFourierTransform (nx.get());
  49. // Step 2: Multiply {X^*(f) * NX(f)} / |X(f)|^2
  50. // Compute Avg (ln |X(f)|) as Avg (ln |X(f)|^2) / 2.
  51. // Treat i=1 separately: x [1] * nx [1] / |x [1]|^2
  52. double lnxa = 0.0;
  53. if (x [1] != 0.0) {
  54. lnxa = 2.0 * log (fabs (x [1]));
  55. x [1] = nx [1] / x [1];
  56. }
  57. if (x [2] != 0.0) {
  58. lnxa = 2.0 * log (fabs (x [2]));
  59. x [2] = nx [2] / x [2];
  60. }
  61. for (integer i = 3; i < nfft; i += 2) {
  62. double xr = x [i], nxr = nx [i];
  63. double xi = x [i + 1], nxi = nx [i + 1];
  64. double xa = xr * xr + xi * xi;
  65. if (xa > 0.0) {
  66. x [i] = (xr * nxr + xi * nxi) / xa;
  67. x [i + 1] = (xr * nxi - xi * nxr) / xa;
  68. lnxa += log (xa);
  69. } else {
  70. x [i] = x [i + 1] = 0.0;
  71. }
  72. }
  73. lnxa /= 2.0 * nfft / 2.0; // TODO
  74. // Step 4: Inverse transform of complex array x
  75. // results in: n * xhat (n)
  76. NUMreverseRealFastFourierTransform (x.get());
  77. // Step 5: Inverse fft-correction factor: 1/nfftd2
  78. // Divide n * xhat (n) by n
  79. for (integer i = 2; i <= my nx; i++) {
  80. thy z [1] [i] = x [i] / ( (i - 1) * nfft);
  81. }
  82. // Step 6: xhat [0] = Avg (ln |X(f)|)
  83. thy z [1] [1] = lnxa;
  84. return thee;
  85. } catch (MelderError) {
  86. Melder_throw (me, U": no Cepstrum created.");
  87. }
  88. }
  89. /* Zijn nog niet elkaars inverse!!!!*/
  90. autoCepstrum Sound_to_Cepstrum (Sound me) {
  91. try {
  92. autoSpectrum spectrum = Sound_to_Spectrum (me, true);
  93. autoCepstrum thee = Spectrum_to_Cepstrum (spectrum.get());
  94. return thee;
  95. } catch (MelderError) {
  96. Melder_throw (me, U": no Cepstrum calculated.");
  97. }
  98. }
  99. autoSound Cepstrum_to_Sound (Cepstrum me) {
  100. try {
  101. autoSpectrum sx = Cepstrum_to_Spectrum (me);
  102. autoSound thee = Spectrum_to_Sound (sx.get());
  103. return thee;
  104. } catch (MelderError) {
  105. Melder_throw (me, U": no Sound calculated.");
  106. }
  107. }
  108. /* End of file Sound_and_Cepstrum.cpp */