content: /**
   * Copyright 2007 DFKI GmbH.
   * All Rights Reserved.  Use is subject to license terms.
   * <p/>
   * This file is part of MARY TTS.
   * <p/>
   * MARY TTS is free software: you can redistribute it and/or modify
   * it under the terms of the GNU Lesser General Public License as published by
   * the Free Software Foundation, version 3 of the License.
   * <p/>
   * This program is distributed in the hope that it will be useful,
   * but WITHOUT ANY WARRANTY; without even the implied warranty of
   * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
   * GNU Lesser General Public License for more details.
   * <p/>
   * You should have received a copy of the GNU Lesser General Public License
   * along with this program.  If not, see <http://www.gnu.org/licenses/>.
   */
  package marytts.signalproc.analysis.distance;
  
  import java.io.File;
  import java.io.IOException;
  
  import lib.sound.sampled.AudioInputStream;
  import lib.sound.sampled.AudioSystem;
  import lib.sound.sampled.UnsupportedAudioFileException;
  import marytts.signalproc.analysis.LpcAnalyser;
  import marytts.signalproc.analysis.LpcAnalyser.LpCoeffs;
  import marytts.signalproc.analysis.LsfAnalyser;
  import marytts.signalproc.window.HammingWindow;
  import marytts.util.data.audio.AudioDoubleDataSource;
  import marytts.util.signal.SignalProcUtils;
  
  /**
   * A class for general purpose distance metrics to compare speech/audio signals objectively
   *
   * @author Oytun T&uumlrk
   */
  public class SpectralDistanceMeasures {
      public static double MAX_SPECTRAL_DISTANCE = 1e+10;
  
      /*
       * Inverse harmonic weighting based LSF distance
       * Performs LP analysis by first preemphasizing and windowing speech frames <speechFrame1> and <speechFrame2>
       * Then, computes LSFs for the two speech frames
       * Finally, the perceptual distance between two lsf vectors are estimated.
       * Note that the input speech frames are not changed during distance computation
       * LP order is automatically decided depending on the sampling rate
       * 
       * speechFrame1: First speech frame (not windowed)
       * speechFrame2: Second speech frame (not windowed)
       * samplingRate: Sampling rate in Hz
       * 
       */
      public static double lsfDist(double[] speechFrame1, double[] speechFrame2, int samplingRate) {
          return lsfDist(speechFrame1, speechFrame2, samplingRate, SignalProcUtils.getLPOrder(samplingRate));
      }
  
      /*
       * Inverse harmonic weighting based LSF distance
       * Performs LP analysis by first preemphasizing and windowing speech frames <speechFrame1> and <speechFrame2>
       * Then, computes LSFs for the two speech frames
       * Finally, the perceptual distance between two lsf vectors are estimated.
       * Note that the input speech frames are not changed during distance computation
       * 
       * speechFrame1: First speech frame (not windowed)
       * speechFrame2: Second speech frame (not windowed)
       * samplingRate: Sampling rate in Hz
       * lpOrder: Desired LP analysis order
       * 
       */
      public static double lsfDist(double[] speechFrame1, double[] speechFrame2, int samplingRate, int lpOrder) {
          //Preemphasis
          double[] windowedSpeechFrame1 = new double[speechFrame1.length];
          System.arraycopy(speechFrame1, 0, windowedSpeechFrame1, 0, speechFrame1.length);
          windowedSpeechFrame1 = SignalProcUtils.applyPreemphasis(windowedSpeechFrame1, 0.97);
  
          double[] windowedSpeechFrame2 = new double[speechFrame2.length];
          System.arraycopy(speechFrame2, 0, windowedSpeechFrame2, 0, speechFrame2.length);
          windowedSpeechFrame2 = SignalProcUtils.applyPreemphasis(windowedSpeechFrame2, 0.97);
          //
  
          //Windowing
          HammingWindow w1 = new HammingWindow(speechFrame1.length);
          w1.apply(windowedSpeechFrame1, 0);
          HammingWindow w2 = new HammingWindow(speechFrame2.length);
          w2.apply(windowedSpeechFrame2, 0);
          //
  
          //LP analysis
          LpCoeffs lpcs1 = LpcAnalyser.calcLPC(windowedSpeechFrame1, lpOrder);
          LpCoeffs lpcs2 = LpcAnalyser.calcLPC(windowedSpeechFrame2, lpOrder);
          //
  
          //LPC to LSF conversion
          double[] lsfs1 = LsfAnalyser.lpc2lsfInHz(lpcs1.getOneMinusA(), samplingRate);
          double[] lsfs2 = LsfAnalyser.lpc2lsfInHz(lpcs2.getOneMinusA(), samplingRate);
          //
  
          return getLsfDist(lsfs1, lsfs2, samplingRate);
      }
  
      public static double getLsfDist(double[] lsfs1, double[] lsfs2, int samplingRate) {
          double[] lsfWgt = getInverseHarmonicLSFWeights(lsfs1);
  
          return getLsfDist(lsfs1, lsfs2, samplingRate, lsfWgt);
      }
  
      public static double getLsfDist(double[] lsfs1, double[] lsfs2, int samplingRate, double[] lsfWgt) {
          assert lsfs1.length == lsfs2.length;
          assert lsfs1.length == lsfWgt.length;
  
          double dist = 0.0;
  
          for (int i = 0; i < lsfs1.length; i++)
              dist += 0.1 * lsfWgt[i] * Math.abs(lsfs1[i] - lsfs2[i]);
  
          dist = dist * 16000 / samplingRate;
          dist = Math.min(20.0, dist / lsfs1.length);
          dist = 10 * (dist + 1e-36);
  
          return dist;
      }
  
      public static double[] getInverseHarmonicLSFWeights(double[] lsfs) {
          assert lsfs.length > 1;
  
          int P = lsfs.length;
  
          int i;
  
          double[] lsfWgt = new double[P];
  
          lsfWgt[0] = 1.0 / Math.abs(lsfs[1] - lsfs[0]);
          lsfWgt[P - 1] = 0.5 / Math.abs(lsfs[P - 1] - lsfs[P - 2]);
  
          for (i = 1; i <= P - 2; i++)
              lsfWgt[i] = 1.0 / Math.min(Math.abs(lsfs[i] - lsfs[i - 1]), Math.abs(lsfs[i + 1] - lsfs[i]));
  
          //Emphasize low frequency LSFs 
          double tmp = 0.0;
          for (i = 0; i < P; i++) {
              lsfWgt[i] = Math.exp(-0.05 * i) * lsfWgt[i];
              tmp += lsfWgt[i];
          }
  
          //Normalize
          for (i = 0; i < P; i++)
              lsfWgt[i] /= tmp;
  
          //Compress dynamic range
          tmp = 0.0;
          for (i = 0; i <= P - 1; i++) {
              lsfWgt[i] = Math.sqrt(lsfWgt[i]);
              tmp += lsfWgt[i];
          }
  
          //Normalize
          for (i = 0; i < P; i++)
              lsfWgt[i] /= tmp;
  
          return lsfWgt;
      }
  
      public static double rmsLogSpectralDist(double[] speechFrame1, double[] speechFrame2, int fftSize, int lpOrder) {
          //Windowing
          double[] windowedSpeechFrame1 = new double[speechFrame1.length];
          System.arraycopy(speechFrame1, 0, windowedSpeechFrame1, 0, speechFrame1.length);
  
          double[] windowedSpeechFrame2 = new double[speechFrame2.length];
          System.arraycopy(speechFrame2, 0, windowedSpeechFrame2, 0, speechFrame2.length);
  
          HammingWindow w1 = new HammingWindow(speechFrame1.length);
          w1.apply(windowedSpeechFrame1, 0);
          HammingWindow w2 = new HammingWindow(speechFrame2.length);
          w2.apply(windowedSpeechFrame2, 0);
          //
  
          double[] Xabs1 = LpcAnalyser.calcSpecFrameLinear(speechFrame1, lpOrder, fftSize);
          double[] Xabs2 = LpcAnalyser.calcSpecFrameLinear(speechFrame2, lpOrder, fftSize);
  
          int w;
          double dist = 0.0;
          for (w = 0; w < Xabs1.length; w++)
              dist += 10 * Math.log10(Xabs1[w] * Xabs1[w]) - 10 * Math.log10(Xabs2[w] * Xabs2[w]);
  
          return Math.min(dist, MAX_SPECTRAL_DISTANCE);
      }
  
      public static double kullbackLeiblerSpectralDist(double[] speechFrame1, double[] speechFrame2, int fftSize, int lpOrder) {
          //Windowing
          double[] windowedSpeechFrame1 = new double[speechFrame1.length];
          System.arraycopy(speechFrame1, 0, windowedSpeechFrame1, 0, speechFrame1.length);
  
          double[] windowedSpeechFrame2 = new double[speechFrame2.length];
          System.arraycopy(speechFrame2, 0, windowedSpeechFrame2, 0, speechFrame2.length);
  
          HammingWindow w1 = new HammingWindow(speechFrame1.length);
          w1.apply(windowedSpeechFrame1, 0);
          HammingWindow w2 = new HammingWindow(speechFrame2.length);
          w2.apply(windowedSpeechFrame2, 0);
          //
  
          double[] Xabs1 = LpcAnalyser.calcSpecFrameLinear(speechFrame1, lpOrder, fftSize);
          double[] Xabs2 = LpcAnalyser.calcSpecFrameLinear(speechFrame2, lpOrder, fftSize);
  
          int w;
          double klDist = 0.0;
          for (w = 0; w < Xabs1.length; w++)
              klDist += Xabs1[w] * Math.log(Xabs1[w] / (Xabs2[w] + 1e-20) + 1e-20);
  
          return Math.min(klDist, MAX_SPECTRAL_DISTANCE);
      }
  
      public static double kullbackLeiblerSymmetricSpectralDist(double[] speechFrame1, double[] speechFrame2, int fftSize, int lpOrder) {
          //Windowing
          double[] windowedSpeechFrame1 = new double[speechFrame1.length];
          System.arraycopy(speechFrame1, 0, windowedSpeechFrame1, 0, speechFrame1.length);
  
          double[] windowedSpeechFrame2 = new double[speechFrame2.length];
          System.arraycopy(speechFrame2, 0, windowedSpeechFrame2, 0, speechFrame2.length);
  
          HammingWindow w1 = new HammingWindow(speechFrame1.length);
          w1.apply(windowedSpeechFrame1, 0);
          HammingWindow w2 = new HammingWindow(speechFrame2.length);
          w2.apply(windowedSpeechFrame2, 0);
          //
  
          double[] Xabs1 = LpcAnalyser.calcSpecFrameLinear(speechFrame1, lpOrder, fftSize);
          double[] Xabs2 = LpcAnalyser.calcSpecFrameLinear(speechFrame2, lpOrder, fftSize);
  
          int w;
  
          double klDist12 = 0.0;
          double klDist21 = 0.0;
          for (w = 0; w < Xabs1.length; w++)
              klDist12 += Xabs1[w] * Math.log(Xabs1[w] / (Xabs2[w] + 1e-20) + 1e-20);
          for (w = 0; w < Xabs2.length; w++)
              klDist21 += Xabs2[w] * Math.log(Xabs2[w] / (Xabs1[w] + 1e-20) + 1e-20);
  
          return Math.min(0.5 * (klDist12 + klDist21), MAX_SPECTRAL_DISTANCE);
      }
  
      public static double itakuraSaitoDistSymmetric(double[] speechFrame1, double[] speechFrame2, int fftSize, int lpOrder) {
          return 0.5 * itakuraSaitoDist(speechFrame1, speechFrame2, fftSize, lpOrder) + 0.5 * itakuraSaitoDist(speechFrame2, speechFrame1, fftSize, lpOrder);
      }
  
      public static double itakuraSaitoDist(double[] speechFrame1, double[] speechFrame2, int fftSize, int lpOrder) {
          double[] preemphasizedFrame1 = SignalProcUtils.applyPreemphasis(speechFrame1, 0.97);
          double[] preemphasizedFrame2 = SignalProcUtils.applyPreemphasis(speechFrame2, 0.97);
  
          //Windowing
          double[] windowedSpeechFrame1 = new double[preemphasizedFrame1.length];
          System.arraycopy(preemphasizedFrame1, 0, windowedSpeechFrame1, 0, preemphasizedFrame1.length);
  
          double[] windowedSpeechFrame2 = new double[preemphasizedFrame2.length];
          System.arraycopy(preemphasizedFrame2, 0, windowedSpeechFrame2, 0, preemphasizedFrame2.length);
  
          HammingWindow w1 = new HammingWindow(speechFrame1.length);
          w1.apply(windowedSpeechFrame1, 0);
          HammingWindow w2 = new HammingWindow(speechFrame2.length);
          w2.apply(windowedSpeechFrame2, 0);
          //
  
          int w;
          double[] Xabs1 = LpcAnalyser.calcSpecFrameLinear(speechFrame1, lpOrder, fftSize);
          double[] Xabs2 = LpcAnalyser.calcSpecFrameLinear(speechFrame2, lpOrder, fftSize);
  
          //Itakura-Saito distance using power spectrum: pf1/pf2 - log(pf1/pf2) - 1
          double dist = 0.0;
          for (w = 0; w < Xabs1.length; w++)
              Xabs1[w] = Xabs1[w] * Xabs1[w];
          for (w = 0; w < Xabs2.length; w++)
              Xabs2[w] = Xabs2[w] * Xabs2[w];
          for (w = 0; w < Xabs1.length; w++)
              dist += Xabs1[w] / Math.max(Xabs2[w], 1e-20) - Math.log(Math.max(Xabs1[w], 1e-20)) + Math.log(Math.max(Xabs2[w], 1e-20)) - 1.0;
  
          return Math.min(dist, MAX_SPECTRAL_DISTANCE);
      }
  
      //This is in fact the symmetrical version of Itakura-Saito distance
      public static double coshDist(double[] speechFrame1, double[] speechFrame2, int fftSize, int lpOrder) {
          double[] preemphasizedFrame1 = SignalProcUtils.applyPreemphasis(speechFrame1, 0.97);
          double[] preemphasizedFrame2 = SignalProcUtils.applyPreemphasis(speechFrame2, 0.97);
  
          //Windowing
          double[] windowedSpeechFrame1 = new double[preemphasizedFrame1.length];
          System.arraycopy(preemphasizedFrame1, 0, windowedSpeechFrame1, 0, preemphasizedFrame1.length);
  
          double[] windowedSpeechFrame2 = new double[preemphasizedFrame2.length];
          System.arraycopy(preemphasizedFrame2, 0, windowedSpeechFrame2, 0, preemphasizedFrame2.length);
  
          HammingWindow w1 = new HammingWindow(speechFrame1.length);
          w1.apply(windowedSpeechFrame1, 0);
          HammingWindow w2 = new HammingWindow(speechFrame2.length);
          w2.apply(windowedSpeechFrame2, 0);
          //
  
          int w;
          double[] Xabs1 = LpcAnalyser.calcSpecFrameLinear(speechFrame1, lpOrder, fftSize);
          double[] Xabs2 = LpcAnalyser.calcSpecFrameLinear(speechFrame2, lpOrder, fftSize);
  
          //COSH distance using power spectrum: 
          // dcosh = 0.5(d12+d21) where
          // dij = Xi/Xj - log(Xi/Xj) - 1 and Xi, Xj are the power spectra under comparison
          double dist12 = 0.0;
          double dist21 = 0.0;
          for (w = 0; w < Xabs1.length; w++)
              Xabs1[w] = Xabs1[w] * Xabs1[w];
          for (w = 0; w < Xabs2.length; w++)
              Xabs2[w] = Xabs2[w] * Xabs2[w];
          for (w = 0; w < Xabs1.length; w++)
              dist12 += Xabs1[w] / (Xabs2[w] + 1e-20) - Math.log(Xabs1[w] / (Xabs2[w] + 1e-20) + 1e-20) - 1;
          for (w = 0; w < Xabs2.length; w++)
              dist21 += Xabs2[w] / (Xabs1[w] + 1e-20) - Math.log(Xabs2[w] / (Xabs1[w] + 1e-20) + 1e-20) - 1;
  
          return Math.min(0.5 * (dist12 + dist21), MAX_SPECTRAL_DISTANCE);
      }
  
      public static void main(String[] args) {
          AudioInputStream inputAudio = null;
          try {
              inputAudio = AudioSystem.getAudioInputStream(new File(args[0]));
          } catch (UnsupportedAudioFileException e) {
              // TODO Auto-generated catch block
              e.printStackTrace();
          } catch (IOException e) {
              // TODO Auto-generated catch block
              e.printStackTrace();
          }
  
          if (inputAudio != null) {
              int samplingRate = (int) inputAudio.getFormat().getSampleRate();
              int ws = (int) (samplingRate * 0.020);
              int ss = (int) (samplingRate * 0.010);
  
              AudioDoubleDataSource signal = new AudioDoubleDataSource(inputAudio);
              double[] x = signal.getAllData();
  
              int numfrm = (int) ((x.length - (double) ws) / ss - 2.0);
              double[] frm1 = new double[ws];
              double[] frm2 = new double[ws];
              double lsfDist;
  
              for (int i = 0; i < numfrm; i++) {
                  System.arraycopy(x, i * ss, frm1, 0, ws);
                  System.arraycopy(x, (i + 1) * ss, frm2, 0, ws);
  
                  lsfDist = SpectralDistanceMeasures.lsfDist(frm1, frm2, samplingRate);
                  System.out.println("Distance between frame " + String.valueOf(i + 1) + " and frame " + String.valueOf(i + 2) + " = " + String.valueOf(lsfDist));
              }
          }
      }
  }
  
  