Added support for producing approximate mean opinion scores.

CMakeLists.txt

@@ -87,6 +87,8 @@ add_library(zimtohrli_base STATIC
     cpp/zimt/loudness.h
     cpp/zimt/masking.cc
     cpp/zimt/masking.h
+    cpp/zimt/mos.cc
+    cpp/zimt/mos.h
     cpp/zimt/zimtohrli.cc
     cpp/zimt/zimtohrli.h
 )
@@ -175,6 +177,7 @@ add_executable(zimtohrli_test
     cpp/zimt/filterbank_test.cc
     cpp/zimt/loudness_test.cc
     cpp/zimt/masking_test.cc
+    cpp/zimt/mos_test.cc
     cpp/zimt/zimtohrli_test.cc
     cpp/zimt/test_file_paths.cc
 )

cpp/zimt/compare.cc

@@ -55,6 +55,7 @@
 #include "sndfile.h"
 #include "zimt/audio.h"
 #include "zimt/cam.h"
+#include "zimt/mos.h"
 #include "zimt/ux.h"
 #include "zimt/zimtohrli.h"
 
@@ -79,6 +80,12 @@ ABSL_FLAG(float, time_norm_order, zimtohrli::Zimtohrli{}.time_norm_order,
 ABSL_FLAG(bool, normalize_amplitude, true,
           "whether to normalize the amplitude of all B sounds to the same max "
           "amplitude as the A sound");
+ABSL_FLAG(bool, output_zimtohrli_distance, false,
+          "Whether to output the raw Zimtohrli distance instead of a mapped "
+          "mean opinion score.");
+ABSL_FLAG(bool, per_channel, false,
+          "Whether to output the produced metric per channel instead of a "
+          "single value for all channels.");
 
 namespace zimtohrli {
 
@@ -172,6 +179,13 @@ std::ostream& operator<<(std::ostream& outs, const DistanceData& data) {
   return outs;
 }
 
+float GetMetric(float zimtohrli_score) {
+  if (absl::GetFlag(FLAGS_output_zimtohrli_distance)) {
+    return zimtohrli_score;
+  }
+  return MOSFromZimtohrli(zimtohrli_score);
+}
+
 int Main(int argc, char* argv[]) {
   absl::ParseCommandLine(argc, argv);
   const std::string path_a = absl::GetFlag(FLAGS_path_a);
@@ -291,6 +305,7 @@ int Main(int argc, char* argv[]) {
   }
 
   const bool ux = absl::GetFlag(FLAGS_ux);
+  const bool per_channel = absl::GetFlag(FLAGS_per_channel);
   if (!ux && !verbose) {
     const size_t num_downscaled_samples_a = static_cast<size_t>(
         std::ceil(static_cast<float>(file_a->Frames().shape()[1]) *
@@ -301,32 +316,53 @@ int Main(int argc, char* argv[]) {
         {num_downscaled_samples_a, z.cam_filterbank->filter.Size()});
     hwy::AlignedNDArray<float, 2> partial_energy_channels_db_a(
         {num_downscaled_samples_a, z.cam_filterbank->filter.Size()});
-    hwy::AlignedNDArray<float, 2> spectrogram_a(
-        {num_downscaled_samples_a, z.cam_filterbank->filter.Size()});
+    std::vector<hwy::AlignedNDArray<float, 2>> file_a_spectrograms;
     for (size_t channel_index = 0; channel_index < file_a->Info().channels;
          ++channel_index) {
+      hwy::AlignedNDArray<float, 2> spectrogram(
+          {num_downscaled_samples_a, z.cam_filterbank->filter.Size()});
       z.Spectrogram(file_a->Frames()[{channel_index}], channels_a,
                     energy_channels_db_a, partial_energy_channels_db_a,
-                    spectrogram_a);
-      for (const AudioFile& file_b : file_b_vector) {
-        const size_t num_downscaled_samples_b = static_cast<size_t>(std::ceil(
-            static_cast<float>(file_b.Frames().shape()[1]) *
-            time_resolution_frequency / z.cam_filterbank->sample_rate));
-        hwy::AlignedNDArray<float, 2> channels_b(
-            {file_b.Frames().shape()[1], z.cam_filterbank->filter.Size()});
-        hwy::AlignedNDArray<float, 2> energy_channels_db_b(
-            {num_downscaled_samples_b, z.cam_filterbank->filter.Size()});
-        hwy::AlignedNDArray<float, 2> partial_energy_channels_db_b(
-            {num_downscaled_samples_b, z.cam_filterbank->filter.Size()});
-        hwy::AlignedNDArray<float, 2> spectrogram_b(
-            {num_downscaled_samples_b, z.cam_filterbank->filter.Size()});
+                    spectrogram);
+      file_a_spectrograms.push_back(std::move(spectrogram));
+    }
+    for (int file_b_index = 0; file_b_index < file_b_vector.size();
+         ++file_b_index) {
+      const AudioFile& file_b = file_b_vector[file_b_index];
+      const size_t num_downscaled_samples_b = static_cast<size_t>(
+          std::ceil(static_cast<float>(file_b.Frames().shape()[1]) *
+                    time_resolution_frequency / z.cam_filterbank->sample_rate));
+      hwy::AlignedNDArray<float, 2> channels_b(
+          {file_b.Frames().shape()[1], z.cam_filterbank->filter.Size()});
+      hwy::AlignedNDArray<float, 2> energy_channels_db_b(
+          {num_downscaled_samples_b, z.cam_filterbank->filter.Size()});
+      hwy::AlignedNDArray<float, 2> partial_energy_channels_db_b(
+          {num_downscaled_samples_b, z.cam_filterbank->filter.Size()});
+      hwy::AlignedNDArray<float, 2> spectrogram_b(
+          {num_downscaled_samples_b, z.cam_filterbank->filter.Size()});
+      float sum_of_squares = 0;
+      for (size_t channel_index = 0; channel_index < file_a->Info().channels;
+           ++channel_index) {
         z.Spectrogram(file_b.Frames()[{channel_index}], channels_b,
                       energy_channels_db_b, partial_energy_channels_db_b,
                       spectrogram_b);
-        std::cout << z.Distance(false, spectrogram_a, spectrogram_b,
-                                unwarp_window_samples)
-                         .value
-                  << std::endl;
+        const float distance =
+            z.Distance(false, file_a_spectrograms[channel_index], spectrogram_b,
+                       unwarp_window_samples)
+                .value;
+        if (per_channel) {
+          std::cout << GetMetric(distance) << std::endl;
+        } else {
+          sum_of_squares += distance * distance;
+        }
+      }
+      if (!per_channel) {
+        for (int file_b_index = 0; file_b_index < file_b_vector.size();
+             ++file_b_index) {
+          std::cout << GetMetric(std::sqrt(sum_of_squares /
+                                           float(file_a->Info().channels)))
+                    << std::endl;
+        }
       }
     }
     return 0;
@@ -358,6 +394,7 @@ int Main(int argc, char* argv[]) {
       const AudioFile& file_b = file_b_vector[b_index];
       std::cout << "A (" << file_a->Path() << ") vs B (" << file_b.Path() << ")"
                 << std::endl;
+      float sum_of_squares = 0;
       for (size_t channel_index = 0;
            channel_index < comparison.analysis_a.size(); ++channel_index) {
         std::cout << "  Channel " << channel_index << std::endl;
@@ -385,7 +422,18 @@ int Main(int argc, char* argv[]) {
             time_resolution_frequency, unwarp_window_samples);
         std::cout << "    Phons channel distance: " << phons_channel_distance
                   << std::endl;
+
+        const float distance = phons_channel_distance.distance.value;
+        sum_of_squares += distance * distance;
+
+        std::cout << "    Channel MOS: " << MOSFromZimtohrli(distance)
+                  << std::endl;
       }
+      const float zimtohrli_file_distance =
+          std::sqrt(sum_of_squares / float(comparison.analysis_a.size()));
+      std::cout << "  File distance: " << zimtohrli_file_distance << std::endl;
+      std::cout << "  File MOS: " << MOSFromZimtohrli(zimtohrli_file_distance)
+                << std::endl;
     }
     return 0;
   }

cpp/zimt/goohrli.cc

@@ -21,6 +21,7 @@
 #include "hwy/aligned_allocator.h"
 #include "hwy/base.h"
 #include "zimt/cam.h"
+#include "zimt/mos.h"
 #include "zimt/zimtohrli.h"
 
 EnergyAndMaxAbsAmplitude Measure(const float* signal, int size) {
@@ -33,17 +34,22 @@ EnergyAndMaxAbsAmplitude Measure(const float* signal, int size) {
       .MaxAbsAmplitude = measurements.max_abs_amplitude};
 }
 
-EnergyAndMaxAbsAmplitude NormalizeAmplitudes(float max_abs_amplitude,
-                                             float* signal, int size) {
+EnergyAndMaxAbsAmplitude NormalizeAmplitude(float max_abs_amplitude,
+                                            float* signal, int size) {
   hwy::AlignedNDArray<float, 1> signal_array({static_cast<size_t>(size)});
   hwy::CopyBytes(signal, signal_array.data(), size * sizeof(float));
   const zimtohrli::EnergyAndMaxAbsAmplitude measurements =
       zimtohrli::NormalizeAmplitude(max_abs_amplitude, signal_array[{}]);
+  hwy::CopyBytes(signal_array.data(), signal, size * sizeof(float));
   return EnergyAndMaxAbsAmplitude{
       .EnergyDBFS = measurements.energy_db_fs,
       .MaxAbsAmplitude = measurements.max_abs_amplitude};
 }
 
+float MOSFromZimtohrli(float zimtohrli_distance) {
+  return zimtohrli::MOSFromZimtohrli(zimtohrli_distance);
+}
+
 Zimtohrli CreateZimtohrli(float sample_rate, float frequency_resolution) {
   zimtohrli::Cam cam{.minimum_bandwidth_hz = frequency_resolution};
   cam.high_threshold_hz = std::min(cam.high_threshold_hz, sample_rate);

cpp/zimt/mos.cc

@@ -0,0 +1,37 @@
+// Copyright 2024 The Zimtohrli Authors. All Rights Reserved.
+//
+// Licensed under the Apache License, Version 2.0 (the "License");
+// you may not use this file except in compliance with the License.
+// You may obtain a copy of the License at
+//
+//      http://www.apache.org/licenses/LICENSE-2.0
+//
+// Unless required by applicable law or agreed to in writing, software
+// distributed under the License is distributed on an "AS IS" BASIS,
+// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+// See the License for the specific language governing permissions and
+// limitations under the License.
+
+#include "zimt/mos.h"
+
+#include <array>
+#include <cmath>
+
+namespace zimtohrli {
+
+const std::array<float, 4> params = {3.439e+00, -4.138e-02, 3.008e+00,
+                                     -1.354e-01};
+
+namespace {
+
+float sigmoid(float x) { return 1 / (1 + std::exp(-x)); }
+
+}  // namespace
+
+// Optimized using `mos_mapping.ipynb`.
+float MOSFromZimtohrli(float zimtohrli_distance) {
+  return 1 + 2 * (sigmoid(params[0] + params[1] * zimtohrli_distance) +
+                  sigmoid(params[2] + params[3] * zimtohrli_distance));
+}
+
+}  // namespace zimtohrli

cpp/zimt/mos.h

@@ -0,0 +1,30 @@
+// Copyright 2024 The Zimtohrli Authors. All Rights Reserved.
+//
+// Licensed under the Apache License, Version 2.0 (the "License");
+// you may not use this file except in compliance with the License.
+// You may obtain a copy of the License at
+//
+//      http://www.apache.org/licenses/LICENSE-2.0
+//
+// Unless required by applicable law or agreed to in writing, software
+// distributed under the License is distributed on an "AS IS" BASIS,
+// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+// See the License for the specific language governing permissions and
+// limitations under the License.
+
+#ifndef CPP_ZIMT_MOS_H_
+#define CPP_ZIMT_MOS_H_
+
+namespace zimtohrli {
+
+// Returns a _very_approximate_ mean opinion score based on the
+// provided Zimtohrli distance.
+// This is calibrated using default settings of v0.1.5, with a
+// minimum channel bandwidth (zimtohrli::Cam.minimum_bandwidth_hz)
+// of 5Hz and perceptual sample rate
+// (zimtohrli::Distance(..., perceptual_sample_rate, ...) of 100Hz.
+float MOSFromZimtohrli(float zimtohrli_distance);
+
+}  // namespace zimtohrli
+
+#endif  // CPP_ZIMT_MOS_H_

cpp/zimt/mos_test.cc

@@ -0,0 +1,34 @@
+// Copyright 2024 The Zimtohrli Authors. All Rights Reserved.
+//
+// Licensed under the Apache License, Version 2.0 (the "License");
+// you may not use this file except in compliance with the License.
+// You may obtain a copy of the License at
+//
+//      http://www.apache.org/licenses/LICENSE-2.0
+//
+// Unless required by applicable law or agreed to in writing, software
+// distributed under the License is distributed on an "AS IS" BASIS,
+// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+// See the License for the specific language governing permissions and
+// limitations under the License.
+
+#include "zimt/mos.h"
+
+#include "gtest/gtest.h"
+
+namespace zimtohrli {
+
+namespace {
+
+TEST(MOS, MOSFromZimtohrli) {
+  const std::vector<float> zimt_scores = {5, 20, 40, 80};
+  const std::vector<float> mos = {4.746790024702545, 4.01181593706087,
+                                  2.8773086764995064, 2.0648331964917945};
+  for (size_t index = 0; index < zimt_scores.size(); ++index) {
+    ASSERT_NEAR(MOSFromZimtohrli(zimt_scores[index]), mos[index], 1e-2);
+  }
+}
+
+}  // namespace
+
+}  // namespace zimtohrli

go/bin/compare/compare.go

@@ -28,6 +28,8 @@ import (
 func main() {
 	pathA := flag.String("path_a", "", "Path to ffmpeg-decodable file with signal A.")
 	pathB := flag.String("path_b", "", "Path to ffmpeg-decodable file with signal B.")
+	outputZimtohrliDistance := flag.Bool("output_zimtohrli_distance", false, "Whether to output the raw Zimtohrli distance instead of a mapped mean opinion score.")
+	perChannel := flag.Bool("per_channel", false, "Whether to output the produced metric per channel instead of a single value for all channels.")
 	frequencyResolution := flag.Float64("frequency_resolution", 5.0, "Band width of smallest filter, i.e. expected frequency resolution of human hearing.")
 	flag.Parse()
 
@@ -53,10 +55,25 @@ func main() {
 		log.Panic(fmt.Errorf("%q has %v channels, and %q has %v channels", *pathA, len(signalA.Samples), *pathB, len(signalB.Samples)))
 	}
 
+	getMetric := func(f float32) float32 {
+		if *outputZimtohrliDistance {
+			return f
+		}
+		return goohrli.MOSFromZimtohrli(f)
+	}
+
 	g := goohrli.New(signalA.Rate, *frequencyResolution)
-	for channelIndex := range signalA.Samples {
-		measurement := goohrli.Measure(signalA.Samples[channelIndex])
-		goohrli.NormalizeAmplitude(measurement.MaxAbsAmplitude, signalB.Samples[channelIndex])
-		fmt.Println(g.Distance(signalA.Samples[channelIndex], signalB.Samples[channelIndex]))
+	if *perChannel {
+		for channelIndex := range signalA.Samples {
+			measurement := goohrli.Measure(signalA.Samples[channelIndex])
+			goohrli.NormalizeAmplitude(measurement.MaxAbsAmplitude, signalB.Samples[channelIndex])
+			fmt.Println(getMetric(g.Distance(signalA.Samples[channelIndex], signalB.Samples[channelIndex])))
+		}
+	} else {
+		dist, err := g.NormalizedAudioDistance(signalA, signalB)
+		if err != nil {
+			log.Panic(err)
+		}
+		fmt.Println(getMetric(float32(dist)))
 	}
 }

go/goohrli/goohrli.go

@@ -48,13 +48,18 @@ func Measure(signal []float32) EnergyAndMaxAbsAmplitude {
 // NormalizeAmplitude normalizes the amplitudes of the signal so that it has the provided max
 // amplitude, and returns the new energ in dB FS, and the new maximum absolute amplitude.
 func NormalizeAmplitude(maxAbsAmplitude float32, signal []float32) EnergyAndMaxAbsAmplitude {
-	measurements := C.NormalizeAmplitudes(C.float(maxAbsAmplitude), (*C.float)(&signal[0]), C.int(len(signal)))
+	measurements := C.NormalizeAmplitude(C.float(maxAbsAmplitude), (*C.float)(&signal[0]), C.int(len(signal)))
 	return EnergyAndMaxAbsAmplitude{
 		EnergyDBFS:      float32(measurements.EnergyDBFS),
 		MaxAbsAmplitude: float32(measurements.MaxAbsAmplitude),
 	}
 }
 
+// MOSFromZimtohrli returns an approximate mean opinion score for a given zimtohrli distance.
+func MOSFromZimtohrli(zimtohrliDistance float32) float32 {
+	return float32(C.MOSFromZimtohrli(C.float(zimtohrliDistance)))
+}
+
 // Goohrli is a Go wrapper around zimtohrli::Zimtohrli.
 type Goohrli struct {
 	zimtohrli           C.Zimtohrli

go/goohrli/goohrli.h

@@ -56,8 +56,11 @@ EnergyAndMaxAbsAmplitude Measure(const float* signal, int size);
 // Normalizes the amplitudes of the signal so that it has the provided max
 // amplitude, and returns the new energ in dB FS, and the new maximum absolute
 // amplitude.
-EnergyAndMaxAbsAmplitude NormalizeAmplitudes(float max_abs_amplitude,
-                                             float* signal_data, int size);
+EnergyAndMaxAbsAmplitude NormalizeAmplitude(float max_abs_amplitude,
+                                            float* signal_data, int size);
+
+// Returns an approximate MOS score for a given Zimtohrli distance.
+float MOSFromZimtohrli(float zimtohrli_distance);
 
 // Deletes a zimtohrli::Analysis.
 void FreeAnalysis(Analysis a);