diff options
Diffstat (limited to 'src/3rdparty/resonance-audio/resonance_audio/ambisonics')
24 files changed, 3932 insertions, 0 deletions
diff --git a/src/3rdparty/resonance-audio/resonance_audio/ambisonics/ambisonic_binaural_decoder.cc b/src/3rdparty/resonance-audio/resonance_audio/ambisonics/ambisonic_binaural_decoder.cc new file mode 100644 index 000000000..2a278389b --- /dev/null +++ b/src/3rdparty/resonance-audio/resonance_audio/ambisonics/ambisonic_binaural_decoder.cc @@ -0,0 +1,80 @@ +/* +Copyright 2018 Google Inc. 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 "ambisonics/ambisonic_binaural_decoder.h" + +#include "ambisonics/utils.h" +#include "base/constants_and_types.h" + + +namespace vraudio { + +AmbisonicBinauralDecoder::AmbisonicBinauralDecoder(const AudioBuffer& sh_hrirs, + size_t frames_per_buffer, + FftManager* fft_manager) + : fft_manager_(fft_manager), + freq_input_(kNumMonoChannels, NextPowTwo(frames_per_buffer) * 2), + filtered_input_(kNumMonoChannels, frames_per_buffer) { + CHECK(fft_manager_); + CHECK_NE(frames_per_buffer, 0U); + const size_t num_channels = sh_hrirs.num_channels(); + const size_t filter_size = sh_hrirs.num_frames(); + CHECK_NE(num_channels, 0U); + CHECK_NE(filter_size, 0U); + sh_hrir_filters_.reserve(num_channels); + for (size_t i = 0; i < num_channels; ++i) { + sh_hrir_filters_.emplace_back( + new PartitionedFftFilter(filter_size, frames_per_buffer, fft_manager_)); + sh_hrir_filters_[i]->SetTimeDomainKernel(sh_hrirs[i]); + } +} + +void AmbisonicBinauralDecoder::Process(const AudioBuffer& input, + AudioBuffer* output) { + + DCHECK(output); + DCHECK_EQ(kNumStereoChannels, output->num_channels()); + DCHECK_EQ(input.num_frames(), output->num_frames()); + DCHECK_EQ(input.num_channels(), sh_hrir_filters_.size()); + + output->Clear(); + + AudioBuffer::Channel* freq_input_channel = &freq_input_[0]; + AudioBuffer::Channel* filtered_input_channel = &filtered_input_[0]; + AudioBuffer::Channel* output_channel_0 = &(*output)[0]; + AudioBuffer::Channel* output_channel_1 = &(*output)[1]; + for (size_t channel = 0; channel < input.num_channels(); ++channel) { + const int degree = GetPeriphonicAmbisonicDegreeForChannel(channel); + fft_manager_->FreqFromTimeDomain(input[channel], freq_input_channel); + sh_hrir_filters_[channel]->Filter(*freq_input_channel); + sh_hrir_filters_[channel]->GetFilteredSignal(filtered_input_channel); + if (degree < 0) { + // Degree is negative: spherical harmonic is asymetric. + // So add contributions to the left channel and subtract from the right + // channel. + *output_channel_0 += *filtered_input_channel; + *output_channel_1 -= *filtered_input_channel; + + } else { + // Degree is zero or positive: spherical harmonic is symetric. + // So add contributions to both left and right channels. + *output_channel_0 += *filtered_input_channel; + *output_channel_1 += *filtered_input_channel; + } + } +} + +} // namespace vraudio diff --git a/src/3rdparty/resonance-audio/resonance_audio/ambisonics/ambisonic_binaural_decoder.h b/src/3rdparty/resonance-audio/resonance_audio/ambisonics/ambisonic_binaural_decoder.h new file mode 100644 index 000000000..98a861234 --- /dev/null +++ b/src/3rdparty/resonance-audio/resonance_audio/ambisonics/ambisonic_binaural_decoder.h @@ -0,0 +1,71 @@ +/* +Copyright 2018 Google Inc. 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 RESONANCE_AUDIO_AMBISONICS_AMBISONIC_BINAURAL_DECODER_H_ +#define RESONANCE_AUDIO_AMBISONICS_AMBISONIC_BINAURAL_DECODER_H_ + +#include <vector> + +#include "base/audio_buffer.h" +#include "dsp/fft_manager.h" +#include "dsp/partitioned_fft_filter.h" + +namespace vraudio { + +// Decodes an Ambisonic sound field, of an arbitrary order, to binaural stereo, +// by performing convolution in the spherical harmonics domain. The order (hence +// the number of channels) of the input must match the order (hence the number +// of channels) of the spherical harmonic-encoded Head Related Impulse Responses +// (HRIRs). Assumes that HRIRs are symmetric with respect to the sagittal plane. +class AmbisonicBinauralDecoder { + public: + // Constructs an |AmbisonicBinauralDecoder| from an |AudioBuffer| containing + // spherical harmonic representation of HRIRs. The order of spherical + // harmonic-encoded HRIRs (hence the number of channels) must match the order + // of the Ambisonic sound field input. + // + // @param sh_hrirs |AudioBuffer| containing time-domain spherical harmonic + // encoded symmetric HRIRs. + // @param frames_per_buffer Number of frames in each input/output buffer. + // @param fft_manager Pointer to a manager to perform FFT transformations. + AmbisonicBinauralDecoder(const AudioBuffer& sh_hrirs, + size_t frames_per_buffer, FftManager* fft_manager); + + // Processes an Ambisonic sound field input and outputs a binaurally decoded + // stereo buffer. + // + // @param input Input buffer to be processed. + // @param output Pointer to a stereo output buffer. + void Process(const AudioBuffer& input, AudioBuffer* output); + + private: + // Manager for all FFT related functionality (not owned). + FftManager* const fft_manager_; + + // Spherical Harmonic HRIR filter kernels. + std::vector<std::unique_ptr<PartitionedFftFilter>> sh_hrir_filters_; + + // Frequency domain representation of the input signal. + PartitionedFftFilter::FreqDomainBuffer freq_input_; + + // Temporary audio buffer to store the convolution output for asymmetric or + // symmetric spherical harmonic HRIR. + AudioBuffer filtered_input_; +}; + +} // namespace vraudio + +#endif // RESONANCE_AUDIO_AMBISONICS_AMBISONIC_BINAURAL_DECODER_H_ diff --git a/src/3rdparty/resonance-audio/resonance_audio/ambisonics/ambisonic_binaural_decoder_test.cc b/src/3rdparty/resonance-audio/resonance_audio/ambisonics/ambisonic_binaural_decoder_test.cc new file mode 100644 index 000000000..ca2b0dfca --- /dev/null +++ b/src/3rdparty/resonance-audio/resonance_audio/ambisonics/ambisonic_binaural_decoder_test.cc @@ -0,0 +1,162 @@ +/* +Copyright 2018 Google Inc. 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 "ambisonics/ambisonic_binaural_decoder.h" + +#include "third_party/googletest/googletest/include/gtest/gtest.h" +#include "base/audio_buffer.h" +#include "base/constants_and_types.h" +#include "dsp/fft_manager.h" + +namespace vraudio { + +namespace { + +// Generates sample audio data where the first sample is 0 and each consecutive +// sample is incremented by 0.001. Then it moves to the next channel and +// continues to increment the samples. +std::vector<std::vector<float>> GenerateAudioData(size_t num_channels, + size_t num_samples) { + std::vector<std::vector<float>> input_data(num_channels, + std::vector<float>(num_samples)); + float sample_value = 0.0f; + const float increment = 0.001f; + for (size_t channel = 0; channel < num_channels; ++channel) { + for (size_t sample = 0; sample < num_samples; ++sample) { + input_data[channel][sample] = sample_value; + sample_value += increment; + } + } + return input_data; +} + +// Number of frames in each audio input/output buffer. +const size_t kFramesPerBuffer = 18; + +} // namespace + +// Tests whether binaural docoding of Ambisonic input using short HRIR filters +// (shorter than the number of frames per buffer) gives correct results. +TEST(AmbisonicBinauralDecoderTest, ShortFilterTest) { + const std::vector<std::vector<float>> kInputData = + GenerateAudioData(kNumFirstOrderAmbisonicChannels, kFramesPerBuffer); + const std::vector<float> kExpectedOutputLeft = { + 0.0042840000f, 0.0087780003f, 0.013486000f, 0.018412000f, 0.023560001f, + 0.028934000f, 0.034538001f, 0.040376000f, 0.046452001f, 0.052770000f, + 0.059333999f, 0.066147998f, 0.073215999f, 0.080541998f, 0.088129997f, + 0.095983997f, 0.10410800f, 0.10638600f}; + const std::vector<float> kExpectedOutputRight = { + 0.0036720000f, 0.0074840002f, 0.011438000f, 0.015536000f, 0.019780001f, + 0.024172001f, 0.028713999f, 0.033408001f, 0.038256001f, 0.043260001f, + 0.048422001f, 0.053743999f, 0.059227999f, 0.064875998f, 0.070689999f, + 0.076672003f, 0.082823999f, 0.084252000f}; + const std::vector<std::vector<float>> kHrirDataShort = + GenerateAudioData(kNumFirstOrderAmbisonicChannels, kFramesPerBuffer - 1); + + AudioBuffer sh_hrirs(kHrirDataShort.size(), kHrirDataShort[0].size()); + sh_hrirs = kHrirDataShort; + + AudioBuffer input(kInputData.size(), kInputData[0].size()); + input = kInputData; + + AudioBuffer output(kNumStereoChannels, kInputData[0].size()); + output.Clear(); + + FftManager fft_manager(kFramesPerBuffer); + AmbisonicBinauralDecoder decoder(sh_hrirs, kFramesPerBuffer, &fft_manager); + decoder.Process(input, &output); + + for (size_t sample = 0; sample < kFramesPerBuffer; ++sample) { + EXPECT_NEAR(kExpectedOutputLeft[sample], output[0][sample], kEpsilonFloat); + EXPECT_NEAR(kExpectedOutputRight[sample], output[1][sample], kEpsilonFloat); + } +} + +// Tests whether binaural docoding of Ambisonic input using HRIR filters +// with the same size as frames per buffer gives correct results. +TEST(AmbisonicBinauralDecoderTest, SameSizeFilterTest) { + const std::vector<std::vector<float>> kInputData = + GenerateAudioData(kNumFirstOrderAmbisonicChannels, kFramesPerBuffer); + const std::vector<float> kExpectedOutputLeft = { + 0.0045360001f, 0.0092879999f, 0.014260001f, 0.019455999f, 0.024879999f, + 0.030536000f, 0.036428001f, 0.042560000f, 0.048935998f, 0.055560000f, + 0.062436000f, 0.069568001f, 0.076959997f, 0.084615998f, 0.092540003f, + 0.10073600f, 0.10920800f, 0.11796000f}; + const std::vector<float> kExpectedOutputRight = { + 0.0038880000f, 0.0079199998f, 0.012098000f, 0.016424000f, 0.020900000f, + 0.025528001f, 0.030309999f, 0.035248000f, 0.040344000f, 0.045600001f, + 0.051018000f, 0.056600001f, 0.062348001f, 0.068264000f, 0.074349999f, + 0.080608003f, 0.087040000f, 0.093648002f}; + const std::vector<std::vector<float>> kHrirDataSame = + GenerateAudioData(kNumFirstOrderAmbisonicChannels, kFramesPerBuffer); + + AudioBuffer sh_hrirs(kHrirDataSame.size(), kHrirDataSame[0].size()); + sh_hrirs = kHrirDataSame; + + AudioBuffer input(kInputData.size(), kInputData[0].size()); + input = kInputData; + + AudioBuffer output(kNumStereoChannels, kInputData[0].size()); + output.Clear(); + + FftManager fft_manager(kFramesPerBuffer); + AmbisonicBinauralDecoder decoder(sh_hrirs, kFramesPerBuffer, &fft_manager); + decoder.Process(input, &output); + + for (size_t sample = 0; sample < kFramesPerBuffer; ++sample) { + EXPECT_NEAR(kExpectedOutputLeft[sample], output[0][sample], kEpsilonFloat); + EXPECT_NEAR(kExpectedOutputRight[sample], output[1][sample], kEpsilonFloat); + } +} + +// Tests whether binaural docoding of Ambisonic input using long HRIR filters +// (longer than the number of frames per buffer) gives correct results. +TEST(AmbisonicBinauralDecoderTest, LongSizeFilterTest) { + const std::vector<std::vector<float>> kInputData = + GenerateAudioData(kNumFirstOrderAmbisonicChannels, kFramesPerBuffer); + const std::vector<float> kExpectedOutputLeft = { + 0.0047880001f, 0.0097980006f, 0.015034000f, 0.020500001f, 0.026200000f, + 0.032138001f, 0.038318001f, 0.044744000f, 0.051419999f, 0.058350001f, + 0.065537997f, 0.072987996f, 0.080704004f, 0.088689998f, 0.096950002f, + 0.10548800f, 0.11430800f, 0.12341400f}; + const std::vector<float> kExpectedOutputRight = { + 00.0041040001f, 0.0083560003f, 0.012758000f, 0.017312000f, 0.022020001f, + 0.026883999f, 0.031906001f, 0.037087999f, 0.042431999f, 0.047940001f, + 0.053613998f, 0.059455998f, 0.065467998f, 0.071652003f, 0.078010000f, + 0.084543996f, 0.091256000f, 0.098148003f}; + const std::vector<std::vector<float>> kHrirDataLong = + GenerateAudioData(kNumFirstOrderAmbisonicChannels, kFramesPerBuffer + 1); + + AudioBuffer sh_hrirs(kHrirDataLong.size(), kHrirDataLong[0].size()); + sh_hrirs = kHrirDataLong; + + AudioBuffer input(kInputData.size(), kInputData[0].size()); + input = kInputData; + + AudioBuffer output(kNumStereoChannels, kInputData[0].size()); + output.Clear(); + + FftManager fft_manager(kFramesPerBuffer); + AmbisonicBinauralDecoder decoder(sh_hrirs, kFramesPerBuffer, &fft_manager); + decoder.Process(input, &output); + + for (size_t sample = 0; sample < kFramesPerBuffer; ++sample) { + EXPECT_NEAR(kExpectedOutputLeft[sample], output[0][sample], kEpsilonFloat); + EXPECT_NEAR(kExpectedOutputRight[sample], output[1][sample], kEpsilonFloat); + } +} + +} // namespace vraudio diff --git a/src/3rdparty/resonance-audio/resonance_audio/ambisonics/ambisonic_codec.h b/src/3rdparty/resonance-audio/resonance_audio/ambisonics/ambisonic_codec.h new file mode 100644 index 000000000..546a88470 --- /dev/null +++ b/src/3rdparty/resonance-audio/resonance_audio/ambisonics/ambisonic_codec.h @@ -0,0 +1,62 @@ +/* +Copyright 2018 Google Inc. 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 RESONANCE_AUDIO_AMBISONICS_AMBISONIC_CODEC_H_ +#define RESONANCE_AUDIO_AMBISONICS_AMBISONIC_CODEC_H_ + +#include <vector> + +#include "base/audio_buffer.h" +#include "base/spherical_angle.h" + +namespace vraudio { +// An encoder/decoder for ambisonic sound fields. It supports variable ambisonic +// order, ACN channel sequencing and SN3D normalization. +class AmbisonicCodec { + public: + virtual ~AmbisonicCodec() {} + // Encodes the given buffer of decoded vectors (unencoded data). + // + // @param input |AudioBuffer| of decoded (unencoded) data. + // @param output |AudioBuffer| to store the encoded data. + virtual void EncodeBuffer(const AudioBuffer& input, AudioBuffer* output) = 0; + + // Decodes the given |AudioBuffer| of planar (non-interleaved) encoded data. + // + // @param input |AudioBuffer| of encoded data. + // @param output |AudioBuffer| to store the decoded data. + virtual void DecodeBuffer(const AudioBuffer& input, AudioBuffer* output) = 0; + + // Returns the maximum periphonic ambisonic order that this codec supports. + virtual int ambisonic_order() const = 0; + + // Returns the number of angles for this codec. + virtual size_t num_angles() const = 0; + + // Returns the maximum number of spherical harmonics computed by this codec. + virtual size_t num_spherical_harmonics() const = 0; + + // Returns the spherical angles currently used to define the encoder/decoder + // matrices. + virtual const std::vector<SphericalAngle>& angles() const = 0; + + // Sets the spherical angles used to build the encoder/decoder matrices. + virtual void set_angles(const std::vector<SphericalAngle>& angles) = 0; +}; + +} // namespace vraudio + +#endif // RESONANCE_AUDIO_AMBISONICS_AMBISONIC_CODEC_H_ diff --git a/src/3rdparty/resonance-audio/resonance_audio/ambisonics/ambisonic_codec_impl.h b/src/3rdparty/resonance-audio/resonance_audio/ambisonics/ambisonic_codec_impl.h new file mode 100644 index 000000000..03a83140f --- /dev/null +++ b/src/3rdparty/resonance-audio/resonance_audio/ambisonics/ambisonic_codec_impl.h @@ -0,0 +1,295 @@ +/* +Copyright 2018 Google Inc. 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 RESONANCE_AUDIO_AMBISONICS_AMBISONIC_CODEC_IMPL_H_ +#define RESONANCE_AUDIO_AMBISONICS_AMBISONIC_CODEC_IMPL_H_ + +#include <cmath> +#include <vector> + +#include "Eigen/Dense" +#include "ambisonics/ambisonic_codec.h" +#include "ambisonics/associated_legendre_polynomials_generator.h" +#include "ambisonics/utils.h" +#include "base/audio_buffer.h" +#include "base/constants_and_types.h" +#include "base/logging.h" +#include "base/spherical_angle.h" +#include "utils/pseudoinverse.h" + +namespace vraudio { +// An encoder/decoder for ambisonic sound fields. It supports variable ambisonic +// order, ACN channel sequencing and SN3D normalization. +// +// @tparam NumAngles Used to fix the number of angles to be encoded/decoded at +// compile-time; use |Eigen::Dynamic| to indicate run-time variability. +// @tparam NumSphericalHarmonics Used to fix the number of spherical harmonic +// components at compile time; use |Eigen::Dynamic| to indicate run-time +// variability. +template <int NumAngles = Eigen::Dynamic, + int NumSphericalHarmonics = Eigen::Dynamic> +class AmbisonicCodecImpl : public AmbisonicCodec { + public: + // Spherical harmonics encoder matrix. + typedef Eigen::Matrix<float, NumSphericalHarmonics, NumAngles> EncoderMatrix; + // Spherical harmonics decoder matrix. + typedef Eigen::Matrix<float, NumAngles, NumSphericalHarmonics> DecoderMatrix; + // Spherical harmonics encoding of a frame or collection of frames (i.e., a + // vector of spherical harmonics). + typedef Eigen::Matrix<float, NumSphericalHarmonics, 1> EncodedVector; + // Decoded sequence of values for each angle / mono frame (i.e., a vector with + // a decoded mono frame for each angle). + typedef Eigen::Matrix<float, NumAngles, 1> DecodedVector; + + // Creates a codec with the given |ambisonic_order| and spherical |angles| to + // compute encoder/decoder matrices. + AmbisonicCodecImpl(int ambisonic_order, + const std::vector<SphericalAngle>& angles); + + // Implements |AmbisonicCodec|. + void EncodeBuffer(const AudioBuffer& input, AudioBuffer* output) override; + void DecodeBuffer(const AudioBuffer& input, AudioBuffer* output) override; + int ambisonic_order() const override; + size_t num_angles() const override; + size_t num_spherical_harmonics() const override; + const std::vector<SphericalAngle>& angles() const override; + void set_angles(const std::vector<SphericalAngle>& angles) override; + + // Encodes the given vector. + + void Encode(const Eigen::Ref<const Eigen::MatrixXf> decoded_vector, + Eigen::Ref<Eigen::MatrixXf> encoded_vector); + + // Decodes the given vector. + + void Decode(const Eigen::Ref<const Eigen::MatrixXf> encoded_vector, + Eigen::Ref<Eigen::MatrixXf> decoded_vector); + + // Gets the ambisonic sound field encoder matrix. + const Eigen::Ref<const Eigen::MatrixXf> GetEncoderMatrix(); + + // Gets the ambisonic sound field decoder matrix. + const Eigen::Ref<const Eigen::MatrixXf> GetDecoderMatrix(); + + // Necessary due to Eigen's alignment requirements on some platforms. + EIGEN_MAKE_ALIGNED_OPERATOR_NEW + + private: + // Returns the unnormalized spherical harmonic + // Y_degree^order(azimuth, elevation). + float UnnormalizedSphericalHarmonic(int degree, int order, + const SphericalAngle& angle) const; + + // The maximum-ordered ambisonic sound field handled by this codec. In the + // case of a periphonic codec, this is the order of the ambisonic sound field. + const int ambisonic_order_; + // Spherical angles used to compute spherical harmonics. For example, for a + // decoder, virtual loudspeaker positions; for an encoder, the position(s) of + // virtual sources relative to the listener. + std::vector<SphericalAngle> angles_; + // Current spherical harmonics encoder matrix if encoder_matrix_invalid_ is + // false. + EncoderMatrix encoder_matrix_; + // True if encoder_matrix_ needs to be recomputed. + bool encoder_matrix_invalid_; + // Current spherical harmonics decoder matrix if encoder_matrix_invalid_ is + // false. + DecoderMatrix decoder_matrix_; + // True if decoder_matrix_ needs to be recomputed. + bool decoder_matrix_invalid_; + // The associated Legendre polynomial generator for this codec. + AssociatedLegendrePolynomialsGenerator alp_generator_; + // Temporary storage for associated Legendre polynomials generated. + std::vector<float> associated_legendre_polynomials_temp_; +}; + +template <int NumAngles, int NumSphericalHarmonics> +AmbisonicCodecImpl<NumAngles, NumSphericalHarmonics>::AmbisonicCodecImpl( + int ambisonic_order, const std::vector<SphericalAngle>& angles) + : ambisonic_order_(ambisonic_order), + alp_generator_(ambisonic_order, false, false) { + DCHECK_GE(ambisonic_order_, 0); + set_angles(angles); +} + +template <int NumAngles, int NumSphericalHarmonics> +void AmbisonicCodecImpl<NumAngles, NumSphericalHarmonics>::Encode( + const Eigen::Ref<const Eigen::MatrixXf> decoded_vector, + Eigen::Ref<Eigen::MatrixXf> encoded_vector) { + encoded_vector.noalias() = GetEncoderMatrix() * decoded_vector; +} + +template <int NumAngles, int NumSphericalHarmonics> +void AmbisonicCodecImpl<NumAngles, NumSphericalHarmonics>::Decode( + const Eigen::Ref<const Eigen::MatrixXf> encoded_vector, + Eigen::Ref<Eigen::MatrixXf> decoded_vector) { + decoded_vector.noalias() = GetDecoderMatrix() * encoded_vector; +} + +template <int NumAngles, int NumSphericalHarmonics> +void AmbisonicCodecImpl<NumAngles, NumSphericalHarmonics>::EncodeBuffer( + const AudioBuffer& input, AudioBuffer* output) { + CHECK(output); + CHECK_EQ(input.num_channels(), num_angles()); + CHECK_EQ(output->num_channels(), num_spherical_harmonics()); + CHECK_EQ(input.num_frames(), output->num_frames()); + + Eigen::Map<const Eigen::Matrix<float, Eigen::Dynamic, Eigen::Dynamic, + Eigen::RowMajor>, + Eigen::Aligned> + unencoded_buffer(&input[0][0], num_angles(), output->GetChannelStride()); + + Eigen::Map< + Eigen::Matrix<float, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>, + Eigen::Aligned> + encoded_buffer(&(*output)[0][0], num_spherical_harmonics(), + input.GetChannelStride()); + + encoded_buffer.noalias() = GetEncoderMatrix() * unencoded_buffer; +} + +template <int NumAngles, int NumSphericalHarmonics> +void AmbisonicCodecImpl<NumAngles, NumSphericalHarmonics>::DecodeBuffer( + const AudioBuffer& input, AudioBuffer* output) { + CHECK(output); + CHECK_EQ(input.num_channels(), num_spherical_harmonics()); + CHECK_EQ(output->num_channels(), num_angles()); + CHECK_EQ(input.num_frames(), output->num_frames()); + + Eigen::Map<const Eigen::Matrix<float, Eigen::Dynamic, Eigen::Dynamic, + Eigen::RowMajor>, + Eigen::Aligned> + encoded_buffer(&input[0][0], num_spherical_harmonics(), + input.GetChannelStride()); + + Eigen::Map< + Eigen::Matrix<float, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>, + Eigen::Aligned> + decoded_buffer(&(*output)[0][0], num_angles(), + output->GetChannelStride()); + + decoded_buffer.noalias() = GetDecoderMatrix() * encoded_buffer; +} + +template <int NumAngles, int NumSphericalHarmonics> +const Eigen::Ref<const Eigen::MatrixXf> +AmbisonicCodecImpl<NumAngles, NumSphericalHarmonics>::GetEncoderMatrix() { + if (encoder_matrix_invalid_) { + encoder_matrix_ = EncoderMatrix( + GetNumPeriphonicComponents(ambisonic_order_), angles_.size()); + for (int col = 0; col < encoder_matrix_.cols(); col++) { + const SphericalAngle& angle = angles_[col]; + associated_legendre_polynomials_temp_ = + alp_generator_.Generate(std::sin(angle.elevation())); + // Compute the actual spherical harmonics using the generated polynomials. + for (int degree = 0; degree <= ambisonic_order_; degree++) { + for (int order = -degree; order <= degree; order++) { + const int row = AcnSequence(degree, order); + if (row == -1) { + // Skip this spherical harmonic. + continue; + } + encoder_matrix_(row, col) = + Sn3dNormalization(degree, order) * + UnnormalizedSphericalHarmonic(degree, order, angle); + } + } + } + encoder_matrix_invalid_ = false; + } + return encoder_matrix_; +} + +template <int NumAngles, int NumSphericalHarmonics> +const Eigen::Ref<const Eigen::MatrixXf> +AmbisonicCodecImpl<NumAngles, NumSphericalHarmonics>::GetDecoderMatrix() { + if (decoder_matrix_invalid_) { + decoder_matrix_ = Pseudoinverse<Eigen::MatrixXf>(GetEncoderMatrix()); + // Condition number of the encoding/decoding matrices. We use the fact that + // the decoding matrix is already a (pseudo)-inverse of the encoding matrix. + const float condition_number = + static_cast<float>(GetEncoderMatrix().norm() * decoder_matrix_.norm()); + const float num_rows = static_cast<float>(GetEncoderMatrix().rows()); + const float num_cols = static_cast<float>(GetEncoderMatrix().cols()); + if (condition_number > + 1.0f / (std::max(num_rows, num_cols) * kEpsilonFloat)) { + LOG(WARNING) << "Ambisonic decoding matrix is ill-conditioned. Results " + << "may be inaccurate."; + } + decoder_matrix_invalid_ = false; + } + return decoder_matrix_; +} + +template <int NumAngles, int NumSphericalHarmonics> +int AmbisonicCodecImpl<NumAngles, NumSphericalHarmonics>::ambisonic_order() + const { + return ambisonic_order_; +} + +template <int NumAngles, int NumSphericalHarmonics> +size_t AmbisonicCodecImpl<NumAngles, NumSphericalHarmonics>::num_angles() + const { + return angles_.size(); +} + +template <int NumAngles, int NumSphericalHarmonics> +size_t AmbisonicCodecImpl< + NumAngles, NumSphericalHarmonics>::num_spherical_harmonics() const { + // Return the worst-case scenario (the number of coefficients for a + // periphonic sound field). + return GetNumPeriphonicComponents(ambisonic_order_); +} + +template <int NumAngles, int NumSphericalHarmonics> +const std::vector<SphericalAngle>& +AmbisonicCodecImpl<NumAngles, NumSphericalHarmonics>::angles() const { + return angles_; +} + +template <int NumAngles, int NumSphericalHarmonics> +void AmbisonicCodecImpl<NumAngles, NumSphericalHarmonics>::set_angles( + const std::vector<SphericalAngle>& angles) { + CHECK_GT(angles.size(), 0); + + angles_ = angles; + encoder_matrix_invalid_ = decoder_matrix_invalid_ = true; +} + +template <int NumAngles, int NumSphericalHarmonics> +float AmbisonicCodecImpl<NumAngles, NumSphericalHarmonics>:: + UnnormalizedSphericalHarmonic(int degree, int order, + const SphericalAngle& angle) const { + const float last_term = + (order >= 0) ? std::cos(static_cast<float>(order) * angle.azimuth()) + : std::sin(static_cast<float>(-order) * angle.azimuth()); + return associated_legendre_polynomials_temp_[alp_generator_.GetIndex( + degree, std::abs(order))] * + last_term; +} + +// Codec for a single source. +template <int NumSphericalHarmonics = Eigen::Dynamic> +using MonoAmbisonicCodec = AmbisonicCodecImpl<1, NumSphericalHarmonics>; + +// Codec for a N-speaker first-order periphonic setup. +template <int NumAngles> +using FirstOrderPeriphonicAmbisonicCodec = + AmbisonicCodecImpl<NumAngles, GetNumPeriphonicComponentsStatic<1>::value>; + +} // namespace vraudio + +#endif // RESONANCE_AUDIO_AMBISONICS_AMBISONIC_CODEC_IMPL_H_ diff --git a/src/3rdparty/resonance-audio/resonance_audio/ambisonics/ambisonic_codec_impl_test.cc b/src/3rdparty/resonance-audio/resonance_audio/ambisonics/ambisonic_codec_impl_test.cc new file mode 100644 index 000000000..c7f0900be --- /dev/null +++ b/src/3rdparty/resonance-audio/resonance_audio/ambisonics/ambisonic_codec_impl_test.cc @@ -0,0 +1,354 @@ +/* +Copyright 2018 Google Inc. 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 "ambisonics/ambisonic_codec_impl.h" + +#include <cmath> +#include <memory> + +#include "third_party/googletest/googletest/include/gtest/gtest.h" +#include "base/audio_buffer.h" +#include "base/constants_and_types.h" +#include "base/misc_math.h" + +namespace vraudio { + +namespace { + +const int kCubeNumVertices = 8; +const float kCubeVertices[kCubeNumVertices][3] = { + {-1.0f, 1.0f, -1.0f}, {1.0f, 1.0f, -1.0f}, {-1.0f, 1.0f, 1.0f}, + {1.0f, 1.0f, 1.0f}, {-1.0f, -1.0f, -1.0f}, {1.0f, -1.0f, -1.0f}, + {-1.0f, -1.0f, 1.0f}, {1.0f, -1.0f, 1.0f}}; + +// Generates a third-order B-format encoder matrix. +AmbisonicCodecImpl<>::EncoderMatrix GenerateThirdOrderBFormatEncoderMatrix( + const SphericalAngle& spherical_angle) { + const size_t kNumSphericalHarmonics = GetNumPeriphonicComponents(3); + const size_t kNumAngles = 1; + AmbisonicCodecImpl<>::EncoderMatrix encoder_matrix(kNumSphericalHarmonics, + kNumAngles); + const float azimuth_rad = spherical_angle.azimuth(); + const float elevation_rad = spherical_angle.elevation(); + encoder_matrix << 1.0f, std::cos(elevation_rad) * std::sin(azimuth_rad), + std::sin(elevation_rad), std::cos(elevation_rad) * std::cos(azimuth_rad), + kSqrtThree / 2.0f * std::cos(elevation_rad) * std::cos(elevation_rad) * + std::sin(2.0f * azimuth_rad), + kSqrtThree / 2.0f * std::sin(2.0f * elevation_rad) * + std::sin(azimuth_rad), + 0.5f * (3.0f * std::sin(elevation_rad) * std::sin(elevation_rad) - 1.0f), + kSqrtThree / 2.0f * std::sin(2.0f * elevation_rad) * + std::cos(azimuth_rad), + kSqrtThree / 2.0f * std::cos(elevation_rad) * std::cos(elevation_rad) * + std::cos(2.0f * azimuth_rad), + std::sqrt(5.0f / 8.0f) * IntegerPow(std::cos(elevation_rad), 3) * + std::sin(3.0f * azimuth_rad), + std::sqrt(15.0f) / 2.0f * std::sin(elevation_rad) * + std::cos(elevation_rad) * std::cos(elevation_rad) * + std::sin(2.0f * azimuth_rad), + std::sqrt(3.0f / 8.0f) * std::cos(elevation_rad) * + (5.0f * std::sin(elevation_rad) * std::sin(elevation_rad) - 1.0f) * + std::sin(azimuth_rad), + 0.5f * std::sin(elevation_rad) * + (5.0f * std::sin(elevation_rad) * std::sin(elevation_rad) - 3.0f), + std::sqrt(3.0f / 8.0f) * std::cos(elevation_rad) * + (5.0f * std::sin(elevation_rad) * std::sin(elevation_rad) - 1.0f) * + std::cos(azimuth_rad), + std::sqrt(15.0f) / 2.0f * std::sin(elevation_rad) * + std::cos(elevation_rad) * std::cos(elevation_rad) * + std::cos(2.0f * azimuth_rad), + std::sqrt(5.0f / 8.0f) * IntegerPow(std::cos(elevation_rad), 3) * + std::cos(3.0f * azimuth_rad); + return encoder_matrix; +} + +class AmbisonicCodecTest : public ::testing::Test { + protected: + AmbisonicCodecTest() { + for (int i = 0; i < kCubeNumVertices; i++) { + const float azimuth_rad = + atan2f(kCubeVertices[i][2], kCubeVertices[i][0]); + const float elevation_rad = + std::asin(kCubeVertices[i][1] / + std::sqrt(kCubeVertices[i][0] * kCubeVertices[i][0] + + kCubeVertices[i][1] * kCubeVertices[i][1] + + kCubeVertices[i][2] * kCubeVertices[i][2])); + + cube_angles_.push_back(SphericalAngle(azimuth_rad, elevation_rad)); + } + + codec_first_order_cube_ = std::unique_ptr<AmbisonicCodecImpl<>>( + new AmbisonicCodecImpl<>(1, cube_angles_)); + } + + std::unique_ptr<AmbisonicCodecImpl<>> codec_first_order_cube_; + std::vector<SphericalAngle> cube_angles_; +}; + +// Tests that encoder and decoder matrices are ~inverses of each other. +TEST_F(AmbisonicCodecTest, EncoderDecoderIdentity) { + const auto encoder_matrix = codec_first_order_cube_->GetEncoderMatrix(); + const auto decoder_matrix = codec_first_order_cube_->GetDecoderMatrix(); + + const auto should_be_identity = encoder_matrix * decoder_matrix; + EXPECT_TRUE(should_be_identity.isApprox( + decltype(should_be_identity)::Identity(should_be_identity.rows(), + should_be_identity.cols()), + kEpsilonFloat)) + << "Matrix should be within " << kEpsilonFloat + << " of an identity matrix: \n" + << should_be_identity; +} + +// Tests that Encode and Decode are ~inverse operations. +TEST_F(AmbisonicCodecTest, EncodeDecodeIsInverse) { + AmbisonicCodecImpl<>::DecodedVector unencoded_vector(8); + unencoded_vector << 1, 2, 3, 4, 5, 6, 7, 8; + AmbisonicCodecImpl<>::EncodedVector encoded_vector(4); + codec_first_order_cube_->Encode(unencoded_vector, encoded_vector); + AmbisonicCodecImpl<>::DecodedVector decoded_vector(8); + codec_first_order_cube_->Decode(encoded_vector, decoded_vector); + + EXPECT_TRUE(unencoded_vector.isApprox(decoded_vector, kEpsilonFloat)) + << "Decoded vector should be within " << kEpsilonFloat << " of: \n" + << unencoded_vector; +} + +// Tests that EncodeBuffer and Decode (over a buffer) are ~inverse operations. +TEST_F(AmbisonicCodecTest, EncodeBufferDecodeVectorsIsInverse) { + const int kNumElements = 256; + const int kNumSphericalHarmonics = GetNumPeriphonicComponentsStatic<1>::value; + + AudioBuffer unencoded_buffer(kCubeNumVertices, kNumElements); + // Produce a buffer of successive numbers between -1 and 1. + for (int element = 0; element < kNumElements; ++element) { + for (int angle = 0; angle < kCubeNumVertices; ++angle) { + unencoded_buffer[angle][element] = + static_cast<float>(element * kCubeNumVertices + angle) / + (0.5f * kNumElements * kCubeNumVertices) - + 1.0f; + } + } + + AudioBuffer encoded_buffer(kNumSphericalHarmonics, kNumElements); + codec_first_order_cube_->EncodeBuffer(unencoded_buffer, &encoded_buffer); + + // Verify the encoded buffer and decoded vectors are ~identical. + for (int element = 0; element < kNumElements; ++element) { + AmbisonicCodecImpl<>::EncodedVector encoded_vector(kNumSphericalHarmonics); + for (int harmonic = 0; harmonic < kNumSphericalHarmonics; ++harmonic) { + encoded_vector[harmonic] = encoded_buffer[harmonic][element]; + } + AmbisonicCodecImpl<>::DecodedVector decoded_vector(kCubeNumVertices); + codec_first_order_cube_->Decode(encoded_vector, decoded_vector); + for (int angle = 0; angle < kCubeNumVertices; ++angle) { + EXPECT_NEAR(unencoded_buffer[angle][element], decoded_vector[angle], + kEpsilonFloat); + } + } +} + +// Tests that Encode (over a buffer) and DecodeBuffer are ~inverse operations. +TEST_F(AmbisonicCodecTest, EncodeVectorsDecodeBufferIsInverse) { + const int kNumElements = 256; + const int kNumSphericalHarmonics = GetNumPeriphonicComponentsStatic<1>::value; + + AudioBuffer unencoded_buffer(kCubeNumVertices, kNumElements); + // Produce a buffer of successive numbers between -1 and 1. + for (int element = 0; element < kNumElements; ++element) { + for (int angle = 0; angle < kCubeNumVertices; ++angle) { + unencoded_buffer[angle][element] = + static_cast<float>(element * kCubeNumVertices + angle) / + (0.5f * kNumElements * kCubeNumVertices) - + 1.0f; + } + } + + AudioBuffer encoded_buffer(kNumSphericalHarmonics, kNumElements); + // Produce the encoded version of unencoded_buffer. + for (int element = 0; element < kNumElements; ++element) { + AmbisonicCodecImpl<>::DecodedVector unencoded_vector(kCubeNumVertices); + for (int angle = 0; angle < kCubeNumVertices; ++angle) { + unencoded_vector[angle] = unencoded_buffer[angle][element]; + } + AmbisonicCodecImpl<>::EncodedVector encoded_vector(kNumSphericalHarmonics); + codec_first_order_cube_->Encode(unencoded_vector, encoded_vector); + for (int harmonic = 0; harmonic < kNumSphericalHarmonics; ++harmonic) { + encoded_buffer[harmonic][element] = encoded_vector[harmonic]; + } + } + + AudioBuffer decoded_buffer(kCubeNumVertices, kNumElements); + codec_first_order_cube_->DecodeBuffer(encoded_buffer, &decoded_buffer); + // Verify the decoded buffer and unencoded buffer are ~identical. + for (int element = 0; element < kNumElements; ++element) { + for (int angle = 0; angle < kCubeNumVertices; ++angle) { + EXPECT_NEAR(unencoded_buffer[angle][element], + decoded_buffer[angle][element], kEpsilonFloat); + } + } +} + +// Tests that EncodeBuffer and DecodeBuffer are ~inverse operations. +TEST_F(AmbisonicCodecTest, EncodeBufferDecodeBufferIsInverse) { + const int kNumElements = 256; + const int kNumSphericalHarmonics = GetNumPeriphonicComponentsStatic<1>::value; + + AudioBuffer unencoded_buffer(kCubeNumVertices, kNumElements); + // Produce a buffer of successive numbers between -1 and 1. + for (int element = 0; element < kNumElements; ++element) { + for (int angle = 0; angle < kCubeNumVertices; ++angle) { + unencoded_buffer[angle][element] = + static_cast<float>(element * kCubeNumVertices + angle) / + (0.5f * kNumElements * kCubeNumVertices) - + 1.0f; + } + } + + + + AudioBuffer encoded_buffer(kNumSphericalHarmonics, kNumElements); + // Produce the encoded version of unencoded_buffer. + for (int element = 0; element < kNumElements; ++element) { + AmbisonicCodecImpl<>::DecodedVector unencoded_vector(kCubeNumVertices); + for (int angle = 0; angle < kCubeNumVertices; ++angle) { + unencoded_vector[angle] = unencoded_buffer[angle][element]; + } + AmbisonicCodecImpl<>::EncodedVector encoded_vector(kNumSphericalHarmonics); + codec_first_order_cube_->Encode(unencoded_vector, encoded_vector); + for (int harmonic = 0; harmonic < kNumSphericalHarmonics; ++harmonic) { + encoded_buffer[harmonic][element] = encoded_vector[harmonic]; + } + } + + AudioBuffer decoded_buffer(kCubeNumVertices, kNumElements); + codec_first_order_cube_->DecodeBuffer(encoded_buffer, &decoded_buffer); + // Verify the decoded buffer and unencoded buffer are ~identical. + for (int element = 0; element < kNumElements; ++element) { + for (int angle = 0; angle < kCubeNumVertices; ++angle) { + EXPECT_NEAR(unencoded_buffer[angle][element], + decoded_buffer[angle][element], kEpsilonFloat); + } + } +} + +// Tests that third-order encoding coefficients produced by Codec are correct. +TEST_F(AmbisonicCodecTest, ThirdOrderEncodingCoefficients) { + const int kNumSphericalHarmonics = GetNumPeriphonicComponentsStatic<3>::value; + + const float kAzimuthStart = -180.0f; + const float kAzimuthStop = 180.0f; + const float kAzimuthStep = 10.0f; + for (float azimuth = kAzimuthStart; azimuth <= kAzimuthStop; + azimuth += kAzimuthStep) { + const float kElevationStart = -90.0f; + const float kElevationStop = 90.0f; + const float kElevationStep = 10.0f; + for (float elevation = kElevationStart; elevation <= kElevationStop; + elevation += kElevationStep) { + const int kAmbisonicOrder = 3; + + const SphericalAngle kSphericalAngle = + SphericalAngle::FromDegrees(azimuth, elevation); + AmbisonicCodecImpl<> mono_codec(kAmbisonicOrder, {kSphericalAngle}); + + AmbisonicCodecImpl<>::DecodedVector kUnencodedVector(1); + kUnencodedVector << 0.12345f; + AmbisonicCodecImpl<>::EncodedVector codec_encoded_vector( + kNumSphericalHarmonics); + mono_codec.Encode(kUnencodedVector, codec_encoded_vector); + const AmbisonicCodecImpl<>::EncoderMatrix + expected_b_format_encoder_matrix = + GenerateThirdOrderBFormatEncoderMatrix(kSphericalAngle); + const AmbisonicCodecImpl<>::EncodedVector expected_encoded_vector = + static_cast<AmbisonicCodecImpl<>::EncodedVector>( + expected_b_format_encoder_matrix * kUnencodedVector); + ASSERT_TRUE( + codec_encoded_vector.isApprox(expected_encoded_vector, kEpsilonFloat)) + << "Codec encoded vector: \n" + << codec_encoded_vector << "\n should be within " << kEpsilonFloat + << " of expected: \n" + << expected_encoded_vector; + } + } +} + +// Tests that Sphere16 decoding matrix (ACN/SN3D) coefficients produced by +// Codec +// are correct. +TEST_F(AmbisonicCodecTest, Sphere16DecoderCoefficients) { + // Expected decoder coefficients for the Sphere16 setup (ACN/SN3D) obtained + // using //matlab/ambix/ambdecodematrix.m. + const std::vector<float> kExpectedDecoderCoefficients{ + 0.0625000000000001f, 0.0625000000000000f, 0.0625000000000000f, + 0.0625000000000000f, 0.0625000000000000f, 0.0625000000000000f, + 0.0625000000000000f, 0.0625000000000000f, 0.0625000000000000f, + 0.0625000000000000f, 0.0625000000000000f, 0.0625000000000000f, + 0.0625000000000000f, 0.0625000000000000f, 0.0625000000000000f, + 0.0625000000000000f, -1.94257951433956e-17f, 0.128971506904330f, + 0.182393254223798f, 0.128971506904330f, 1.73262494006158e-16f, + -0.128971506904330f, -0.182393254223798f, -0.128971506904330f, + -2.79581233842125e-16f, 0.116817928199481f, -1.55203460683255e-16f, + -0.116817928199482f, 0.0826027491940168f, 0.0826027491940169f, + -0.0826027491940164f, -0.0826027491940165f, -0.0440164745323702f, + 0.0440164745323702f, -0.0440164745323703f, 0.0440164745323702f, + -0.0440164745323703f, 0.0440164745323702f, -0.0440164745323702f, + 0.0440164745323703f, 0.146497161016369f, 0.146497161016369f, + 0.146497161016369f, 0.146497161016369f, -0.146497161016369f, + -0.146497161016369f, -0.146497161016369f, -0.146497161016369f, + 0.182393254223799f, 0.128971506904330f, 7.26758979552257e-17f, + -0.128971506904330f, -0.182393254223798f, -0.128971506904330f, + -5.02371557522155e-17f, 0.128971506904330f, 0.116817928199482f, + 6.83999538850465e-17f, -0.116817928199482f, -3.47728677633385e-17f, + 0.0826027491940167f, -0.0826027491940166f, -0.0826027491940166f, + 0.0826027491940166f}; + const int kAmbisonicOrder = 1; + const int kNumChannels = 4; + const int kNumVirtualSpeakers = 16; + std::vector<SphericalAngle> sphere16_angles{ + SphericalAngle::FromDegrees(0.0f, -13.6f), + SphericalAngle::FromDegrees(45.0f, 13.6f), + SphericalAngle::FromDegrees(90.0f, -13.6f), + SphericalAngle::FromDegrees(135.0f, 13.6f), + SphericalAngle::FromDegrees(180.0f, -13.6f), + SphericalAngle::FromDegrees(-135.0f, 13.6f), + SphericalAngle::FromDegrees(-90.0f, -13.6f), + SphericalAngle::FromDegrees(-45.0f, 13.6f), + SphericalAngle::FromDegrees(0.0f, 51.5f), + SphericalAngle::FromDegrees(90.0f, 51.5f), + SphericalAngle::FromDegrees(180.0f, 51.5f), + SphericalAngle::FromDegrees(-90.0f, 51.5f), + SphericalAngle::FromDegrees(45.0f, -51.5f), + SphericalAngle::FromDegrees(135.0f, -51.5f), + SphericalAngle::FromDegrees(-135.0f, -51.5f), + SphericalAngle::FromDegrees(-45.0f, -51.5f)}; + + std::unique_ptr<AmbisonicCodecImpl<>> codec_first_order_sphere16( + new AmbisonicCodecImpl<>(kAmbisonicOrder, sphere16_angles)); + const AmbisonicCodecImpl<>::DecoderMatrix decoder_matrix = + codec_first_order_sphere16->GetDecoderMatrix(); + // Check if the size of the decoding matrix is correct. + ASSERT_EQ(decoder_matrix.size(), kNumVirtualSpeakers * kNumChannels); + // Check each coefficient against MATLAB. + for (size_t i = 0; i < kExpectedDecoderCoefficients.size(); ++i) { + EXPECT_NEAR(kExpectedDecoderCoefficients[i], decoder_matrix(i), + kEpsilonFloat); + } +} + +} // namespace + +} // namespace vraudio diff --git a/src/3rdparty/resonance-audio/resonance_audio/ambisonics/ambisonic_lookup_table.cc b/src/3rdparty/resonance-audio/resonance_audio/ambisonics/ambisonic_lookup_table.cc new file mode 100644 index 000000000..2a2bc1b68 --- /dev/null +++ b/src/3rdparty/resonance-audio/resonance_audio/ambisonics/ambisonic_lookup_table.cc @@ -0,0 +1,231 @@ +/* +Copyright 2018 Google Inc. 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 "ambisonics/ambisonic_lookup_table.h" + +#include "ambisonics/ambisonic_spread_coefficients.h" +#include "ambisonics/associated_legendre_polynomials_generator.h" +#include "ambisonics/utils.h" +#include "base/constants_and_types.h" +#include "base/misc_math.h" + + +namespace vraudio { + +namespace { + +// Number of azimuth angles to store in the pre-computed encoder lookup table +// (for 0 - 90 degrees using 1 degree increments). +const size_t kNumAzimuths = 91; + +// Number of elevation angles to store in the pre-computed encoder lookup table +// (for 0 - 90 degrees using 1 degree increments). +const size_t kNumElevations = 91; + +// Number of Cartesian axes for which to compute spherical harmonic symmetries. +const size_t kNumAxes = 3; + +// Minimum angular source spreads at different orders for which we start to +// apply spread gain correction coefficients. Array index corresponds to +// ambisonic order. For more information about sound source spread control, +// please refer to the Matlab code and the corresponding paper. +const int kMinSpreads[kMaxSupportedAmbisonicOrder + 1] = {361, 54, 40, 31}; + +// Spread coefficients are stored sequentially in a 1-d table. Therefore, to +// access coefficients for the required ambisonic order we need to apply offsets +// which are equal to the total number of coefficients for previous orders. For +// example, the total number of coefficient in each ambisonic order 'n' is +// equal to the number of unique orders multiplied by number of unique spreads, +// i.e.: (n + 1) * (360 - kMinSpreads[n] + 1). +const int kSpreadCoeffOffsets[kMaxSupportedAmbisonicOrder + 1] = {0, 1, 615, + 1578}; + +size_t GetEncoderTableIndex(size_t i, size_t j, size_t k, size_t width, + size_t depth) { + return i * width * depth + j * depth + k; +} + +size_t GetSymmetriesTableIndex(size_t i, size_t j, size_t width) { + return i * width + j; +} + +int GetSpreadTableIndex(int ambisonic_order, float source_spread_deg) { + // Number of unique ambisonic orders in the sound field. For example, in a + // first order sound field we have components of order 0 and 1. + const int num_unique_orders = ambisonic_order + 1; + return kSpreadCoeffOffsets[ambisonic_order] + + num_unique_orders * (static_cast<int>(source_spread_deg) - + kMinSpreads[ambisonic_order]); +} + +} // namespace + +AmbisonicLookupTable::AmbisonicLookupTable(int max_ambisonic_order) + : max_ambisonic_order_(max_ambisonic_order), + max_num_coeffs_in_table_( + GetNumPeriphonicComponents(max_ambisonic_order_) - 1), + encoder_table_(kNumAzimuths * kNumElevations * max_num_coeffs_in_table_), + symmetries_table_(kNumAxes * max_num_coeffs_in_table_) { + DCHECK_GE(max_ambisonic_order_, 0); + DCHECK_LE(max_ambisonic_order_, kMaxSupportedAmbisonicOrder); + ComputeEncoderTable(); + ComputeSymmetriesTable(); +} + +void AmbisonicLookupTable::GetEncodingCoeffs( + int ambisonic_order, const SphericalAngle& source_direction, + float source_spread_deg, std::vector<float>* encoding_coeffs) const { + + DCHECK_GE(ambisonic_order, 0); + DCHECK_LE(ambisonic_order, kMaxSupportedAmbisonicOrder); + DCHECK_GE(source_spread_deg, 0.0f); + DCHECK_LE(source_spread_deg, 360.0f); + DCHECK(encoding_coeffs); + // Raw coefficients are stored only for ambisonic orders 1 and up, since 0th + // order raw coefficient is always 1. + const size_t num_raw_coeffs = GetNumPeriphonicComponents(ambisonic_order) - 1; + // The actual number of returned Ambisonic coefficients is therefore + // |num_raw_coeffs + 1|. + DCHECK_EQ(encoding_coeffs->size(), num_raw_coeffs + 1); + DCHECK_GE(source_direction.azimuth(), -kPi); + DCHECK_LE(source_direction.azimuth(), kTwoPi); + DCHECK_GE(source_direction.elevation(), -kHalfPi); + DCHECK_LE(source_direction.elevation(), kHalfPi); + const int azimuth_deg = + source_direction.azimuth() < kPi + ? static_cast<int>(source_direction.azimuth() * kDegreesFromRadians) + : static_cast<int>(source_direction.azimuth() * kDegreesFromRadians) - + 360; + const int elevation_deg = + static_cast<int>(source_direction.elevation() * kDegreesFromRadians); + const size_t abs_azimuth_deg = static_cast<size_t>(std::abs(azimuth_deg)); + const size_t azimuth_idx = + abs_azimuth_deg > 90 ? 180 - abs_azimuth_deg : abs_azimuth_deg; + const size_t elevation_idx = static_cast<size_t>(std::abs(elevation_deg)); + (*encoding_coeffs)[0] = 1.0f; + for (size_t raw_coeff_idx = 0; raw_coeff_idx < num_raw_coeffs; + ++raw_coeff_idx) { + // Variable to hold information about spherical harmonic phase flip. 1 means + // no flip; -1 means 180 degrees flip. + float flip = 1.0f; + // The following three 'if' statements implement the logic to correct the + // phase of the current spherical harmonic, depending on which quadrant the + // sound source is located in. For more information, please see the Matlab + // code and the corresponding paper. + if (azimuth_deg < 0) { + flip = symmetries_table_[GetSymmetriesTableIndex( + 0, raw_coeff_idx, max_num_coeffs_in_table_)]; + } + if (elevation_deg < 0) { + flip *= symmetries_table_[GetSymmetriesTableIndex( + 1, raw_coeff_idx, max_num_coeffs_in_table_)]; + } + if (abs_azimuth_deg > 90) { + flip *= symmetries_table_[GetSymmetriesTableIndex( + 2, raw_coeff_idx, max_num_coeffs_in_table_)]; + } + const size_t encoder_table_idx = + GetEncoderTableIndex(azimuth_idx, elevation_idx, raw_coeff_idx, + kNumElevations, max_num_coeffs_in_table_); + (*encoding_coeffs)[raw_coeff_idx + 1] = + encoder_table_[encoder_table_idx] * flip; + } + + // If the spread is more than min. theoretical spread for the given + // |ambisonic_order|, multiply the encoding coefficients by the required + // spread control gains from the |kSpreadCoeffs| lookup table. + if (source_spread_deg >= kMinSpreads[ambisonic_order]) { + const int spread_table_idx = + GetSpreadTableIndex(ambisonic_order, source_spread_deg); + (*encoding_coeffs)[0] *= kSpreadCoeffs[spread_table_idx]; + for (size_t coeff = 1; coeff < encoding_coeffs->size(); ++coeff) { + const int current_coefficient_degree = + GetPeriphonicAmbisonicOrderForChannel(coeff); + (*encoding_coeffs)[coeff] *= + kSpreadCoeffs[spread_table_idx + current_coefficient_degree]; + } + } +} + +void AmbisonicLookupTable::ComputeEncoderTable() { + + // Associated Legendre polynomial generator. + AssociatedLegendrePolynomialsGenerator alp_generator( + max_ambisonic_order_, /*condon_shortley_phase=*/false, + /*compute_negative_order=*/false); + // Temporary storage for associated Legendre polynomials generated. + std::vector<float> temp_associated_legendre_polynomials; + for (size_t azimuth_idx = 0; azimuth_idx < kNumAzimuths; ++azimuth_idx) { + for (size_t elevation_idx = 0; elevation_idx < kNumElevations; + ++elevation_idx) { + const SphericalAngle angle( + static_cast<float>(azimuth_idx) * kRadiansFromDegrees, + static_cast<float>(elevation_idx) * kRadiansFromDegrees); + temp_associated_legendre_polynomials = + alp_generator.Generate(std::sin(angle.elevation())); + // First spherical harmonic is always equal 1 for all angles so we do not + // need to compute and store it. + for (int degree = 1; degree <= max_ambisonic_order_; ++degree) { + for (int order = -degree; order <= degree; ++order) { + const size_t alp_index = + alp_generator.GetIndex(degree, std::abs(order)); + const float alp_value = + temp_associated_legendre_polynomials[alp_index]; + const size_t raw_coeff_idx = AcnSequence(degree, order) - 1; + const size_t encoder_table_idx = + GetEncoderTableIndex(azimuth_idx, elevation_idx, raw_coeff_idx, + kNumElevations, max_num_coeffs_in_table_); + encoder_table_[encoder_table_idx] = + Sn3dNormalization(degree, order) * + UnnormalizedSphericalHarmonic(alp_value, order, angle.azimuth()); + } + } + } + } +} + +void AmbisonicLookupTable::ComputeSymmetriesTable() { + + for (int degree = 1; degree <= max_ambisonic_order_; ++degree) { + for (int order = -degree; order <= degree; ++order) { + const size_t raw_coeff_idx = AcnSequence(degree, order) - 1; + // Symmetry wrt the left-right axis (Y). + symmetries_table_[GetSymmetriesTableIndex(0, raw_coeff_idx, + max_num_coeffs_in_table_)] = + order < 0 ? -1.0f : 1.0f; + // Symmetry wrt the up-down axis (Z). + symmetries_table_[GetSymmetriesTableIndex(1, raw_coeff_idx, + max_num_coeffs_in_table_)] = + static_cast<float>(IntegerPow(-1, degree + order)); + // Symmetry wrt the front-back axis (X). + symmetries_table_[GetSymmetriesTableIndex(2, raw_coeff_idx, + max_num_coeffs_in_table_)] = + order < 0 ? static_cast<float>(-IntegerPow(-1, std::abs(order))) + : static_cast<float>(IntegerPow(-1, order)); + } + } +} + +float AmbisonicLookupTable::UnnormalizedSphericalHarmonic(float alp_value, + int order, + float azimuth_rad) { + const float horizontal_term = + (order >= 0) ? std::cos(static_cast<float>(order) * azimuth_rad) + : std::sin(static_cast<float>(-order) * azimuth_rad); + return alp_value * horizontal_term; +} + +} // namespace vraudio diff --git a/src/3rdparty/resonance-audio/resonance_audio/ambisonics/ambisonic_lookup_table.h b/src/3rdparty/resonance-audio/resonance_audio/ambisonics/ambisonic_lookup_table.h new file mode 100644 index 000000000..096068d0a --- /dev/null +++ b/src/3rdparty/resonance-audio/resonance_audio/ambisonics/ambisonic_lookup_table.h @@ -0,0 +1,92 @@ +/* +Copyright 2018 Google Inc. 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 RESONANCE_AUDIO_AMBISONICS_AMBISONIC_LOOKUP_TABLE_H_ +#define RESONANCE_AUDIO_AMBISONICS_AMBISONIC_LOOKUP_TABLE_H_ + +#include <vector> + +#include "base/spherical_angle.h" + +namespace vraudio { + +// Represents a lookup table for encoding of Ambisonic periphonic sound fields. +// Supports arbitrary Ambisonic order and uses AmbiX convention (ACN channel +// sequencing, SN3D normalization). +class AmbisonicLookupTable { + public: + // Creates Ambisonic (AmbiX) encoder lookup table given the + // |max_ambisonic_order| used by the client. AmbiX convention uses ACN channel + // sequencing and SN3D normalization. + explicit AmbisonicLookupTable(int max_ambisonic_order); + + // Gets spherical harmonic encoding coefficients for a given order and + // writes them to |encoding_coeffs|. + // + // @param ambisonic_order Ambisonic order of the encoded sound source. + // @param source_direction Direction of a sound source in spherical + // coordinates. + // @param source_spread_deg Encoded sound source angular spread in degrees. + // @param encoding_coeffs Pointer to a vector of Ambisonic encoding + // coefficients. + void GetEncodingCoeffs(int ambisonic_order, + const SphericalAngle& source_direction, + float source_spread_deg, + std::vector<float>* encoding_coeffs) const; + + private: + // Computes a lookup table of encoding coefficients for one quadrant of the + // sphere. + void ComputeEncoderTable(); + + // Computes a table of spherical harmonics symmetry information for all + // cartesian axes. Value of 1 indicates that the current spherical harmonic is + // symmetric with respect to the current axis. Value of -1 indicates that the + // current spherical harmonic is anti-symmetric with respect to the current + // axis. + void ComputeSymmetriesTable(); + + // Returns the unnormalized spherical harmonic coefficient: + // Y_degree^order(azimuth, elevation). + // + // @param alp_value Associated Legendre polynomial for the given degree and + // order evaluated at sin elevation angle: P_degree^order(sin(elevation)). + // @param order Order of the Associated Legendre polynomial. + // @param azimuth_rad Azimuth angle in radians. + // @return Unnormalized spherical harmonic coefficient. + float UnnormalizedSphericalHarmonic(float alp_value, int order, + float azimuth_rad); + + // Ambisonic order. + const int max_ambisonic_order_; + + // Maximum number of coefficients to be stored in the lookup table equal to + // the number of Ambisonic channels for |max_ambisonic_order_| - 1. This is + // because we do not need to store the coefficient for the first spherical + // harmonic coefficient as it is always 1. + const size_t max_num_coeffs_in_table_; + + // Lookup table for storing Ambisonic encoding coefficients for one quadrant + // of the sphere. + std::vector<float> encoder_table_; + + // Lookup table for storing information about spherical harmonic symmetries. + std::vector<float> symmetries_table_; +}; + +} // namespace vraudio + +#endif // RESONANCE_AUDIO_AMBISONICS_AMBISONIC_LOOKUP_TABLE_H_ diff --git a/src/3rdparty/resonance-audio/resonance_audio/ambisonics/ambisonic_lookup_table_test.cc b/src/3rdparty/resonance-audio/resonance_audio/ambisonics/ambisonic_lookup_table_test.cc new file mode 100644 index 000000000..8f0c9c9a6 --- /dev/null +++ b/src/3rdparty/resonance-audio/resonance_audio/ambisonics/ambisonic_lookup_table_test.cc @@ -0,0 +1,240 @@ +/* +Copyright 2018 Google Inc. 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 "ambisonics/ambisonic_lookup_table.h" + +#include <cmath> + +#include "third_party/googletest/googletest/include/gtest/gtest.h" +#include "ambisonics/utils.h" +#include "base/constants_and_types.h" +#include "base/misc_math.h" + +namespace vraudio { + +namespace { + +const int kSourceAmbisonicOrder0 = 0; +const int kSourceAmbisonicOrder1 = 1; +const int kSourceAmbisonicOrder2 = 2; +const int kSourceAmbisonicOrder3 = 3; + +// Minimum angular source spread of 0 ensures that no gain correction +// coefficients should be applied to the Ambisonic encoding coefficients. +const float kMinSpreadDeg = 0.0f; + +} // namespace + +class AmbisonicLookupTableTest : public ::testing::TestWithParam<int> { + protected: + AmbisonicLookupTableTest() + : lookup_table_(kMaxSupportedAmbisonicOrder), + source_ambisonic_order_(GetParam()), + encoding_coeffs_(GetNumPeriphonicComponents(source_ambisonic_order_)) {} + + // Tests whether GetEncodingCoeffs() method returns correct coefficients. + void TestEncodingCoefficients() { + for (size_t i = 0; i < kSourceDirections.size(); ++i) { + const std::vector<float> kExpectedCoeffs = + GenerateExpectedCoeffs(kSourceDirections[i]); + lookup_table_.GetEncodingCoeffs(source_ambisonic_order_, + kSourceDirections[i], kMinSpreadDeg, + &encoding_coeffs_); + for (size_t j = 0; j < encoding_coeffs_.size(); ++j) { + EXPECT_NEAR(kExpectedCoeffs[j], encoding_coeffs_[j], kEpsilonFloat); + } + } + } + + // Tests whether changing the source spread results in correct gain changes in + // pressure and velocity ambisonic encoding coefficients. + void TestSpreadEnergy() { + // Choose an artbitrary source direction from |kSourceDirections|. + const SphericalAngle source_direction = kSourceDirections.back(); + lookup_table_.GetEncodingCoeffs(source_ambisonic_order_, source_direction, + kMinSpreadDeg, &encoding_coeffs_); + float current_pressure_energy = GetPressureEnergy(); + float current_velocity_energy = GetVelocityEnergy(); + for (int spread_deg = 1; spread_deg <= 360; ++spread_deg) { + lookup_table_.GetEncodingCoeffs(source_ambisonic_order_, source_direction, + static_cast<float>(spread_deg), + &encoding_coeffs_); + float new_pressure_energy = GetPressureEnergy(); + float new_velocity_energy = GetVelocityEnergy(); + EXPECT_TRUE(new_pressure_energy >= current_pressure_energy); + EXPECT_TRUE(new_velocity_energy <= current_velocity_energy); + current_pressure_energy = new_pressure_energy; + current_velocity_energy = new_velocity_energy; + } + } + + private: + // Generates expected ambisonic encoding coefficients for ambisonic orders 0 + // to 3, according to http://ambisonics.ch/standards/channels/index. + std::vector<float> GenerateExpectedCoeffs(const SphericalAngle& angle) { + const float a = angle.azimuth(); + const float e = angle.elevation(); + return std::vector<float>{ + 1.0f, + std::cos(e) * std::sin(a), + std::sin(e), + std::cos(e) * std::cos(a), + kSqrtThree / 2.0f * std::cos(e) * std::cos(e) * std::sin(2.0f * a), + kSqrtThree / 2.0f * std::sin(2.0f * e) * std::sin(a), + 0.5f * (3.0f * std::sin(e) * std::sin(e) - 1.0f), + kSqrtThree / 2.0f * std::sin(2.0f * e) * std::cos(a), + kSqrtThree / 2.0f * std::cos(e) * std::cos(e) * std::cos(2.0f * a), + std::sqrt(5.0f / 8.0f) * IntegerPow(std::cos(e), 3) * + std::sin(3.0f * a), + std::sqrt(15.0f) / 2.0f * std::sin(e) * std::cos(e) * std::cos(e) * + std::sin(2.0f * a), + std::sqrt(3.0f / 8.0f) * std::cos(e) * + (5.0f * std::sin(e) * std::sin(e) - 1.0f) * std::sin(a), + 0.5f * std::sin(e) * (5.0f * std::sin(e) * std::sin(e) - 3.0f), + std::sqrt(3.0f / 8.0f) * std::cos(e) * + (5.0f * std::sin(e) * std::sin(e) - 1.0f) * std::cos(a), + std::sqrt(15.0f) / 2.0f * std::sin(e) * std::cos(e) * std::cos(e) * + std::cos(2.0f * a), + std::sqrt(5.0f / 8.0f) * IntegerPow(std::cos(e), 3) * + std::cos(3.0f * a)}; + } + + // Computes energy of the pressure channel (Ambisonic channel 0). + float GetPressureEnergy() { + return encoding_coeffs_[0] * encoding_coeffs_[0]; + } + + // Computes total energy of all the velocity channels (Ambisonic channel 1 and + // above). + float GetVelocityEnergy() { + float velocity_energy = 0.0f; + for (size_t i = 1; i < encoding_coeffs_.size(); ++i) { + velocity_energy += encoding_coeffs_[i] * encoding_coeffs_[i]; + } + return velocity_energy; + } + + // Source angles corresponding to each Cartesian axis direction as well as + // some arbitrary directions (full degrees) in each quadrant of the sphere. + const std::vector<SphericalAngle> kSourceDirections = { + {SphericalAngle::FromDegrees(0.0f, 0.0f)} /* Front */, + {SphericalAngle::FromDegrees(90.0f, 0.0f)} /* Left */, + {SphericalAngle::FromDegrees(180.0f, 0.0f)} /* Back */, + {SphericalAngle::FromDegrees(-90.0f, 0.0f)} /* Right */, + {SphericalAngle::FromDegrees(0.0f, 90.0f)} /* Up */, + {SphericalAngle::FromDegrees(0.0f, -90.0f)} /* Down */, + {SphericalAngle::FromDegrees(10.0f, 20.0f)} /* Left-Top-Front */, + {SphericalAngle::FromDegrees(-30.0f, 40.0f)} /* Right-Top-Front */, + {SphericalAngle::FromDegrees(50.0f, -60.0f)} /* Left-Bottom-Front */, + {SphericalAngle::FromDegrees(290.0f, -80.0f)} /* Right-Bottom-Front */, + {SphericalAngle::FromDegrees(110.0f, 5.0f)} /* Left-Top-Back */, + {SphericalAngle::FromDegrees(-120.0f, 15.0f)} /* Right-Top-Back */, + {SphericalAngle::FromDegrees(130.0f, -25.0f)} /* Left-Bottom-Back */, + {SphericalAngle::FromDegrees(220.0f, -35.0f)} /* Right-Bottom-Back */, + }; + + AmbisonicLookupTable lookup_table_; + const int source_ambisonic_order_; + std::vector<float> encoding_coeffs_; +}; + +// Tests whether GetEncodingCoeffs() method returns correct coefficients. +TEST_P(AmbisonicLookupTableTest, GetEncodingCoeffsTest) { + TestEncodingCoefficients(); +} + +// Tests whether changing the source spread results in correct gain changes in +// pressure and velocity ambisonic encoding coefficients. For example, +// increasing of the source spread should result in overall monotonic increase +// of energy in the pressure channel and overall monotonic decrease in the +// velocity channels. +TEST_P(AmbisonicLookupTableTest, SpreadEnergyTest) { TestSpreadEnergy(); } + +INSTANTIATE_TEST_CASE_P(TestParameters, AmbisonicLookupTableTest, + testing::Values(kSourceAmbisonicOrder0, + kSourceAmbisonicOrder1, + kSourceAmbisonicOrder2, + kSourceAmbisonicOrder3)); + +class PreComputedCoeffsTest : public ::testing::Test { + protected: + // Tests whether the lookup table returns correct coefficients for sources + // with arbitrary ambisonic order, direction and spread. + void TestSpreadCoefficients() { + // Test spread coefficients for each config. + for (const auto& config : kConfigs) { + const int source_ambisonic_order = config.source_ambisonic_order; + const SphericalAngle& source_direction = config.source_direction; + const float source_spread_deg = config.source_spread_deg; + const std::vector<float>& expected_coeffs = config.expected_coefficients; + std::vector<float> encoding_coeffs( + GetNumPeriphonicComponents(source_ambisonic_order)); + AmbisonicLookupTable lookup_table(kMaxSupportedAmbisonicOrder); + lookup_table.GetEncodingCoeffs(source_ambisonic_order, source_direction, + source_spread_deg, &encoding_coeffs); + for (size_t i = 0; i < encoding_coeffs.size(); ++i) { + EXPECT_NEAR(expected_coeffs[i], encoding_coeffs[i], kEpsilonFloat); + } + } + } + + private: + struct TestConfig { + int source_ambisonic_order; + SphericalAngle source_direction; + float source_spread_deg; + std::vector<float> expected_coefficients; + }; + + const std::vector<TestConfig> kConfigs = { + // Zeroth order sound source. + {0 /* ambisonic order */, + SphericalAngle::FromDegrees(0.0f, 0.0f) /* source direction */, + 120.0f /* source spread */, + {1.0f} /* expected coefficients */}, + + // First order sound source. + {1 /* ambisonic order */, + SphericalAngle::FromDegrees(36.0f, 45.0f) /* source direction */, + 70.0f /* source spread */, + {1.20046f, 0.310569f, 0.528372f, 0.427462f} + /* expected coefficients */}, + + // Second order sound source. + {2 /* ambisonic order */, + SphericalAngle::FromDegrees(55.0f, -66.0f) /* source direction */, + 41.0f /* source spread */, + {1.038650f, 0.337027f, -0.924096f, 0.2359899f, 0.124062f, -0.485807f, + 0.6928289f, -0.340166f, -0.045155f} + /* expected coefficients */}, + + // Third order sound source. + {3 /* ambisonic order */, + SphericalAngle::FromDegrees(-13.0f, 90.0f) /* source direction */, + 32.0f /* source spread */, + {1.03237f, 0.0f, 1.02119f, 0.0f, 0.0f, 0.0f, 0.990433f, 0.0f, 0.0f, 0.0f, + 0.0f, 0.0f, 0.898572f, 0.0f, 0.0f, 0.0f} + /* expected coefficients */}}; +}; + +// Tests whether the lookup table returns correct coefficients for sources with +// arbitrary ambisonic order, direction and spread. The expected coefficients +// have been pre-computed using Matlab. +TEST_F(PreComputedCoeffsTest, SpreadCoefficientsTest) { + TestSpreadCoefficients(); +} + +} // namespace vraudio diff --git a/src/3rdparty/resonance-audio/resonance_audio/ambisonics/ambisonic_spread_coefficients.h b/src/3rdparty/resonance-audio/resonance_audio/ambisonics/ambisonic_spread_coefficients.h new file mode 100644 index 000000000..0c245b87f --- /dev/null +++ b/src/3rdparty/resonance-audio/resonance_audio/ambisonics/ambisonic_spread_coefficients.h @@ -0,0 +1,610 @@ +/* +Copyright 2018 Google Inc. 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 RESONANCE_AUDIO_AMBISONICS_AMBISONIC_SPREAD_COEFFICIENTS_H_ +#define RESONANCE_AUDIO_AMBISONICS_AMBISONIC_SPREAD_COEFFICIENTS_H_ + +namespace vraudio { + +// Lookup table with gain coefficients to control Ambisonic sound source spread +// at different orders, up to and including order 3. For information on how the +// spread control coefficients have been derived please refer to Matlab code as +// well as the corresponding paper. +static const float kSpreadCoeffs[5089] = { + 1.0f, 1.0019f, 0.998092f, 1.01441f, 0.985357f, + 1.02704f, 0.972176f, 1.03975f, 0.958555f, 1.05252f, + 0.944499f, 1.06532f, 0.930017f, 1.07813f, 0.915119f, + 1.09091f, 0.899819f, 1.10365f, 0.884131f, 1.11631f, + 0.868073f, 1.12886f, 0.851665f, 1.14128f, 0.834929f, + 1.15354f, 0.817888f, 1.16561f, 0.800569f, 1.17747f, + 0.782999f, 1.1891f, 0.76521f, 1.20046f, 0.747231f, + 1.21155f, 0.729095f, 1.22234f, 0.710837f, 1.23281f, + 0.692489f, 1.24296f, 0.674086f, 1.25276f, 0.655664f, + 1.26221f, 0.637256f, 1.27131f, 0.618896f, 1.28003f, + 0.600619f, 1.28839f, 0.582456f, 1.29637f, 0.564439f, + 1.30399f, 0.546597f, 1.31123f, 0.528959f, 1.31811f, + 0.511551f, 1.32463f, 0.494397f, 1.3308f, 0.47752f, + 1.33663f, 0.46094f, 1.34212f, 0.444677f, 1.34729f, + 0.428746f, 1.35214f, 0.413161f, 1.3567f, 0.397935f, + 1.36096f, 0.383079f, 1.36495f, 0.3686f, 1.36867f, + 0.354505f, 1.37215f, 0.340799f, 1.37538f, 0.327486f, + 1.37839f, 0.314566f, 1.38118f, 0.302039f, 1.38378f, + 0.289906f, 1.38618f, 0.278163f, 1.38841f, 0.266808f, + 1.39047f, 0.255835f, 1.39238f, 0.24524f, 1.39414f, + 0.235017f, 1.39576f, 0.225159f, 1.39726f, 0.21566f, + 1.39864f, 0.206511f, 1.39991f, 0.197704f, 1.40108f, + 0.189232f, 1.40215f, 0.181084f, 1.40314f, 0.173254f, + 1.40405f, 0.165731f, 1.40488f, 0.158507f, 1.40565f, + 0.151572f, 1.40635f, 0.144918f, 1.40699f, 0.138535f, + 1.40758f, 0.132414f, 1.40812f, 0.126546f, 1.40861f, + 0.120923f, 1.40906f, 0.115536f, 1.40947f, 0.110376f, + 1.40985f, 0.105436f, 1.4102f, 0.100705f, 1.41051f, + 0.0961781f, 1.4108f, 0.0918459f, 1.41107f, 0.087701f, + 1.41131f, 0.0837362f, 1.41153f, 0.0799443f, 1.41173f, + 0.0763185f, 1.41191f, 0.0728519f, 1.41208f, 0.069538f, + 1.41223f, 0.0663707f, 1.41237f, 0.0633438f, 1.41249f, + 0.0604515f, 1.41261f, 0.0576881f, 1.41271f, 0.0550481f, + 1.41281f, 0.0525264f, 1.4129f, 0.0501178f, 1.41298f, + 0.0478176f, 1.41305f, 0.045621f, 1.41312f, 0.0435236f, + 1.41318f, 0.0415211f, 1.41323f, 0.0396092f, 1.41328f, + 0.0377841f, 1.41333f, 0.0360419f, 1.41337f, 0.034379f, + 1.4134f, 0.0327918f, 1.41344f, 0.031277f, 1.41347f, + 0.0298315f, 1.4135f, 0.028452f, 1.41352f, 0.0271356f, + 1.41355f, 0.0258796f, 1.41357f, 0.0246811f, 1.41359f, + 0.0235377f, 1.41361f, 0.0224468f, 1.41362f, 0.0214061f, + 1.41364f, 0.0204132f, 1.41365f, 0.0194661f, 1.41366f, + 0.0185626f, 1.41367f, 0.0177008f, 1.41368f, 0.0168787f, + 1.41369f, 0.0160946f, 1.4137f, 0.0153468f, 1.41371f, + 0.0146335f, 1.41372f, 0.0139531f, 1.41372f, 0.0133043f, + 1.41373f, 0.0126855f, 1.41373f, 0.0120953f, 1.41374f, + 0.0115325f, 1.41374f, 0.0109958f, 1.41375f, 0.0104839f, + 1.41375f, 0.00999584f, 1.41375f, 0.00953039f, 1.41376f, + 0.00908654f, 1.41376f, 0.0086633f, 1.41376f, 0.00825972f, + 1.41376f, 0.00787488f, 1.41376f, 0.00750793f, 1.41377f, + 0.00715803f, 1.41377f, 0.0068244f, 1.41377f, 0.00650628f, + 1.41377f, 0.00620296f, 1.41377f, 0.00591376f, 1.41377f, + 0.00563801f, 1.41377f, 0.00537509f, 1.41378f, 0.00512441f, + 1.41378f, 0.0048854f, 1.41378f, 0.00465752f, 1.41378f, + 0.00444026f, 1.41378f, 0.00423311f, 1.41378f, 0.00403562f, + 1.41378f, 0.00384732f, 1.41378f, 0.0036678f, 1.41378f, + 0.00349665f, 1.41378f, 0.00333348f, 1.41378f, 0.00317791f, + 1.41378f, 0.00302959f, 1.41378f, 0.00288819f, 1.41378f, + 0.00275338f, 1.41378f, 0.00262486f, 1.41378f, 0.00250233f, + 1.41378f, 0.00238552f, 1.41378f, 0.00227415f, 1.41378f, + 0.00216798f, 1.41378f, 0.00206677f, 1.41378f, 0.00197027f, + 1.41378f, 0.00187828f, 1.41378f, 0.00179058f, 1.41378f, + 0.00170697f, 1.41378f, 0.00162727f, 1.41378f, 0.00155128f, + 1.41378f, 0.00147884f, 1.41378f, 0.00140979f, 1.41378f, + 0.00134395f, 1.41378f, 0.00128119f, 1.41378f, 0.00122136f, + 1.41378f, 0.00116432f, 1.41378f, 0.00110994f, 1.41378f, + 0.0010581f, 1.41378f, 0.00100869f, 1.41378f, 0.000961575f, + 1.41378f, 0.000916665f, 1.41378f, 0.000873851f, 1.41378f, + 0.000833036f, 1.41378f, 0.000794127f, 1.41378f, 0.000757035f, + 1.41378f, 0.000721675f, 1.41378f, 0.000687966f, 1.41378f, + 0.000655831f, 1.41378f, 0.000625197f, 1.41378f, 0.000595994f, + 1.41378f, 0.000568154f, 1.41378f, 0.000541615f, 1.41378f, + 0.000516315f, 1.41378f, 0.000492197f, 1.41378f, 0.000469205f, + 1.41378f, 0.000447287f, 1.41378f, 0.000426393f, 1.41378f, + 0.000406474f, 1.41378f, 0.000387486f, 1.41378f, 0.000369385f, + 1.41378f, 0.000352129f, 1.41378f, 0.00033568f, 1.41378f, + 0.000319998f, 1.41378f, 0.00030505f, 1.41378f, 0.000290799f, + 1.41378f, 0.000277214f, 1.41378f, 0.000264264f, 1.41378f, + 0.000251918f, 1.41378f, 0.00024015f, 1.41378f, 0.000228931f, + 1.41378f, 0.000218236f, 1.41378f, 0.000208041f, 1.41378f, + 0.000198322f, 1.41378f, 0.000189056f, 1.41378f, 0.000180224f, + 1.41378f, 0.000171805f, 1.41378f, 0.000163778f, 1.41378f, + 0.000156127f, 1.41378f, 0.000148833f, 1.41378f, 0.00014188f, + 1.41378f, 0.000135252f, 1.41378f, 0.000128933f, 1.41378f, + 0.000122909f, 1.41378f, 0.000117167f, 1.41378f, 0.000111693f, + 1.41378f, 0.000106475f, 1.41378f, 0.000101501f, 1.41378f, + 9.67588e-05f, 1.41378f, 9.22384e-05f, 1.41378f, 8.79291e-05f, + 1.41378f, 8.38211e-05f, 1.41378f, 7.99051e-05f, 1.41378f, + 7.6172e-05f, 1.41378f, 7.26133e-05f, 1.41378f, 6.92209e-05f, + 1.41378f, 6.5987e-05f, 1.41378f, 6.29041e-05f, 1.41378f, + 5.99653e-05f, 1.41378f, 5.71637e-05f, 1.41378f, 5.44931e-05f, + 1.41378f, 5.19472e-05f, 1.41378f, 4.95202e-05f, 1.41378f, + 4.72067e-05f, 1.41378f, 4.50012e-05f, 1.41378f, 4.28988e-05f, + 1.41378f, 4.08946e-05f, 1.41378f, 3.8984e-05f, 1.41378f, + 3.71627e-05f, 1.41378f, 3.54264e-05f, 1.41378f, 3.37713e-05f, + 1.41378f, 3.21935e-05f, 1.41378f, 3.06895e-05f, 1.41378f, + 2.92557e-05f, 1.41378f, 2.78888e-05f, 1.41378f, 2.65859e-05f, + 1.41378f, 2.53438e-05f, 1.41378f, 2.41597e-05f, 1.41378f, + 2.3031e-05f, 1.41378f, 2.1955e-05f, 1.41378f, 2.09292e-05f, + 1.41378f, 1.99514e-05f, 1.41378f, 1.90193e-05f, 1.41378f, + 1.81307e-05f, 1.41378f, 1.72837e-05f, 1.41378f, 1.64762e-05f, + 1.41378f, 1.57064e-05f, 1.41378f, 1.49726e-05f, 1.41378f, + 1.42731e-05f, 1.41378f, 1.36062e-05f, 1.41378f, 1.29706e-05f, + 1.41378f, 1.23646e-05f, 1.41378f, 1.17869e-05f, 1.41378f, + 1.12362e-05f, 1.41378f, 1.07113e-05f, 1.41378f, 1.02108e-05f, + 1.41378f, 9.73377e-06f, 1.41378f, 9.27901e-06f, 1.41378f, + 8.84549e-06f, 1.41378f, 8.43223e-06f, 1.41378f, 8.03827e-06f, + 1.41378f, 7.66273e-06f, 1.41378f, 7.30472e-06f, 1.41378f, + 6.96344e-06f, 1.41378f, 6.63811e-06f, 1.41378f, 6.32798e-06f, + 1.41378f, 6.03233e-06f, 1.41378f, 5.7505e-06f, 1.41378f, + 5.48184e-06f, 1.41378f, 5.22572e-06f, 1.41378f, 4.98158e-06f, + 1.41378f, 4.74884e-06f, 1.41378f, 4.52697e-06f, 1.41378f, + 4.31547e-06f, 1.41378f, 4.11385e-06f, 1.41378f, 3.92165e-06f, + 1.41378f, 3.73843e-06f, 1.41378f, 3.56377e-06f, 1.41378f, + 3.39727e-06f, 1.41378f, 3.23855e-06f, 1.41378f, 3.08724e-06f, + 1.41378f, 2.94301e-06f, 1.41378f, 2.80551e-06f, 1.41378f, + 2.67443e-06f, 1.41378f, 2.54948e-06f, 1.41378f, 2.43037e-06f, + 1.41378f, 2.31682e-06f, 1.41378f, 2.20858e-06f, 1.41378f, + 2.1054e-06f, 1.41378f, 2.00703e-06f, 1.41378f, 1.91326e-06f, + 1.41378f, 1.82388e-06f, 1.41378f, 1.73866e-06f, 1.41378f, + 1.65743e-06f, 1.41378f, 1.58e-06f, 1.41378f, 1.50618e-06f, + 1.41378f, 1.43581e-06f, 1.41378f, 1.36873e-06f, 1.41378f, + 1.30478e-06f, 1.41378f, 1.24382e-06f, 1.41378f, 1.18571e-06f, + 1.01483f, 1.00465f, 0.970383f, 1.03865f, 1.01155f, + 0.9215f, 1.06314f, 1.0179f, 0.869372f, 1.08819f, + 1.02356f, 0.813901f, 1.11364f, 1.02839f, 0.755004f, + 1.13935f, 1.03225f, 0.692625f, 1.16513f, 1.03497f, + 0.626735f, 1.19076f, 1.0364f, 0.557346f, 1.21602f, + 1.03638f, 0.484517f, 1.24065f, 1.03474f, 0.408362f, + 1.2644f, 1.03133f, 0.329059f, 1.28697f, 1.02601f, + 0.246853f, 1.3081f, 1.01866f, 0.162061f, 1.3275f, + 1.00918f, 0.0750718f, 1.34495f, 0.99753f, 0.0f, + 1.36175f, 0.984802f, 0.0f, 1.37869f, 0.971629f, + 0.0f, 1.39575f, 0.958015f, 0.0f, 1.4129f, + 0.943967f, 0.0f, 1.43008f, 0.929493f, 0.0f, + 1.44728f, 0.914603f, 0.0f, 1.46444f, 0.899312f, + 0.0f, 1.48154f, 0.883633f, 0.0f, 1.49853f, + 0.867584f, 0.0f, 1.51538f, 0.851185f, 0.0f, + 1.53206f, 0.834458f, 0.0f, 1.54851f, 0.817427f, + 0.0f, 1.56472f, 0.800117f, 0.0f, 1.58064f, + 0.782558f, 0.0f, 1.59624f, 0.764779f, 0.0f, + 1.6115f, 0.74681f, 0.0f, 1.62638f, 0.728685f, + 0.0f, 1.64086f, 0.710436f, 0.0f, 1.65492f, + 0.692098f, 0.0f, 1.66854f, 0.673706f, 0.0f, + 1.6817f, 0.655294f, 0.0f, 1.69439f, 0.636896f, + 0.0f, 1.7066f, 0.618548f, 0.0f, 1.71831f, + 0.600281f, 0.0f, 1.72953f, 0.582128f, 0.0f, + 1.74025f, 0.564121f, 0.0f, 1.75047f, 0.546289f, + 0.0f, 1.76019f, 0.528661f, 0.0f, 1.76943f, + 0.511262f, 0.0f, 1.77819f, 0.494118f, 0.0f, + 1.78647f, 0.477251f, 0.0f, 1.79429f, 0.46068f, + 0.0f, 1.80166f, 0.444426f, 0.0f, 1.8086f, + 0.428504f, 0.0f, 1.81511f, 0.412928f, 0.0f, + 1.82123f, 0.397711f, 0.0f, 1.82695f, 0.382863f, + 0.0f, 1.8323f, 0.368392f, 0.0f, 1.8373f, + 0.354306f, 0.0f, 1.84196f, 0.340607f, 0.0f, + 1.84631f, 0.327301f, 0.0f, 1.85035f, 0.314388f, + 0.0f, 1.8541f, 0.301869f, 0.0f, 1.85758f, + 0.289743f, 0.0f, 1.86081f, 0.278006f, 0.0f, + 1.8638f, 0.266657f, 0.0f, 1.86657f, 0.255691f, + 0.0f, 1.86912f, 0.245102f, 0.0f, 1.87149f, + 0.234885f, 0.0f, 1.87367f, 0.225033f, 0.0f, + 1.87568f, 0.215538f, 0.0f, 1.87753f, 0.206394f, + 0.0f, 1.87924f, 0.197593f, 0.0f, 1.8808f, + 0.189125f, 0.0f, 1.88225f, 0.180982f, 0.0f, + 1.88357f, 0.173156f, 0.0f, 1.88479f, 0.165638f, + 0.0f, 1.88591f, 0.158418f, 0.0f, 1.88694f, + 0.151487f, 0.0f, 1.88788f, 0.144836f, 0.0f, + 1.88874f, 0.138457f, 0.0f, 1.88953f, 0.132339f, + 0.0f, 1.89025f, 0.126475f, 0.0f, 1.89091f, + 0.120855f, 0.0f, 1.89152f, 0.115471f, 0.0f, + 1.89207f, 0.110314f, 0.0f, 1.89258f, 0.105376f, + 0.0f, 1.89305f, 0.100649f, 0.0f, 1.89347f, + 0.0961239f, 0.0f, 1.89386f, 0.0917941f, 0.0f, + 1.89421f, 0.0876516f, 0.0f, 1.89453f, 0.083689f, + 0.0f, 1.89483f, 0.0798993f, 0.0f, 1.8951f, + 0.0762755f, 0.0f, 1.89534f, 0.0728108f, 0.0f, + 1.89557f, 0.0694989f, 0.0f, 1.89577f, 0.0663333f, + 0.0f, 1.89596f, 0.0633081f, 0.0f, 1.89613f, + 0.0604174f, 0.0f, 1.89628f, 0.0576555f, 0.0f, + 1.89642f, 0.0550171f, 0.0f, 1.89655f, 0.0524968f, + 0.0f, 1.89667f, 0.0500896f, 0.0f, 1.89678f, + 0.0477907f, 0.0f, 1.89687f, 0.0455953f, 0.0f, + 1.89696f, 0.0434991f, 0.0f, 1.89704f, 0.0414977f, + 0.0f, 1.89712f, 0.0395869f, 0.0f, 1.89718f, + 0.0377628f, 0.0f, 1.89724f, 0.0360216f, 0.0f, + 1.8973f, 0.0343596f, 0.0f, 1.89735f, 0.0327733f, + 0.0f, 1.8974f, 0.0312594f, 0.0f, 1.89744f, + 0.0298146f, 0.0f, 1.89748f, 0.0284359f, 0.0f, + 1.89751f, 0.0271203f, 0.0f, 1.89754f, 0.025865f, + 0.0f, 1.89757f, 0.0246672f, 0.0f, 1.8976f, + 0.0235244f, 0.0f, 1.89762f, 0.0224342f, 0.0f, + 1.89764f, 0.021394f, 0.0f, 1.89766f, 0.0204017f, + 0.0f, 1.89768f, 0.0194551f, 0.0f, 1.8977f, + 0.0185521f, 0.0f, 1.89771f, 0.0176908f, 0.0f, + 1.89773f, 0.0168692f, 0.0f, 1.89774f, 0.0160856f, + 0.0f, 1.89775f, 0.0153381f, 0.0f, 1.89776f, + 0.0146252f, 0.0f, 1.89777f, 0.0139453f, 0.0f, + 1.89778f, 0.0132968f, 0.0f, 1.89778f, 0.0126783f, + 0.0f, 1.89779f, 0.0120885f, 0.0f, 1.8978f, + 0.011526f, 0.0f, 1.8978f, 0.0109896f, 0.0f, + 1.89781f, 0.010478f, 0.0f, 1.89781f, 0.00999021f, + 0.0f, 1.89782f, 0.00952502f, 0.0f, 1.89782f, + 0.00908142f, 0.0f, 1.89782f, 0.00865842f, 0.0f, + 1.89783f, 0.00825506f, 0.0f, 1.89783f, 0.00787044f, + 0.0f, 1.89783f, 0.00750369f, 0.0f, 1.89784f, + 0.00715399f, 0.0f, 1.89784f, 0.00682055f, 0.0f, + 1.89784f, 0.00650262f, 0.0f, 1.89784f, 0.00619947f, + 0.0f, 1.89784f, 0.00591042f, 0.0f, 1.89785f, + 0.00563483f, 0.0f, 1.89785f, 0.00537206f, 0.0f, + 1.89785f, 0.00512152f, 0.0f, 1.89785f, 0.00488265f, + 0.0f, 1.89785f, 0.0046549f, 0.0f, 1.89785f, + 0.00443776f, 0.0f, 1.89785f, 0.00423073f, 0.0f, + 1.89785f, 0.00403334f, 0.0f, 1.89785f, 0.00384516f, + 0.0f, 1.89785f, 0.00366574f, 0.0f, 1.89785f, + 0.00349468f, 0.0f, 1.89786f, 0.0033316f, 0.0f, + 1.89786f, 0.00317612f, 0.0f, 1.89786f, 0.00302788f, + 0.0f, 1.89786f, 0.00288656f, 0.0f, 1.89786f, + 0.00275183f, 0.0f, 1.89786f, 0.00262338f, 0.0f, + 1.89786f, 0.00250092f, 0.0f, 1.89786f, 0.00238417f, + 0.0f, 1.89786f, 0.00227287f, 0.0f, 1.89786f, + 0.00216676f, 0.0f, 1.89786f, 0.0020656f, 0.0f, + 1.89786f, 0.00196916f, 0.0f, 1.89786f, 0.00187722f, + 0.0f, 1.89786f, 0.00178957f, 0.0f, 1.89786f, + 0.00170601f, 0.0f, 1.89786f, 0.00162635f, 0.0f, + 1.89786f, 0.00155041f, 0.0f, 1.89786f, 0.00147801f, + 0.0f, 1.89786f, 0.00140899f, 0.0f, 1.89786f, + 0.00134319f, 0.0f, 1.89786f, 0.00128047f, 0.0f, + 1.89786f, 0.00122067f, 0.0f, 1.89786f, 0.00116366f, + 0.0f, 1.89786f, 0.00110932f, 0.0f, 1.89786f, + 0.00105751f, 0.0f, 1.89786f, 0.00100812f, 0.0f, + 1.89786f, 0.000961033f, 0.0f, 1.89786f, 0.000916148f, + 0.0f, 1.89786f, 0.000873358f, 0.0f, 1.89786f, + 0.000832566f, 0.0f, 1.89786f, 0.000793679f, 0.0f, + 1.89786f, 0.000756608f, 0.0f, 1.89786f, 0.000721268f, + 0.0f, 1.89786f, 0.000687578f, 0.0f, 1.89786f, + 0.000655462f, 0.0f, 1.89786f, 0.000624845f, 0.0f, + 1.89786f, 0.000595658f, 0.0f, 1.89786f, 0.000567834f, + 0.0f, 1.89786f, 0.00054131f, 0.0f, 1.89786f, + 0.000516024f, 0.0f, 1.89786f, 0.000491919f, 0.0f, + 1.89786f, 0.000468941f, 0.0f, 1.89786f, 0.000447035f, + 0.0f, 1.89786f, 0.000426152f, 0.0f, 1.89786f, + 0.000406245f, 0.0f, 1.89786f, 0.000387268f, 0.0f, + 1.89786f, 0.000369177f, 0.0f, 1.89786f, 0.000351931f, + 0.0f, 1.89786f, 0.000335491f, 0.0f, 1.89786f, + 0.000319818f, 0.0f, 1.89786f, 0.000304878f, 0.0f, + 1.89786f, 0.000290635f, 0.0f, 1.89786f, 0.000277058f, + 0.0f, 1.89786f, 0.000264115f, 0.0f, 1.89786f, + 0.000251776f, 0.0f, 1.89786f, 0.000240014f, 0.0f, + 1.89786f, 0.000228802f, 0.0f, 1.89786f, 0.000218113f, + 0.0f, 1.89786f, 0.000207923f, 0.0f, 1.89786f, + 0.00019821f, 0.0f, 1.89786f, 0.00018895f, 0.0f, + 1.89786f, 0.000180123f, 0.0f, 1.89786f, 0.000171708f, + 0.0f, 1.89786f, 0.000163686f, 0.0f, 1.89786f, + 0.000156039f, 0.0f, 1.89786f, 0.000148749f, 0.0f, + 1.89786f, 0.0001418f, 0.0f, 1.89786f, 0.000135175f, + 0.0f, 1.89786f, 0.00012886f, 0.0f, 1.89786f, + 0.00012284f, 0.0f, 1.89786f, 0.000117101f, 0.0f, + 1.89786f, 0.00011163f, 0.0f, 1.89786f, 0.000106415f, + 0.0f, 1.89786f, 0.000101444f, 0.0f, 1.89786f, + 9.67043e-05f, 0.0f, 1.89786f, 9.21864e-05f, 0.0f, + 1.89786f, 8.78795e-05f, 0.0f, 1.89786f, 8.37739e-05f, + 0.0f, 1.89786f, 7.98601e-05f, 0.0f, 1.89786f, + 7.61291e-05f, 0.0f, 1.89786f, 7.25724e-05f, 0.0f, + 1.89786f, 6.91819e-05f, 0.0f, 1.89786f, 6.59498e-05f, + 0.0f, 1.89786f, 6.28686e-05f, 0.0f, 1.89786f, + 5.99315e-05f, 0.0f, 1.89786f, 5.71315e-05f, 0.0f, + 1.89786f, 5.44624e-05f, 0.0f, 1.89786f, 5.19179e-05f, + 0.0f, 1.89786f, 4.94923e-05f, 0.0f, 1.89786f, + 4.71801e-05f, 0.0f, 1.89786f, 4.49758e-05f, 0.0f, + 1.89786f, 4.28746e-05f, 0.0f, 1.89786f, 4.08715e-05f, + 0.0f, 1.89786f, 3.8962e-05f, 0.0f, 1.89786f, + 3.71417e-05f, 0.0f, 1.89786f, 3.54065e-05f, 0.0f, + 1.89786f, 3.37523e-05f, 0.0f, 1.89786f, 3.21754e-05f, + 0.0f, 1.89786f, 3.06722e-05f, 0.0f, 1.89786f, + 2.92392e-05f, 0.0f, 1.89786f, 2.78731e-05f, 0.0f, + 1.89786f, 2.65709e-05f, 0.0f, 1.89786f, 2.53295e-05f, + 0.0f, 1.89786f, 2.41461e-05f, 0.0f, 1.89786f, + 2.3018e-05f, 0.0f, 1.89786f, 2.19426e-05f, 0.0f, + 1.89786f, 2.09175e-05f, 0.0f, 1.89786f, 1.99402e-05f, + 0.0f, 1.89786f, 1.90086e-05f, 0.0f, 1.89786f, + 1.81205e-05f, 0.0f, 1.89786f, 1.72739e-05f, 0.0f, + 1.89786f, 1.64669e-05f, 0.0f, 1.89786f, 1.56975e-05f, + 0.0f, 1.89786f, 1.49642e-05f, 0.0f, 1.89786f, + 1.4265e-05f, 0.0f, 1.89786f, 1.35986e-05f, 0.0f, + 1.89786f, 1.29632e-05f, 0.0f, 1.89786f, 1.23576e-05f, + 0.0f, 1.89786f, 1.17802e-05f, 0.0f, 1.89786f, + 1.12299e-05f, 0.0f, 1.89786f, 1.07052e-05f, 0.0f, + 1.89786f, 1.02051e-05f, 0.0f, 1.89786f, 9.72828e-06f, + 0.0f, 1.89786f, 9.27378e-06f, 0.0f, 1.89786f, + 8.84051e-06f, 0.0f, 1.89786f, 8.42748e-06f, 0.0f, + 1.89786f, 8.03374e-06f, 0.0f, 1.89786f, 7.65841e-06f, + 0.0f, 1.89786f, 7.3006e-06f, 0.0f, 1.89786f, + 6.95952e-06f, 0.0f, 1.89786f, 6.63437e-06f, 0.0f, + 1.89786f, 6.32441e-06f, 0.0f, 1.89786f, 6.02893e-06f, + 0.0f, 1.89786f, 5.74726e-06f, 0.0f, 1.89786f, + 5.47875e-06f, 0.0f, 1.89786f, 5.22278e-06f, 0.0f, + 1.89786f, 4.97877e-06f, 0.0f, 1.89786f, 4.74616e-06f, + 0.0f, 1.89786f, 4.52442e-06f, 0.0f, 1.89786f, + 4.31304e-06f, 0.0f, 1.89786f, 4.11153e-06f, 0.0f, + 1.89786f, 3.91944e-06f, 0.0f, 1.89786f, 3.73632e-06f, + 0.0f, 1.89786f, 3.56176e-06f, 0.0f, 1.89786f, + 3.39536e-06f, 0.0f, 1.89786f, 3.23672e-06f, 0.0f, + 1.89786f, 3.0855e-06f, 0.0f, 1.89786f, 2.94135e-06f, + 0.0f, 1.89786f, 2.80393e-06f, 0.0f, 1.89786f, + 2.67293e-06f, 0.0f, 1.89786f, 2.54805e-06f, 0.0f, + 1.89786f, 2.429e-06f, 0.0f, 1.89786f, 2.31552e-06f, + 0.0f, 1.89786f, 2.20734e-06f, 0.0f, 1.89786f, + 2.10421e-06f, 0.0f, 1.89786f, 2.0059e-06f, 0.0f, + 1.89786f, 1.91218e-06f, 0.0f, 1.89786f, 1.82285e-06f, + 0.0f, 1.89786f, 1.73768e-06f, 0.0f, 1.89786f, + 1.6565e-06f, 0.0f, 1.89786f, 1.57911e-06f, 0.0f, + 1.89786f, 1.50533e-06f, 0.0f, 1.89786f, 1.435e-06f, + 0.0f, 1.89786f, 1.36796e-06f, 0.0f, 1.89786f, + 1.30405e-06f, 0.0f, 1.89786f, 1.24312e-06f, 0.0f, + 1.89786f, 1.18504e-06f, 0.0f, 1.00324f, 1.00216f, + 0.999152f, 0.990038f, 1.03237f, 1.02119f, 0.990433f, + 0.898572f, 1.06269f, 1.04023f, 0.979161f, 0.799806f, + 1.094f, 1.05895f, 0.964976f, 0.693603f, 1.126f, + 1.07701f, 0.947526f, 0.57989f, 1.15835f, 1.09398f, + 0.926474f, 0.45869f, 1.19059f, 1.10944f, 0.901512f, + 0.330158f, 1.22223f, 1.12289f, 0.87237f, 0.194621f, + 1.25268f, 1.13384f, 0.838839f, 0.0526136f, 1.28199f, + 1.14236f, 0.801199f, 0.0f, 1.31207f, 1.15021f, + 0.760839f, 0.0f, 1.34301f, 1.15742f, 0.717799f, + 0.0f, 1.37465f, 1.16386f, 0.671999f, 0.0f, + 1.40681f, 1.16935f, 0.623371f, 0.0f, 1.43929f, + 1.17374f, 0.571868f, 0.0f, 1.47185f, 1.17684f, + 0.517465f, 0.0f, 1.50423f, 1.17846f, 0.460174f, + 0.0f, 1.53613f, 1.17844f, 0.400043f, 0.0f, + 1.56725f, 1.17657f, 0.337165f, 0.0f, 1.59725f, + 1.1727f, 0.271688f, 0.0f, 1.62577f, 1.16665f, + 0.203815f, 0.0f, 1.65245f, 1.15828f, 0.133806f, + 0.0f, 1.67697f, 1.14751f, 0.0619832f, 0.0f, + 1.69901f, 1.13426f, 0.0f, 0.0f, 1.72022f, + 1.11979f, 0.0f, 0.0f, 1.74163f, 1.10481f, + 0.0f, 0.0f, 1.76318f, 1.08933f, 0.0f, + 0.0f, 1.78484f, 1.07336f, 0.0f, 0.0f, + 1.80655f, 1.0569f, 0.0f, 0.0f, 1.82827f, + 1.03997f, 0.0f, 0.0f, 1.84995f, 1.02258f, + 0.0f, 0.0f, 1.87155f, 1.00475f, 0.0f, + 0.0f, 1.89302f, 0.986504f, 0.0f, 0.0f, + 1.91431f, 0.967857f, 0.0f, 0.0f, 1.93537f, + 0.948837f, 0.0f, 0.0f, 1.95615f, 0.929471f, + 0.0f, 0.0f, 1.97662f, 0.90979f, 0.0f, + 0.0f, 1.99674f, 0.889823f, 0.0f, 0.0f, + 2.01645f, 0.869607f, 0.0f, 0.0f, 2.03572f, + 0.849175f, 0.0f, 0.0f, 2.05452f, 0.828565f, + 0.0f, 0.0f, 2.07282f, 0.807816f, 0.0f, + 0.0f, 2.09058f, 0.786964f, 0.0f, 0.0f, + 2.10779f, 0.766051f, 0.0f, 0.0f, 2.12441f, + 0.745115f, 0.0f, 0.0f, 2.14044f, 0.724196f, + 0.0f, 0.0f, 2.15586f, 0.703332f, 0.0f, + 0.0f, 2.17065f, 0.682561f, 0.0f, 0.0f, + 2.18482f, 0.661921f, 0.0f, 0.0f, 2.19836f, + 0.641445f, 0.0f, 0.0f, 2.21128f, 0.621169f, + 0.0f, 0.0f, 2.22356f, 0.601125f, 0.0f, + 0.0f, 2.23523f, 0.581341f, 0.0f, 0.0f, + 2.24629f, 0.561847f, 0.0f, 0.0f, 2.25675f, + 0.542667f, 0.0f, 0.0f, 2.26663f, 0.523826f, + 0.0f, 0.0f, 2.27594f, 0.505344f, 0.0f, + 0.0f, 2.28471f, 0.487239f, 0.0f, 0.0f, + 2.29294f, 0.469528f, 0.0f, 0.0f, 2.30066f, + 0.452225f, 0.0f, 0.0f, 2.30789f, 0.435342f, + 0.0f, 0.0f, 2.31465f, 0.418888f, 0.0f, + 0.0f, 2.32097f, 0.40287f, 0.0f, 0.0f, + 2.32686f, 0.387294f, 0.0f, 0.0f, 2.33234f, + 0.372164f, 0.0f, 0.0f, 2.33745f, 0.357481f, + 0.0f, 0.0f, 2.34219f, 0.343246f, 0.0f, + 0.0f, 2.34659f, 0.329458f, 0.0f, 0.0f, + 2.35066f, 0.316113f, 0.0f, 0.0f, 2.35444f, + 0.303208f, 0.0f, 0.0f, 2.35794f, 0.290738f, + 0.0f, 0.0f, 2.36117f, 0.278698f, 0.0f, + 0.0f, 2.36415f, 0.26708f, 0.0f, 0.0f, + 2.36691f, 0.255878f, 0.0f, 0.0f, 2.36945f, + 0.245082f, 0.0f, 0.0f, 2.37179f, 0.234685f, + 0.0f, 0.0f, 2.37394f, 0.224677f, 0.0f, + 0.0f, 2.37592f, 0.215048f, 0.0f, 0.0f, + 2.37775f, 0.20579f, 0.0f, 0.0f, 2.37942f, + 0.196891f, 0.0f, 0.0f, 2.38096f, 0.188342f, + 0.0f, 0.0f, 2.38237f, 0.180132f, 0.0f, + 0.0f, 2.38367f, 0.172251f, 0.0f, 0.0f, + 2.38486f, 0.164689f, 0.0f, 0.0f, 2.38595f, + 0.157435f, 0.0f, 0.0f, 2.38694f, 0.150479f, + 0.0f, 0.0f, 2.38786f, 0.143811f, 0.0f, + 0.0f, 2.38869f, 0.137421f, 0.0f, 0.0f, + 2.38946f, 0.131299f, 0.0f, 0.0f, 2.39016f, + 0.125435f, 0.0f, 0.0f, 2.3908f, 0.11982f, + 0.0f, 0.0f, 2.39139f, 0.114445f, 0.0f, + 0.0f, 2.39192f, 0.1093f, 0.0f, 0.0f, + 2.39241f, 0.104376f, 0.0f, 0.0f, 2.39286f, + 0.099666f, 0.0f, 0.0f, 2.39327f, 0.0951603f, + 0.0f, 0.0f, 2.39364f, 0.0908511f, 0.0f, + 0.0f, 2.39398f, 0.0867305f, 0.0f, 0.0f, + 2.39429f, 0.082791f, 0.0f, 0.0f, 2.39457f, + 0.0790251f, 0.0f, 0.0f, 2.39483f, 0.0754256f, + 0.0f, 0.0f, 2.39506f, 0.0719858f, 0.0f, + 0.0f, 2.39528f, 0.0686988f, 0.0f, 0.0f, + 2.39547f, 0.0655584f, 0.0f, 0.0f, 2.39565f, + 0.0625583f, 0.0f, 0.0f, 2.39582f, 0.0596925f, + 0.0f, 0.0f, 2.39596f, 0.0569554f, 0.0f, + 0.0f, 2.3961f, 0.0543413f, 0.0f, 0.0f, + 2.39622f, 0.0518451f, 0.0f, 0.0f, 2.39633f, + 0.0494615f, 0.0f, 0.0f, 2.39644f, 0.0471857f, + 0.0f, 0.0f, 2.39653f, 0.0450131f, 0.0f, + 0.0f, 2.39661f, 0.0429389f, 0.0f, 0.0f, + 2.39669f, 0.0409591f, 0.0f, 0.0f, 2.39676f, + 0.0390693f, 0.0f, 0.0f, 2.39682f, 0.0372656f, + 0.0f, 0.0f, 2.39688f, 0.0355441f, 0.0f, + 0.0f, 2.39694f, 0.0339013f, 0.0f, 0.0f, + 2.39698f, 0.0323336f, 0.0f, 0.0f, 2.39703f, + 0.0308377f, 0.0f, 0.0f, 2.39707f, 0.0294103f, + 0.0f, 0.0f, 2.3971f, 0.0280484f, 0.0f, + 0.0f, 2.39714f, 0.0267489f, 0.0f, 0.0f, + 2.39717f, 0.0255092f, 0.0f, 0.0f, 2.39719f, + 0.0243265f, 0.0f, 0.0f, 2.39722f, 0.0231982f, + 0.0f, 0.0f, 2.39724f, 0.0221218f, 0.0f, + 0.0f, 2.39726f, 0.0210951f, 0.0f, 0.0f, + 2.39728f, 0.0201157f, 0.0f, 0.0f, 2.3973f, + 0.0191815f, 0.0f, 0.0f, 2.39731f, 0.0182904f, + 0.0f, 0.0f, 2.39733f, 0.0174405f, 0.0f, + 0.0f, 2.39734f, 0.0166299f, 0.0f, 0.0f, + 2.39735f, 0.0158567f, 0.0f, 0.0f, 2.39736f, + 0.0151194f, 0.0f, 0.0f, 2.39737f, 0.0144161f, + 0.0f, 0.0f, 2.39738f, 0.0137455f, 0.0f, + 0.0f, 2.39739f, 0.0131059f, 0.0f, 0.0f, + 2.3974f, 0.0124959f, 0.0f, 0.0f, 2.3974f, + 0.0119143f, 0.0f, 0.0f, 2.39741f, 0.0113596f, + 0.0f, 0.0f, 2.39741f, 0.0108306f, 0.0f, + 0.0f, 2.39742f, 0.0103262f, 0.0f, 0.0f, + 2.39742f, 0.00984523f, 0.0f, 0.0f, 2.39743f, + 0.00938658f, 0.0f, 0.0f, 2.39743f, 0.00894924f, + 0.0f, 0.0f, 2.39744f, 0.00853223f, 0.0f, + 0.0f, 2.39744f, 0.00813459f, 0.0f, 0.0f, + 2.39744f, 0.00775545f, 0.0f, 0.0f, 2.39744f, + 0.00739393f, 0.0f, 0.0f, 2.39745f, 0.00704923f, + 0.0f, 0.0f, 2.39745f, 0.00672056f, 0.0f, + 0.0f, 2.39745f, 0.00640719f, 0.0f, 0.0f, + 2.39745f, 0.00610841f, 0.0f, 0.0f, 2.39745f, + 0.00582353f, 0.0f, 0.0f, 2.39745f, 0.00555191f, + 0.0f, 0.0f, 2.39746f, 0.00529295f, 0.0f, + 0.0f, 2.39746f, 0.00504604f, 0.0f, 0.0f, + 2.39746f, 0.00481063f, 0.0f, 0.0f, 2.39746f, + 0.00458619f, 0.0f, 0.0f, 2.39746f, 0.00437221f, + 0.0f, 0.0f, 2.39746f, 0.0041682f, 0.0f, + 0.0f, 2.39746f, 0.0039737f, 0.0f, 0.0f, + 2.39746f, 0.00378826f, 0.0f, 0.0f, 2.39746f, + 0.00361147f, 0.0f, 0.0f, 2.39746f, 0.00344291f, + 0.0f, 0.0f, 2.39746f, 0.00328222f, 0.0f, + 0.0f, 2.39746f, 0.00312902f, 0.0f, 0.0f, + 2.39746f, 0.00298297f, 0.0f, 0.0f, 2.39747f, + 0.00284372f, 0.0f, 0.0f, 2.39747f, 0.00271097f, + 0.0f, 0.0f, 2.39747f, 0.00258441f, 0.0f, + 0.0f, 2.39747f, 0.00246376f, 0.0f, 0.0f, + 2.39747f, 0.00234873f, 0.0f, 0.0f, 2.39747f, + 0.00223908f, 0.0f, 0.0f, 2.39747f, 0.00213453f, + 0.0f, 0.0f, 2.39747f, 0.00203487f, 0.0f, + 0.0f, 2.39747f, 0.00193986f, 0.0f, 0.0f, + 2.39747f, 0.00184928f, 0.0f, 0.0f, 2.39747f, + 0.00176292f, 0.0f, 0.0f, 2.39747f, 0.0016806f, + 0.0f, 0.0f, 2.39747f, 0.00160212f, 0.0f, + 0.0f, 2.39747f, 0.00152731f, 0.0f, 0.0f, + 2.39747f, 0.00145598f, 0.0f, 0.0f, 2.39747f, + 0.00138799f, 0.0f, 0.0f, 2.39747f, 0.00132316f, + 0.0f, 0.0f, 2.39747f, 0.00126137f, 0.0f, + 0.0f, 2.39747f, 0.00120246f, 0.0f, 0.0f, + 2.39747f, 0.0011463f, 0.0f, 0.0f, 2.39747f, + 0.00109276f, 0.0f, 0.0f, 2.39747f, 0.00104172f, + 0.0f, 0.0f, 2.39747f, 0.000993069f, 0.0f, + 0.0f, 2.39747f, 0.000946686f, 0.0f, 0.0f, + 2.39747f, 0.000902469f, 0.0f, 0.0f, 2.39747f, + 0.000860316f, 0.0f, 0.0f, 2.39747f, 0.000820132f, + 0.0f, 0.0f, 2.39747f, 0.000781825f, 0.0f, + 0.0f, 2.39747f, 0.000745306f, 0.0f, 0.0f, + 2.39747f, 0.000710492f, 0.0f, 0.0f, 2.39747f, + 0.000677305f, 0.0f, 0.0f, 2.39747f, 0.000645667f, + 0.0f, 0.0f, 2.39747f, 0.000615507f, 0.0f, + 0.0f, 2.39747f, 0.000586756f, 0.0f, 0.0f, + 2.39747f, 0.000559347f, 0.0f, 0.0f, 2.39747f, + 0.000533218f, 0.0f, 0.0f, 2.39747f, 0.00050831f, + 0.0f, 0.0f, 2.39747f, 0.000484565f, 0.0f, + 0.0f, 2.39747f, 0.000461929f, 0.0f, 0.0f, + 2.39747f, 0.000440351f, 0.0f, 0.0f, 2.39747f, + 0.00041978f, 0.0f, 0.0f, 2.39747f, 0.00040017f, + 0.0f, 0.0f, 2.39747f, 0.000381476f, 0.0f, + 0.0f, 2.39747f, 0.000363656f, 0.0f, 0.0f, + 2.39747f, 0.000346667f, 0.0f, 0.0f, 2.39747f, + 0.000330473f, 0.0f, 0.0f, 2.39747f, 0.000315034f, + 0.0f, 0.0f, 2.39747f, 0.000300317f, 0.0f, + 0.0f, 2.39747f, 0.000286287f, 0.0f, 0.0f, + 2.39747f, 0.000272913f, 0.0f, 0.0f, 2.39747f, + 0.000260164f, 0.0f, 0.0f, 2.39747f, 0.00024801f, + 0.0f, 0.0f, 2.39747f, 0.000236423f, 0.0f, + 0.0f, 2.39747f, 0.000225378f, 0.0f, 0.0f, + 2.39747f, 0.000214849f, 0.0f, 0.0f, 2.39747f, + 0.000204812f, 0.0f, 0.0f, 2.39747f, 0.000195244f, + 0.0f, 0.0f, 2.39747f, 0.000186122f, 0.0f, + 0.0f, 2.39747f, 0.000177427f, 0.0f, 0.0f, + 2.39747f, 0.000169138f, 0.0f, 0.0f, 2.39747f, + 0.000161236f, 0.0f, 0.0f, 2.39747f, 0.000153704f, + 0.0f, 0.0f, 2.39747f, 0.000146523f, 0.0f, + 0.0f, 2.39747f, 0.000139678f, 0.0f, 0.0f, + 2.39747f, 0.000133152f, 0.0f, 0.0f, 2.39747f, + 0.000126932f, 0.0f, 0.0f, 2.39747f, 0.000121001f, + 0.0f, 0.0f, 2.39747f, 0.000115348f, 0.0f, + 0.0f, 2.39747f, 0.00010996f, 0.0f, 0.0f, + 2.39747f, 0.000104822f, 0.0f, 0.0f, 2.39747f, + 9.99252e-05f, 0.0f, 0.0f, 2.39747f, 9.52568e-05f, + 0.0f, 0.0f, 2.39747f, 9.08065e-05f, 0.0f, + 0.0f, 2.39747f, 8.65641e-05f, 0.0f, 0.0f, + 2.39747f, 8.25199e-05f, 0.0f, 0.0f, 2.39747f, + 7.86646e-05f, 0.0f, 0.0f, 2.39747f, 7.49895e-05f, + 0.0f, 0.0f, 2.39747f, 7.1486e-05f, 0.0f, + 0.0f, 2.39747f, 6.81463e-05f, 0.0f, 0.0f, + 2.39747f, 6.49625e-05f, 0.0f, 0.0f, 2.39747f, + 6.19275e-05f, 0.0f, 0.0f, 2.39747f, 5.90343e-05f, + 0.0f, 0.0f, 2.39747f, 5.62762e-05f, 0.0f, + 0.0f, 2.39747f, 5.3647e-05f, 0.0f, 0.0f, + 2.39747f, 5.11407e-05f, 0.0f, 0.0f, 2.39747f, + 4.87514e-05f, 0.0f, 0.0f, 2.39747f, 4.64738e-05f, + 0.0f, 0.0f, 2.39747f, 4.43025e-05f, 0.0f, + 0.0f, 2.39747f, 4.22327e-05f, 0.0f, 0.0f, + 2.39747f, 4.02596e-05f, 0.0f, 0.0f, 2.39747f, + 3.83787e-05f, 0.0f, 0.0f, 2.39747f, 3.65857e-05f, + 0.0f, 0.0f, 2.39747f, 3.48764e-05f, 0.0f, + 0.0f, 2.39747f, 3.3247e-05f, 0.0f, 0.0f, + 2.39747f, 3.16937e-05f, 0.0f, 0.0f, 2.39747f, + 3.0213e-05f, 0.0f, 0.0f, 2.39747f, 2.88014e-05f, + 0.0f, 0.0f, 2.39747f, 2.74558e-05f, 0.0f, + 0.0f, 2.39747f, 2.61731e-05f, 0.0f, 0.0f, + 2.39747f, 2.49503e-05f, 0.0f, 0.0f, 2.39747f, + 2.37846e-05f, 0.0f, 0.0f, 2.39747f, 2.26734e-05f, + 0.0f, 0.0f, 2.39747f, 2.16141e-05f, 0.0f, + 0.0f, 2.39747f, 2.06043e-05f, 0.0f, 0.0f, + 2.39747f, 1.96417e-05f, 0.0f, 0.0f, 2.39747f, + 1.8724e-05f, 0.0f, 0.0f, 2.39747f, 1.78492e-05f, + 0.0f, 0.0f, 2.39747f, 1.70153e-05f, 0.0f, + 0.0f, 2.39747f, 1.62203e-05f, 0.0f, 0.0f, + 2.39747f, 1.54625e-05f, 0.0f, 0.0f, 2.39747f, + 1.47401e-05f, 0.0f, 0.0f, 2.39747f, 1.40515e-05f, + 0.0f, 0.0f, 2.39747f, 1.3395e-05f, 0.0f, + 0.0f, 2.39747f, 1.27692e-05f, 0.0f, 0.0f, + 2.39747f, 1.21726e-05f, 0.0f, 0.0f, 2.39747f, + 1.16039e-05f, 0.0f, 0.0f, 2.39747f, 1.10617e-05f, + 0.0f, 0.0f, 2.39747f, 1.05449e-05f, 0.0f, + 0.0f, 2.39747f, 1.00523e-05f, 0.0f, 0.0f, + 2.39747f, 9.58263e-06f, 0.0f, 0.0f, 2.39747f, + 9.13493e-06f, 0.0f, 0.0f, 2.39747f, 8.70814e-06f, + 0.0f, 0.0f, 2.39747f, 8.3013e-06f, 0.0f, + 0.0f, 2.39747f, 7.91346e-06f, 0.0f, 0.0f, + 2.39747f, 7.54374e-06f, 0.0f, 0.0f, 2.39747f, + 7.1913e-06f, 0.0f, 0.0f, 2.39747f, 6.85532e-06f, + 0.0f, 0.0f, 2.39747f, 6.53504e-06f, 0.0f, + 0.0f, 2.39747f, 6.22972e-06f, 0.0f, 0.0f, + 2.39747f, 5.93867e-06f, 0.0f, 0.0f, 2.39747f, + 5.66121e-06f, 0.0f, 0.0f, 2.39747f, 5.39672e-06f, + 0.0f, 0.0f, 2.39747f, 5.14458e-06f, 0.0f, + 0.0f, 2.39747f, 4.90423e-06f, 0.0f, 0.0f, + 2.39747f, 4.6751e-06f, 0.0f, 0.0f, 2.39747f, + 4.45668e-06f, 0.0f, 0.0f, 2.39747f, 4.24846e-06f, + 0.0f, 0.0f, 2.39747f, 4.04997e-06f, 0.0f, + 0.0f, 2.39747f, 3.86076e-06f, 0.0f, 0.0f, + 2.39747f, 3.68038e-06f, 0.0f, 0.0f, 2.39747f, + 3.50843e-06f, 0.0f, 0.0f, 2.39747f, 3.34452e-06f, + 0.0f, 0.0f, 2.39747f, 3.18826e-06f, 0.0f, + 0.0f, 2.39747f, 3.03931e-06f, 0.0f, 0.0f, + 2.39747f, 2.89731e-06f, 0.0f, 0.0f, 2.39747f, + 2.76195e-06f, 0.0f, 0.0f, 2.39747f, 2.63291e-06f, + 0.0f, 0.0f, 2.39747f, 2.5099e-06f, 0.0f, + 0.0f, 2.39747f, 2.39263e-06f, 0.0f, 0.0f, + 2.39747f, 2.28085e-06f, 0.0f, 0.0f, 2.39747f, + 2.17429e-06f, 0.0f, 0.0f, 2.39747f, 2.0727e-06f, + 0.0f, 0.0f, 2.39747f, 1.97587e-06f, 0.0f, + 0.0f, 2.39747f, 1.88355e-06f, 0.0f, 0.0f, + 2.39747f, 1.79555e-06f, 0.0f, 0.0f, 2.39747f, + 1.71167e-06f, 0.0f, 0.0f, 2.39747f, 1.6317e-06f, + 0.0f, 0.0f, 2.39747f, 1.55546e-06f, 0.0f, + 0.0f, 2.39747f, 1.48279e-06f, 0.0f, 0.0f, + 2.39747f, 1.41352e-06f, 0.0f, 0.0f, 2.39747f, + 1.34748e-06f, 0.0f, 0.0f}; + +} // namespace vraudio + +#endif // RESONANCE_AUDIO_AMBISONICS_AMBISONIC_SPREAD_COEFFICIENTS_H_ diff --git a/src/3rdparty/resonance-audio/resonance_audio/ambisonics/associated_legendre_polynomials_generator.cc b/src/3rdparty/resonance-audio/resonance_audio/ambisonics/associated_legendre_polynomials_generator.cc new file mode 100644 index 000000000..49b7abf89 --- /dev/null +++ b/src/3rdparty/resonance-audio/resonance_audio/ambisonics/associated_legendre_polynomials_generator.cc @@ -0,0 +1,147 @@ +/* +Copyright 2018 Google Inc. 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 "ambisonics/associated_legendre_polynomials_generator.h" + +#include <cmath> + +#include "base/logging.h" +#include "base/misc_math.h" + +namespace vraudio { + +AssociatedLegendrePolynomialsGenerator::AssociatedLegendrePolynomialsGenerator( + int max_degree, bool condon_shortley_phase, bool compute_negative_order) + : max_degree_(max_degree), + condon_shortley_phase_(condon_shortley_phase), + compute_negative_order_(compute_negative_order) { + DCHECK_GE(max_degree_, 0); +} + +std::vector<float> AssociatedLegendrePolynomialsGenerator::Generate( + float x) const { + std::vector<float> values(GetNumValues()); + + // Bases for the recurrence relations. + values[GetIndex(0, 0)] = ComputeValue(0, 0, x, values); + if (max_degree_ >= 1) values[GetIndex(1, 0)] = ComputeValue(1, 0, x, values); + + // Using recurrence relations, we now compute the rest of the values needed. + // (degree, 0), based on (degree - 1, 0) and (degree - 2, 0): + for (int degree = 2; degree <= max_degree_; ++degree) { + const int order = 0; + values[GetIndex(degree, order)] = ComputeValue(degree, order, x, values); + } + // (degree, degree): + for (int degree = 1; degree <= max_degree_; ++degree) { + const int order = degree; + values[GetIndex(degree, order)] = ComputeValue(degree, order, x, values); + } + // (degree, degree - 1): + for (int degree = 2; degree <= max_degree_; ++degree) { + const int order = degree - 1; + values[GetIndex(degree, order)] = ComputeValue(degree, order, x, values); + } + // The remaining positive orders, based on (degree - 1, order) and + // (degree - 2, order): + for (int degree = 3; degree <= max_degree_; ++degree) { + for (int order = 1; order <= degree - 2; ++order) { + values[GetIndex(degree, order)] = ComputeValue(degree, order, x, values); + } + } + // (degree, -order): + if (compute_negative_order_) { + for (int degree = 1; degree <= max_degree_; ++degree) { + for (int order = 1; order <= degree; ++order) { + values[GetIndex(degree, -order)] = + ComputeValue(degree, -order, x, values); + } + } + } + if (!condon_shortley_phase_) { + for (int degree = 1; degree <= max_degree_; ++degree) { + const int start_order = compute_negative_order_ ? -degree : 0; + for (int order = start_order; order <= degree; ++order) { + // Undo the Condon-Shortley phase. + values[GetIndex(degree, order)] *= + static_cast<float>(std::pow(-1, order)); + } + } + } + return values; +} + +size_t AssociatedLegendrePolynomialsGenerator::GetNumValues() const { + if (compute_negative_order_) + return (max_degree_ + 1) * (max_degree_ + 1); + else + return ((max_degree_ + 1) * (max_degree_ + 2)) / 2; +} + +size_t AssociatedLegendrePolynomialsGenerator::GetIndex(int degree, + int order) const { + CheckIndexValidity(degree, order); + size_t result; + if (compute_negative_order_) { + result = static_cast<size_t>(degree * (degree + 1) + order); + } else { + result = static_cast<size_t>((degree * (degree + 1)) / 2 + order); + } + DCHECK_GE(result, 0U); + DCHECK_LT(result, GetNumValues()); + return result; +} + +float AssociatedLegendrePolynomialsGenerator::ComputeValue( + int degree, int order, float x, const std::vector<float>& values) const { + CheckIndexValidity(degree, order); + if (degree == 0 && order == 0) { + return 1; + } else if (degree == 1 && order == 0) { + return x; + } else if (degree == order) { + return std::pow(-1.0f, static_cast<float>(degree)) * + DoubleFactorial(2 * degree - 1) * + std::pow((1.0f - x * x), 0.5f * static_cast<float>(degree)); + } else if (order == degree - 1) { + return x * static_cast<float>(2 * degree - 1) * + values[GetIndex(degree - 1, degree - 1)]; + } else if (order < 0) { + return std::pow(-1.0f, static_cast<float>(order)) * + Factorial(degree + order) / Factorial(degree - order) * + values[GetIndex(degree, -order)]; + } else { + return (static_cast<float>(2 * degree - 1) * x * + values[GetIndex(degree - 1, order)] - + static_cast<float>(degree - 1 + order) * + values[GetIndex(degree - 2, order)]) / + static_cast<float>(degree - order); + } +} + +void AssociatedLegendrePolynomialsGenerator::CheckIndexValidity( + int degree, int order) const { + DCHECK_GE(degree, 0); + DCHECK_LE(degree, max_degree_); + if (compute_negative_order_) { + DCHECK_LE(-degree, order); + } else { + DCHECK_GE(order, 0); + } + DCHECK_LE(order, degree); +} + +} // namespace vraudio diff --git a/src/3rdparty/resonance-audio/resonance_audio/ambisonics/associated_legendre_polynomials_generator.h b/src/3rdparty/resonance-audio/resonance_audio/ambisonics/associated_legendre_polynomials_generator.h new file mode 100644 index 000000000..84eb9d2f6 --- /dev/null +++ b/src/3rdparty/resonance-audio/resonance_audio/ambisonics/associated_legendre_polynomials_generator.h @@ -0,0 +1,89 @@ +/* +Copyright 2018 Google Inc. 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 RESONANCE_AUDIO_AMBISONICS_ASSOCIATED_LEGENDRE_POLYNOMIALS_GENERATOR_H_ +#define RESONANCE_AUDIO_AMBISONICS_ASSOCIATED_LEGENDRE_POLYNOMIALS_GENERATOR_H_ + +#include <stddef.h> +#include <vector> + +namespace vraudio { + +// Generates associated Legendre polynomials. +class AssociatedLegendrePolynomialsGenerator { + public: + // Constructs a generator for associated Legendre polynomials (ALP). + // + // @param max_degree The maximum ALP degree supported by this generator. + // @param condon_shortley_phase Whether the Condon-Shortley phase, (-1)^order, + // should be included in the polynomials generated. + // @param compute_negative_order Whether this generator should compute + // negative-ordered polynomials. + AssociatedLegendrePolynomialsGenerator(int max_degree, + bool condon_shortley_phase, + bool compute_negative_order); + + // Generates the associated Legendre polynomials for the given |x|, returning + // the computed sequence. + // + // @param x The abscissa (the polynomials' variable). + // @return Output vector of computed values. + std::vector<float> Generate(float x) const; + + // Gets the number of associated Legendre polynomials this generator produces. + // + // @return The number of associated Legendre polynomials this generator + // produces. + size_t GetNumValues() const; + + // Gets the index into the output vector for the given |degree| and |order|. + // + // @param degree The polynomial's degree. + // @param order The polynomial's order. + // @return The index into the vector of computed values corresponding to the + // specified ALP. + size_t GetIndex(int degree, int order) const; + + private: + // Computes the ALP for (degree, order) the given |x| using recurrence + // relations. It is assumed that the ALPs necessary for each computation are + // already computed and stored in |values|. + // + // @param degree The degree of the polynomial being computed. + // @param degree The order of the polynomial being computed. + // @param values The previously computed values. + // @return The computed polynomial. + inline float ComputeValue(int degree, int order, float x, + const std::vector<float>& values) const; + + // Checks the validity of the given index. + // + // @param degree The polynomial's degree. + // @param order The polynomial's order. + inline void CheckIndexValidity(int degree, int order) const; + + // The maximum polynomial degree that can be computed; must be >= 0. + int max_degree_ = 0; + // Whether the Condon-Shortley phase, (-1)^order, should be included in the + // polynomials generated. + bool condon_shortley_phase_ = false; + // Whether this generator should compute negative-ordered polynomials. + bool compute_negative_order_ = false; +}; + +} // namespace vraudio + +#endif // RESONANCE_AUDIO_AMBISONICS_ASSOCIATED_LEGENDRE_POLYNOMIALS_GENERATOR_H_ diff --git a/src/3rdparty/resonance-audio/resonance_audio/ambisonics/associated_legendre_polynomials_generator_test.cc b/src/3rdparty/resonance-audio/resonance_audio/ambisonics/associated_legendre_polynomials_generator_test.cc new file mode 100644 index 000000000..b9b6cdec4 --- /dev/null +++ b/src/3rdparty/resonance-audio/resonance_audio/ambisonics/associated_legendre_polynomials_generator_test.cc @@ -0,0 +1,167 @@ +/* +Copyright 2018 Google Inc. 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 "ambisonics/associated_legendre_polynomials_generator.h" + +#include <cmath> + +#include "third_party/googletest/googletest/include/gtest/gtest.h" +#include "base/constants_and_types.h" +#include "base/misc_math.h" + +namespace vraudio { + +namespace { + +// Tolerated test error margin for single-precision floating points. +const float kEpsilon = 1e-5f; + +// Generates associated Legendre polynomials up to and including the 4th degree, +// with the Condon-Shortley phase and negative orders included. +std::vector<float> GenerateExpectedValuesFourthDegree(float x) { + // From http://en.wikipedia.org/wiki/Associated_Legendre_polynomials + // Comments are (degree, order). + const std::vector<float> expected_values = { + 1.0f, // (0, 0) + 0.5f * std::sqrt(1.0f - x * x), // (1, -1) + x, // (1, 0) + -std::sqrt(1.0f - x * x), // (1, 1) + 1.0f / 8.0f * (1.0f - x * x), // (2, -2) + 0.5f * x * std::sqrt(1.0f - x * x), // (2, -1) + 0.5f * (3.0f * x * x - 1.0f), // (2, 0) + -3.0f * x * std::sqrt(1.0f - x * x), // (2, 1) + 3.0f * (1.0f - x * x), // (2, 2) + 15.0f / 720.0f * std::pow(1.0f - x * x, 3.0f / 2.0f), // (3, -3) + 15.0f / 120.0f * x * (1.0f - x * x), // (3, -2) + 3.0f / 24.0f * (5.0f * x * x - 1.0f) * + std::sqrt(1.0f - x * x), // (3, -1) + 0.5f * (5.0f * IntegerPow(x, 3) - 3.0f * x), // (3, 0) + -3.0f / 2.0f * (5.0f * x * x - 1.0f) * std::sqrt(1.0f - x * x), // (3, 1) + 15.0f * x * (1.0f - x * x), // (3, 2) + -15.0f * std::pow(1.0f - x * x, 3.0f / 2.0f), // (3, 3) + 105.0f / 40320.0f * IntegerPow(1.0f - x * x, 2), // (4, -4) + 105.0f / 5040.0f * x * std::pow(1.0f - x * x, 3.0f / 2.0f), // (4, -3) + 15.0f / 720.0f * (7.0f * x * x - 1.0f) * (1.0f - x * x), // (4, -2) + 5.0f / 40.0f * (7.0f * IntegerPow(x, 3) - 3.0f * x) * + std::sqrt(1.0f - x * x), // (4, -1) + 1.0f / 8.0f * + (35.0f * IntegerPow(x, 4) - 30.0f * x * x + 3.0f), // (4, 0) + -5.0f / 2.0f * (7.0f * IntegerPow(x, 3) - 3.0f * x) * + std::sqrt(1.0f - x * x), // (4, 1) + 15.0f / 2.0f * (7.0f * x * x - 1.0f) * (1.0f - x * x), // (4, 2) + -105.0f * x * std::pow(1.0f - x * x, 3.0f / 2.0f), // (4, 3) + 105.0f * IntegerPow(1.0f - x * x, 2) // (4, 4) + }; + + return expected_values; +} + +// Tests that the values given by GetIndex are successive indices (n, n+1, n+2, +// and so on). +TEST(AssociatedLegendrePolynomialsGeneratorTest, GetIndex_SuccessiveIndices) { + const int kMaxDegree = 5; + const AssociatedLegendrePolynomialsGenerator alp_generator(kMaxDegree, false, + true); + int last_index = -1; + for (int degree = 0; degree <= kMaxDegree; ++degree) { + for (int order = -degree; order <= degree; ++order) { + int index = static_cast<int>(alp_generator.GetIndex(degree, order)); + EXPECT_EQ(last_index + 1, index); + last_index = index; + } + } +} + +// Tests that the zeroth-degree, zeroth-order ALP is always 1. +TEST(AssociatedLegendrePolynomialsGeneratorTest, Generate_ZerothElementIsOne) { + const int kMaxDegree = 10; + for (int max_degree = 0; max_degree <= kMaxDegree; ++max_degree) { + for (int condon_shortley_phase = 0; condon_shortley_phase <= 1; + ++condon_shortley_phase) { + for (int compute_negative_order = 0; compute_negative_order <= 1; + ++compute_negative_order) { + AssociatedLegendrePolynomialsGenerator alp_generator( + max_degree, condon_shortley_phase != 0, + compute_negative_order != 0); + + const float kVariableStep = 0.2f; + for (float x = -1.0f; x <= 1.0f; x += kVariableStep) { + const std::vector<float> values = alp_generator.Generate(x); + + EXPECT_NEAR(values[0], 1.0f, kEpsilon); + } + } + } + } +} + +// Tests that the polynomials generated are correct until the 4th degree. +TEST(AssociatedLegendrePolynomialsGeneratorTest, Generate_CorrectFourthDegree) { + const int kMaxDegree = 4; + const bool kCondonShortleyPhase = true; + const bool kComputeNegativeOrder = true; + const AssociatedLegendrePolynomialsGenerator alp_generator( + kMaxDegree, kCondonShortleyPhase, kComputeNegativeOrder); + + const float kVariableStep = 0.05f; + for (float x = -1.0f; x <= 1.0f; x += kVariableStep) { + const std::vector<float> generated_values = alp_generator.Generate(x); + const std::vector<float> expected_values = + GenerateExpectedValuesFourthDegree(x); + ASSERT_EQ(expected_values.size(), generated_values.size()); + for (size_t i = 0; i < expected_values.size(); ++i) { + EXPECT_NEAR(generated_values[i], expected_values[i], kEpsilon) + << " at index " << i; + } + } +} + +// Tests that the Condon-Shortley phase is correctly applied. +TEST(AssociatedLegendrePolynomialsGeneratorTest, Generate_CondonShortleyPhase) { + const int kMaxDegree = 10; + const float kValue = 0.12345f; + for (int max_degree = 0; max_degree <= kMaxDegree; ++max_degree) { + for (int compute_negative_order = 0; compute_negative_order <= 1; + ++compute_negative_order) { + const AssociatedLegendrePolynomialsGenerator alp_generator_without_phase( + max_degree, false, compute_negative_order != 0); + const std::vector<float> values_without_phase = + alp_generator_without_phase.Generate(kValue); + + const AssociatedLegendrePolynomialsGenerator alp_generator_with_phase( + max_degree, true, compute_negative_order != 0); + const std::vector<float> values_with_phase = + alp_generator_with_phase.Generate(kValue); + + ASSERT_EQ(values_with_phase.size(), values_without_phase.size()); + for (int degree = 0; degree <= max_degree; ++degree) { + const int start_order = compute_negative_order ? -degree : 0; + for (int order = start_order; order <= degree; ++order) { + const size_t index = + alp_generator_without_phase.GetIndex(degree, order); + const float expected = values_without_phase[index] * + std::pow(-1.0f, static_cast<float>(order)); + EXPECT_NEAR(values_with_phase[index], expected, kEpsilon) + << " at degree " << degree << " and order " << order; + } + } + } + } +} + +} // namespace + +} // namespace vraudio diff --git a/src/3rdparty/resonance-audio/resonance_audio/ambisonics/foa_rotator.cc b/src/3rdparty/resonance-audio/resonance_audio/ambisonics/foa_rotator.cc new file mode 100644 index 000000000..f1ac68118 --- /dev/null +++ b/src/3rdparty/resonance-audio/resonance_audio/ambisonics/foa_rotator.cc @@ -0,0 +1,117 @@ +/* +Copyright 2018 Google Inc. 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. +*/ + +#if defined(_WIN32) +#define _SCL_SECURE_NO_WARNINGS +#endif + +#include "ambisonics/foa_rotator.h" + +#include <algorithm> + +#include "base/constants_and_types.h" +#include "base/misc_math.h" + +namespace vraudio { + +bool FoaRotator::Process(const WorldRotation& target_rotation, + const AudioBuffer& input, AudioBuffer* output) { + + DCHECK(output); + DCHECK_EQ(input.num_channels(), kNumFirstOrderAmbisonicChannels); + DCHECK_EQ(input.num_channels(), output->num_channels()); + DCHECK_EQ(input.num_frames(), output->num_frames()); + + static const WorldRotation kIdentityRotation; + + if (current_rotation_.AngularDifferenceRad(kIdentityRotation) < + kRotationQuantizationRad && + target_rotation.AngularDifferenceRad(kIdentityRotation) < + kRotationQuantizationRad) { + return false; + } + + if (current_rotation_.AngularDifferenceRad(target_rotation) < + kRotationQuantizationRad) { + // Rotate the whole input buffer frame by frame. + Rotate(current_rotation_, 0, input.num_frames(), input, output); + return true; + } + + // In order to perform a smooth rotation, we divide the buffer into + // chunks of size |kSlerpFrameInterval|. + // + + const size_t kSlerpFrameInterval = 32; + + WorldRotation slerped_rotation; + // Rotate the input buffer at every slerp update interval. Truncate the + // final chunk if the input buffer is not an integer multiple of the + // chunk size. + for (size_t i = 0; i < input.num_frames(); i += kSlerpFrameInterval) { + const size_t duration = + std::min(input.num_frames() - i, kSlerpFrameInterval); + const float interpolation_factor = static_cast<float>(i + duration) / + static_cast<float>(input.num_frames()); + slerped_rotation = + current_rotation_.slerp(interpolation_factor, target_rotation); + // Rotate the input buffer frame by frame within the current chunk. + Rotate(slerped_rotation, i, duration, input, output); + } + + current_rotation_ = target_rotation; + + return true; +} + +void FoaRotator::Rotate(const WorldRotation& target_rotation, + size_t start_location, size_t duration, + const AudioBuffer& input, AudioBuffer* output) { + + const AudioBuffer::Channel& input_channel_audio_space_w = input[0]; + const AudioBuffer::Channel& input_channel_audio_space_y = input[1]; + const AudioBuffer::Channel& input_channel_audio_space_z = input[2]; + const AudioBuffer::Channel& input_channel_audio_space_x = input[3]; + AudioBuffer::Channel* output_channel_audio_space_w = &(*output)[0]; + AudioBuffer::Channel* output_channel_audio_space_y = &(*output)[1]; + AudioBuffer::Channel* output_channel_audio_space_z = &(*output)[2]; + AudioBuffer::Channel* output_channel_audio_space_x = &(*output)[3]; + + for (size_t frame = start_location; frame < start_location + duration; + ++frame) { + // Convert the current audio frame into world space position. + temp_audio_position_(0) = input_channel_audio_space_x[frame]; + temp_audio_position_(1) = input_channel_audio_space_y[frame]; + temp_audio_position_(2) = input_channel_audio_space_z[frame]; + ConvertWorldFromAudioPosition(temp_audio_position_, &temp_world_position_); + // Apply rotation to |world_position| and return to audio space. + temp_rotated_world_position_ = target_rotation * temp_world_position_; + + ConvertAudioFromWorldPosition(temp_rotated_world_position_, + &temp_rotated_audio_position_); + (*output_channel_audio_space_x)[frame] = + temp_rotated_audio_position_(0); // X + (*output_channel_audio_space_y)[frame] = + temp_rotated_audio_position_(1); // Y + (*output_channel_audio_space_z)[frame] = + temp_rotated_audio_position_(2); // Z + } + // Copy W channel. + std::copy_n(&input_channel_audio_space_w[start_location], duration, + &(*output_channel_audio_space_w)[start_location]); +} + +} // namespace vraudio diff --git a/src/3rdparty/resonance-audio/resonance_audio/ambisonics/foa_rotator.h b/src/3rdparty/resonance-audio/resonance_audio/ambisonics/foa_rotator.h new file mode 100644 index 000000000..a91b975d4 --- /dev/null +++ b/src/3rdparty/resonance-audio/resonance_audio/ambisonics/foa_rotator.h @@ -0,0 +1,61 @@ +/* +Copyright 2018 Google Inc. 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 RESONANCE_AUDIO_AMBISONICS_FOA_ROTATOR_H_ +#define RESONANCE_AUDIO_AMBISONICS_FOA_ROTATOR_H_ + +#include "base/audio_buffer.h" +#include "base/spherical_angle.h" + +namespace vraudio { + +// Rotator for first order ambisonic soundfields. It supports ACN channel +// ordering and SN3D normalization (AmbiX). +class FoaRotator { + public: + // @param target_rotation Target rotation to be applied to the input buffer. + // @param input First order soundfield input buffer to be rotated. + // @param output Pointer to output buffer. + // @return True if rotation has been applied. + bool Process(const WorldRotation& target_rotation, const AudioBuffer& input, + AudioBuffer* output); + + private: + // Method which rotates a specified chunk of data in the AudioBuffer. + // + // @param target_rotation Target rotation to be applied to the soundfield. + // @param start_location Sample index in the soundfield where the rotation + // should begin. + // @param duration Number of samples in soundfield to be rotated. + // @param input First order soundfield input buffer to be rotated. + // @param output Pointer to output buffer. + void Rotate(const WorldRotation& target_rotation, size_t start_location, + size_t duration, const AudioBuffer& input, AudioBuffer* output); + + // Current rotation which is used in the interpolation process in order to + // perform a smooth rotation. Initialized with an identity matrix. + WorldRotation current_rotation_; + + // Preallocation of temporary variables used during rotation. + AudioPosition temp_audio_position_; + WorldPosition temp_world_position_; + AudioPosition temp_rotated_audio_position_; + WorldPosition temp_rotated_world_position_; +}; + +} // namespace vraudio + +#endif // RESONANCE_AUDIO_AMBISONICS_FOA_ROTATOR_H_ diff --git a/src/3rdparty/resonance-audio/resonance_audio/ambisonics/foa_rotator_test.cc b/src/3rdparty/resonance-audio/resonance_audio/ambisonics/foa_rotator_test.cc new file mode 100644 index 000000000..e2f957fa2 --- /dev/null +++ b/src/3rdparty/resonance-audio/resonance_audio/ambisonics/foa_rotator_test.cc @@ -0,0 +1,201 @@ +/* +Copyright 2018 Google Inc. 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 "ambisonics/foa_rotator.h" + +#include <algorithm> + +#include "third_party/googletest/googletest/include/gtest/gtest.h" +#include "ambisonics/ambisonic_codec_impl.h" +#include "base/constants_and_types.h" +#include "utils/planar_interleaved_conversion.h" + +namespace vraudio { + +namespace { + +using testing::tuple; +using testing::Values; + +const int kAmbisonicOrder = 1; + +// Rotation angle used in the test. +const float kAngleDegrees = 90.0f; + +// Initial, arbitrary direction of the encoded soundfield source. +const SphericalAngle kInitialSourceAngle = + SphericalAngle::FromDegrees(22.0f, 33.0f); + +// Expected direction of the sound source rotated by |kAngleDegrees| against the +// right facing axis (x in the world coordinate system). +const SphericalAngle kXrotatedSourceAngle = + SphericalAngle::FromDegrees(150.021778639249f, 51.0415207997462f); + +// Expected direction of the sound source rotated by |kAngleDegrees| against the +// upward facing axis (y in the world coordinate system). +const SphericalAngle kYrotatedSourceAngle = + SphericalAngle::FromDegrees(112.0f, 33.0f); + +// Expected direction of the sound source rotated by |kAngleDegrees| against the +// back facing axis (z in the world coordinate system). +const SphericalAngle kZrotatedSourceAngle = + SphericalAngle::FromDegrees(35.0077312770459f, -18.3108067341351f); + +// Rotation interpolation interval in terms of frames (used by the FoaRotator). +const size_t kSlerpFrameInterval = 32; + +class FoaRotatorTest : public ::testing::Test { + protected: + FoaRotatorTest() {} + + // Virtual methods from ::testing::Test + ~FoaRotatorTest() override {} + + void SetUp() override { + // Initialize the first order soundfield rotator. + foa_rotator_ = std::unique_ptr<FoaRotator>(new FoaRotator); + } + + void TearDown() override {} + + // Test method which creates a soundfield buffer, rotates it by + // |rotation_angle| against |rotation_axis|, and compares the output to the + // reference soundfield buffer (where the source is spatialized at the + // |expected_source_angle|). + void CompareRotatedAndReferenceSoundfields( + const std::vector<float>& input_data, float rotation_angle, + const WorldPosition& rotation_axis, + const SphericalAngle& expected_source_angle) { + AudioBuffer input_buffer(1, input_data.size()); + FillAudioBuffer(input_data, 1, &input_buffer); + // Initialize a first order mono codec with the |kInitialSourceAngle|. This + // will be used to obtain the rotated soundfield. + std::unique_ptr<MonoAmbisonicCodec<>> rotated_source_mono_codec( + new MonoAmbisonicCodec<>(kAmbisonicOrder, {kInitialSourceAngle})); + // Initialize a first order mono codec with the expected source angle. This + // will be used as a reference for the rotated soundfield. + std::unique_ptr<MonoAmbisonicCodec<>> reference_source_mono_codec( + new MonoAmbisonicCodec<>(kAmbisonicOrder, {expected_source_angle})); + // Generate soundfield buffers representing sound sources at required + // angles. + AudioBuffer encoded_rotated_buffer(kNumFirstOrderAmbisonicChannels, + input_data.size()); + AudioBuffer encoded_reference_buffer(kNumFirstOrderAmbisonicChannels, + input_data.size()); + rotated_source_mono_codec->EncodeBuffer(input_buffer, + &encoded_rotated_buffer); + reference_source_mono_codec->EncodeBuffer(input_buffer, + &encoded_reference_buffer); + // Rotate the test soundfield by |rotation_angle| degrees wrt the given + // |rotation_axis|. + const WorldRotation rotation = WorldRotation( + AngleAxisf(rotation_angle * kRadiansFromDegrees, rotation_axis)); + const bool result = foa_rotator_->Process(rotation, encoded_rotated_buffer, + &encoded_rotated_buffer); + EXPECT_TRUE(result); + // Check if the sound source in the reference and rotated buffers are in the + // same direction. + // If the buffer size is more than |kSlerpFrameInterval|, due to + // interpolation, we expect that the last |kSlerpFrameInterval| frames have + // reached the target rotation. + // If the buffer size is less than |kSlerpFrameInterval|, because no + // interpolation is applied, the rotated soundfield should result from + // frame 0. + for (size_t channel = 0; channel < encoded_rotated_buffer.num_channels(); + ++channel) { + const int num_frames = + static_cast<int>(encoded_rotated_buffer[channel].size()); + const int interval = static_cast<int>(kSlerpFrameInterval); + const int frames_to_compare = + (num_frames % interval) ? (num_frames % interval) : interval; + const int start_frame = std::max(0, num_frames - frames_to_compare); + ASSERT_LT(start_frame, num_frames); + for (int frame = start_frame; frame < num_frames; ++frame) { + EXPECT_NEAR(encoded_rotated_buffer[channel][frame], + encoded_reference_buffer[channel][frame], kEpsilonFloat); + } + } + } + + // First Order Ambisonic rotator instance. + std::unique_ptr<FoaRotator> foa_rotator_; +}; + +// Tests that no rotation is aplied if |kRotationQuantizationRad| is not +// exceeded. +TEST_F(FoaRotatorTest, RotationThresholdTest) { + const size_t kFramesPerBuffer = 16; + const std::vector<float> kInputVector( + kFramesPerBuffer * kNumFirstOrderAmbisonicChannels, 1.0f); + AudioBuffer input_buffer(kNumFirstOrderAmbisonicChannels, kFramesPerBuffer); + FillAudioBuffer(kInputVector, kNumFirstOrderAmbisonicChannels, &input_buffer); + AudioBuffer output_buffer(kNumFirstOrderAmbisonicChannels, kFramesPerBuffer); + const WorldRotation kSmallRotation = + WorldRotation(1.0f, 0.001f, 0.001f, 0.001f); + const WorldRotation kLargeRotation = WorldRotation(1.0f, 0.1f, 0.1f, 0.1f); + EXPECT_FALSE( + foa_rotator_->Process(kSmallRotation, input_buffer, &output_buffer)); + EXPECT_TRUE( + foa_rotator_->Process(kLargeRotation, input_buffer, &output_buffer)); +} + +typedef tuple<WorldPosition, SphericalAngle> TestParams; +class FoaAxesRotationTest : public FoaRotatorTest, + public ::testing::WithParamInterface<TestParams> {}; + +// Tests first order soundfield rotation against the x, y and z axis using long +// input buffers (>> slerp interval). +TEST_P(FoaAxesRotationTest, CompareWithExpectedAngleLongBuffer) { + const WorldPosition& rotation_axis = ::testing::get<0>(GetParam()); + const SphericalAngle& expected_angle = ::testing::get<1>(GetParam()); + const size_t kLongFramesPerBuffer = 512; + const std::vector<float> kLongInputData(kLongFramesPerBuffer, 1.0f); + CompareRotatedAndReferenceSoundfields(kLongInputData, kAngleDegrees, + rotation_axis, expected_angle); +} + +// Tests first order soundfield rotation against the x, y and z axes using short +// input buffers (< slerp interval). +TEST_P(FoaAxesRotationTest, CompareWithExpectedAngleShortBuffer) { + const WorldPosition& rotation_axis = ::testing::get<0>(GetParam()); + const SphericalAngle& expected_angle = ::testing::get<1>(GetParam()); + const size_t kShortFramesPerBuffer = kSlerpFrameInterval / 2; + const std::vector<float> kShortInputData(kShortFramesPerBuffer, 1.0f); + CompareRotatedAndReferenceSoundfields(kShortInputData, kAngleDegrees, + rotation_axis, expected_angle); +} + +// Tests first order soundfield rotation against the x, y and z axes using +// buffer sizes that are (> slerp interval) and not integer multiples of +// the slerp interval. +TEST_P(FoaAxesRotationTest, CompareWithExpectedAngleOddBufferSizes) { + const WorldPosition& rotation_axis = ::testing::get<0>(GetParam()); + const SphericalAngle& expected_angle = ::testing::get<1>(GetParam()); + const size_t frames_per_buffer = kSlerpFrameInterval + 3U; + const std::vector<float> kShortInputData(frames_per_buffer, 1.0f); + CompareRotatedAndReferenceSoundfields(kShortInputData, kAngleDegrees, + rotation_axis, expected_angle); +} + +INSTANTIATE_TEST_CASE_P( + TestParameters, FoaAxesRotationTest, + Values(TestParams({1.0f, 0.0f, 0.0f}, kXrotatedSourceAngle), + TestParams({0.0f, 1.0f, 0.0f}, kYrotatedSourceAngle), + TestParams({0.0f, 0.0f, 1.0f}, kZrotatedSourceAngle))); + +} // namespace + +} // namespace vraudio diff --git a/src/3rdparty/resonance-audio/resonance_audio/ambisonics/hoa_rotator.cc b/src/3rdparty/resonance-audio/resonance_audio/ambisonics/hoa_rotator.cc new file mode 100644 index 000000000..a76575e22 --- /dev/null +++ b/src/3rdparty/resonance-audio/resonance_audio/ambisonics/hoa_rotator.cc @@ -0,0 +1,308 @@ +/* +Copyright 2018 Google Inc. 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 "ambisonics/hoa_rotator.h" + +#include <algorithm> +#include <cmath> + +#include "ambisonics/utils.h" +#include "base/constants_and_types.h" +#include "base/logging.h" +#include "base/misc_math.h" + +namespace vraudio { + +namespace { + +// Below are the helper methods to compute SH rotation using recursion. The code +// is branched / modified from: + +// maths described in the following papers: +// +// [1] R. Green, "Spherical Harmonic Lighting: The Gritty Details", GDC 2003, +// http://www.research.scea.com/gdc2003/spherical-harmonic-lighting.pdf +// [2] J. Ivanic and K. Ruedenberg, "Rotation Matrices for Real Spherical +// Harmonics. Direct Determination by Recursion", J. Phys. Chem., vol. 100, +// no. 15, pp. 6342-6347, 1996. +// http://pubs.acs.org/doi/pdf/10.1021/jp953350u +// [2b] Corrections to initial publication: +// http://pubs.acs.org/doi/pdf/10.1021/jp9833350 + +// Kronecker Delta function. +inline float KroneckerDelta(int i, int j) { return (i == j) ? 1.0f : 0.0f; } + +// [2] uses an odd convention of referring to the rows and columns using +// centered indices, so the middle row and column are (0, 0) and the upper +// left would have negative coordinates. +// +// This is a convenience function to allow us to access an Eigen::MatrixXf +// in the same manner, assuming r is a (2l+1)x(2l+1) matrix. +float GetCenteredElement(const Eigen::MatrixXf& r, int i, int j) { + // The shift to go from [-l, l] to [0, 2l] is (rows - 1) / 2 = l, + // (since the matrix is assumed to be square, rows == cols). + const int offset = (static_cast<int>(r.rows()) - 1) / 2; + return r(i + offset, j + offset); +} + +// Helper function defined in [2] that is used by the functions U, V, W. +// This should not be called on its own, as U, V, and W (and their coefficients) +// select the appropriate matrix elements to access arguments |a| and |b|. +float P(int i, int a, int b, int l, const std::vector<Eigen::MatrixXf>& r) { + if (b == l) { + return GetCenteredElement(r[1], i, 1) * + GetCenteredElement(r[l - 1], a, l - 1) - + GetCenteredElement(r[1], i, -1) * + GetCenteredElement(r[l - 1], a, -l + 1); + } else if (b == -l) { + return GetCenteredElement(r[1], i, 1) * + GetCenteredElement(r[l - 1], a, -l + 1) + + GetCenteredElement(r[1], i, -1) * + GetCenteredElement(r[l - 1], a, l - 1); + } else { + return GetCenteredElement(r[1], i, 0) * GetCenteredElement(r[l - 1], a, b); + } +} + +// The functions U, V, and W should only be called if the correspondingly +// named coefficient u, v, w from the function ComputeUVWCoeff() is non-zero. +// When the coefficient is 0, these would attempt to access matrix elements that +// are out of bounds. The vector of rotations, |r|, must have the |l - 1| +// previously completed band rotations. These functions are valid for |l >= 2|. + +float U(int m, int n, int l, const std::vector<Eigen::MatrixXf>& r) { + // Although [1, 2] split U into three cases for m == 0, m < 0, m > 0 + // the actual values are the same for all three cases. + return P(0, m, n, l, r); +} + +float V(int m, int n, int l, const std::vector<Eigen::MatrixXf>& r) { + if (m == 0) { + return P(1, 1, n, l, r) + P(-1, -1, n, l, r); + } else if (m > 0) { + const float d = KroneckerDelta(m, 1); + return P(1, m - 1, n, l, r) * std::sqrt(1 + d) - + P(-1, -m + 1, n, l, r) * (1 - d); + } else { + // Note there is apparent errata in [1,2,2b] dealing with this particular + // case. [2b] writes it should be P*(1-d)+P*(1-d)^0.5 + // [1] writes it as P*(1+d)+P*(1-d)^0.5, but going through the math by hand, + // you must have it as P*(1-d)+P*(1+d)^0.5 to form a 2^.5 term, which + // parallels the case where m > 0. + const float d = KroneckerDelta(m, -1); + return P(1, m + 1, n, l, r) * (1 - d) + + P(-1, -m - 1, n, l, r) * std::sqrt(1 + d); + } +} + +float W(int m, int n, int l, const std::vector<Eigen::MatrixXf>& r) { + if (m == 0) { + // Whenever this happens, w is also 0 so W can be anything. + return 0.0f; + } else if (m > 0) { + return P(1, m + 1, n, l, r) + P(-1, -m - 1, n, l, r); + } else { + return P(1, m - 1, n, l, r) - P(-1, -m + 1, n, l, r); + } +} + +// Calculates the coefficients applied to the U, V, and W functions. Because +// their equations share many common terms they are computed simultaneously. +void ComputeUVWCoeff(int m, int n, int l, float* u, float* v, float* w) { + const float d = KroneckerDelta(m, 0); + const float denom = (abs(n) == l ? static_cast<float>(2 * l * (2 * l - 1)) + : static_cast<float>((l + n) * (l - n))); + const float one_over_denom = 1.0f / denom; + + *u = std::sqrt(static_cast<float>((l + m) * (l - m)) * one_over_denom); + *v = 0.5f * + std::sqrt((1.0f + d) * static_cast<float>(l + abs(m) - 1) * + (static_cast<float>(l + abs(m))) * one_over_denom) * + (1.0f - 2.0f * d); + *w = -0.5f * + std::sqrt(static_cast<float>(l - abs(m) - 1) * + (static_cast<float>(l - abs(m))) * one_over_denom) * + (1.0f - d); +} + +// Calculates the (2l+1)x(2l+1) rotation matrix for the band l. +// This uses the matrices computed for band 1 and band l-1 to compute the +// matrix for band l. |rotations| must contain the previously computed l-1 +// rotation matrices. +// +// This implementation comes from p. 5 (6346), Table 1 and 2 in [2] taking +// into account the corrections from [2b]. +void ComputeBandRotation(int l, std::vector<Eigen::MatrixXf>* rotations) { + // The lth band rotation matrix has rows and columns equal to the number of + // coefficients within that band (-l <= m <= l implies 2l + 1 coefficients). + Eigen::MatrixXf rotation(2 * l + 1, 2 * l + 1); + for (int m = -l; m <= l; ++m) { + for (int n = -l; n <= l; ++n) { + float u, v, w; + ComputeUVWCoeff(m, n, l, &u, &v, &w); + + // The functions U, V, W are only safe to call if the coefficients + // u, v, w are not zero. + if (std::abs(u) > 0.0f) u *= U(m, n, l, *rotations); + if (std::abs(v) > 0.0f) v *= V(m, n, l, *rotations); + if (std::abs(w) > 0.0f) w *= W(m, n, l, *rotations); + + rotation(m + l, n + l) = (u + v + w); + } + } + (*rotations)[l] = rotation; +} + +} // namespace + +HoaRotator::HoaRotator(int ambisonic_order) + : ambisonic_order_(ambisonic_order), + rotation_matrices_(ambisonic_order_ + 1), + rotation_matrix_( + static_cast<int>(GetNumPeriphonicComponents(ambisonic_order)), + static_cast<int>(GetNumPeriphonicComponents(ambisonic_order))) { + DCHECK_GE(ambisonic_order_, 2); + + // Initialize rotation sub-matrices. + // Order 0 matrix (first band) is simply the 1x1 identity. + Eigen::MatrixXf r(1, 1); + r(0, 0) = 1.0f; + rotation_matrices_[0] = r; + // All the other ambisonic orders (bands) are set to identity matrices of + // corresponding sizes. + for (int l = 1; l <= ambisonic_order_; ++l) { + const size_t submatrix_size = GetNumNthOrderPeriphonicComponents(l); + r.resize(submatrix_size, submatrix_size); + rotation_matrices_[l] = r.setIdentity(); + } + // Initialize the final rotation matrix to an identity matrix. + rotation_matrix_.setIdentity(); +} + +bool HoaRotator::Process(const WorldRotation& target_rotation, + const AudioBuffer& input, AudioBuffer* output) { + + DCHECK(output); + DCHECK_EQ(input.num_channels(), GetNumPeriphonicComponents(ambisonic_order_)); + DCHECK_EQ(input.num_channels(), output->num_channels()); + DCHECK_EQ(input.num_frames(), output->num_frames()); + + static const WorldRotation kIdentityRotation; + + if (current_rotation_.AngularDifferenceRad(kIdentityRotation) < + kRotationQuantizationRad && + target_rotation.AngularDifferenceRad(kIdentityRotation) < + kRotationQuantizationRad) { + return false; + } + + const size_t channel_stride = input.GetChannelStride(); + + typedef Eigen::Matrix<float, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor> + RowMajorMatrixf; + + const Eigen::Map<const RowMajorMatrixf, Eigen::Aligned, Eigen::OuterStride<>> + input_matrix(input[0].begin(), static_cast<int>(input.num_channels()), + static_cast<int>(input.num_frames()), + Eigen::OuterStride<>(static_cast<int>(channel_stride))); + + Eigen::Map<RowMajorMatrixf, Eigen::Aligned, Eigen::OuterStride<>> + output_matrix((*output)[0].begin(), + static_cast<int>(input.num_channels()), + static_cast<int>(input.num_frames()), + Eigen::OuterStride<>(static_cast<int>(channel_stride))); + + if (current_rotation_.AngularDifferenceRad(target_rotation) < + kRotationQuantizationRad) { + output_matrix = rotation_matrix_ * input_matrix; + return true; + } + + // In order to perform a smooth rotation, we divide the buffer into + // chunks of size |kSlerpFrameInterval|. + // + + const size_t kSlerpFrameInterval = 32; + + WorldRotation slerped_rotation; + // Rotate the input buffer at every slerp update interval. Truncate the + // final chunk if the input buffer is not an integer multiple of the + // chunk size. + for (size_t i = 0; i < input.num_frames(); i += kSlerpFrameInterval) { + const size_t duration = + std::min(input.num_frames() - i, kSlerpFrameInterval); + const float interpolation_factor = static_cast<float>(i + duration) / + static_cast<float>(input.num_frames()); + UpdateRotationMatrix( + current_rotation_.slerp(interpolation_factor, target_rotation)); + output_matrix.block(0 /* first channel */, i, output->num_channels(), + duration) = + rotation_matrix_ * input_matrix.block(0 /* first channel */, i, + input.num_channels(), duration); + } + current_rotation_ = target_rotation; + + return true; +} + +void HoaRotator::UpdateRotationMatrix(const WorldRotation& rotation) { + + + // There is no need to update 0th order 1-element sub-matrix. + // First order sub-matrix can be updated directly from the WorldRotation + // quaternion. However, we must account for the flipped left-right and + // front-back axis in the World coordinates. + AudioRotation rotation_audio_space; + ConvertAudioFromWorldRotation(rotation, &rotation_audio_space); + rotation_matrices_[1] = rotation_audio_space.toRotationMatrix(); + rotation_matrix_.block(1, 1, 3, 3) = rotation_matrices_[1]; + + // Sub-matrices for the remaining orders are updated recursively using the + // equations provided in [2, 2b]. An example final rotation matrix composed of + // sub-matrices of orders 0 to 3 has the following structure: + // + // X | 0 0 0 | 0 0 0 0 0 | 0 0 0 0 0 0 0 + // ------------------------------------- + // 0 | X X X | 0 0 0 0 0 | 0 0 0 0 0 0 0 + // 0 | X X X | 0 0 0 0 0 | 0 0 0 0 0 0 0 + // 0 | X X X | 0 0 0 0 0 | 0 0 0 0 0 0 0 + // ------------------------------------- + // 0 | 0 0 0 | X X X X X | 0 0 0 0 0 0 0 + // 0 | 0 0 0 | X X X X X | 0 0 0 0 0 0 0 + // 0 | 0 0 0 | X X X X X | 0 0 0 0 0 0 0 + // 0 | 0 0 0 | X X X X X | 0 0 0 0 0 0 0 + // 0 | 0 0 0 | X X X X X | 0 0 0 0 0 0 0 + // ------------------------------------- + // 0 | 0 0 0 | 0 0 0 0 0 | X X X X X X X + // 0 | 0 0 0 | 0 0 0 0 0 | X X X X X X X + // 0 | 0 0 0 | 0 0 0 0 0 | X X X X X X X + // 0 | 0 0 0 | 0 0 0 0 0 | X X X X X X X + // 0 | 0 0 0 | 0 0 0 0 0 | X X X X X X X + // 0 | 0 0 0 | 0 0 0 0 0 | X X X X X X X + // 0 | 0 0 0 | 0 0 0 0 0 | X X X X X X X + // + for (int current_order = 2; current_order <= ambisonic_order_; + ++current_order) { + ComputeBandRotation(current_order, &rotation_matrices_); + const int index = current_order * current_order; + const int size = + static_cast<int>(GetNumNthOrderPeriphonicComponents(current_order)); + rotation_matrix_.block(index, index, size, size) = + rotation_matrices_[current_order]; + } +} + +} // namespace vraudio diff --git a/src/3rdparty/resonance-audio/resonance_audio/ambisonics/hoa_rotator.h b/src/3rdparty/resonance-audio/resonance_audio/ambisonics/hoa_rotator.h new file mode 100644 index 000000000..956811df4 --- /dev/null +++ b/src/3rdparty/resonance-audio/resonance_audio/ambisonics/hoa_rotator.h @@ -0,0 +1,69 @@ +/* +Copyright 2018 Google Inc. 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 RESONANCE_AUDIO_AMBISONICS_HOA_ROTATOR_H_ +#define RESONANCE_AUDIO_AMBISONICS_HOA_ROTATOR_H_ + +#include <vector> + +#include "Eigen/Dense" +#include "base/audio_buffer.h" +#include "base/misc_math.h" + +namespace vraudio { + +// Rotator for higher order ambisonic sound fields. It supports ACN channel +// ordering and SN3D normalization (AmbiX). +class HoaRotator { + public: + // Constructs a sound field rotator of an arbitrary ambisonic order. + // + // @param ambisonic_order Order of ambisonic sound field. + explicit HoaRotator(int ambisonic_order); + + // Performs a smooth inplace rotation of a sound field buffer from + // |current_rotation_| to |target_rotation|. + // + // @param target_rotation Target rotation to be applied to the input buffer. + // @param input Higher order sound field input buffer to be rotated. + // @param output Pointer to output buffer. + // @return True if rotation has been applied. + bool Process(const WorldRotation& target_rotation, const AudioBuffer& input, + AudioBuffer* output); + + private: + // Updates the rotation matrix with using supplied WorldRotation. + // + // @param rotation World rotation. + void UpdateRotationMatrix(const WorldRotation& rotation); + + // Order of the ambisonic sound field handled by the rotator. + const int ambisonic_order_; + + // Current rotation which is used in the interpolation process in order to + // compute new rotation matrix. Initialized with an identity rotation. + WorldRotation current_rotation_; + + // Spherical harmonics rotation sub-matrices for each order. + std::vector<Eigen::MatrixXf> rotation_matrices_; + + // Final spherical harmonics rotation matrix. + Eigen::MatrixXf rotation_matrix_; +}; + +} // namespace vraudio + +#endif // RESONANCE_AUDIO_AMBISONICS_HOA_ROTATOR_H_ diff --git a/src/3rdparty/resonance-audio/resonance_audio/ambisonics/hoa_rotator_test.cc b/src/3rdparty/resonance-audio/resonance_audio/ambisonics/hoa_rotator_test.cc new file mode 100644 index 000000000..6472412c9 --- /dev/null +++ b/src/3rdparty/resonance-audio/resonance_audio/ambisonics/hoa_rotator_test.cc @@ -0,0 +1,201 @@ +/* +Copyright 2018 Google Inc. 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 "ambisonics/hoa_rotator.h" + +#include <algorithm> + +#include "third_party/googletest/googletest/include/gtest/gtest.h" +#include "ambisonics/ambisonic_codec_impl.h" +#include "base/constants_and_types.h" +#include "utils/planar_interleaved_conversion.h" + +namespace vraudio { + +namespace { + +using testing::tuple; +using testing::Values; + +const int kAmbisonicOrder = 3; + +// Rotation angle used in the test. +const float kAngleDegrees = 90.0f; + +// Initial, arbitrary direction of the encoded soundfield source. +const SphericalAngle kInitialSourceAngle = + SphericalAngle::FromDegrees(22.0f, 33.0f); + +// Expected direction of the sound source rotated by |kAngleDegrees| against the +// right facing axis (x in the world coordinate system). +const SphericalAngle kXrotatedSourceAngle = + SphericalAngle::FromDegrees(150.021778639249f, 51.0415207997462f); + +// Expected direction of the sound source rotated by |kAngleDegrees| against the +// upward facing axis (y in the world coordinate system). +const SphericalAngle kYrotatedSourceAngle = + SphericalAngle::FromDegrees(112.0f, 33.0f); + +// Expected direction of the sound source rotated by |kAngleDegrees| against the +// back facing axis (z in the world coordinate system). +const SphericalAngle kZrotatedSourceAngle = + SphericalAngle::FromDegrees(35.0077312770459f, -18.3108067341351f); + +// Rotation interpolation interval in terms of frames (used by the HoaRotator). +const size_t kSlerpFrameInterval = 32; + +class HoaRotatorTest : public ::testing::Test { + protected: + HoaRotatorTest() {} + + // Virtual methods from ::testing::Test + ~HoaRotatorTest() override {} + + void SetUp() override { + // Initialize the third order soundfield rotator. + hoa_rotator_ = std::unique_ptr<HoaRotator>(new HoaRotator(kAmbisonicOrder)); + } + + void TearDown() override {} + + // Test method which creates a HOA soundfield buffer, rotates it by + // |rotation_angle| against |rotation_axis|, and compares the output to the + // reference soundfield buffer (where the source is already spatialized at the + // |expected_source_angle|). + void CompareRotatedAndReferenceSoundfields( + const std::vector<float>& input_data, float rotation_angle, + const WorldPosition& rotation_axis, + const SphericalAngle& expected_source_angle) { + AudioBuffer input_buffer(1, input_data.size()); + FillAudioBuffer(input_data, 1, &input_buffer); + // Initialize a third order mono codec with the |kInitialSourceAngle|. This + // will be used to obtain the rotated soundfield. + std::unique_ptr<MonoAmbisonicCodec<>> rotated_source_mono_codec( + new MonoAmbisonicCodec<>(kAmbisonicOrder, {kInitialSourceAngle})); + // Initialize a third order mono codec with the expected source angle. This + // will be used as a reference for the rotated soundfield. + std::unique_ptr<MonoAmbisonicCodec<>> reference_source_mono_codec( + new MonoAmbisonicCodec<>(kAmbisonicOrder, {expected_source_angle})); + // Generate soundfield buffers representing sound sources at required + // angles. + AudioBuffer encoded_rotated_buffer(kNumThirdOrderAmbisonicChannels, + input_data.size()); + AudioBuffer encoded_reference_buffer(kNumThirdOrderAmbisonicChannels, + input_data.size()); + rotated_source_mono_codec->EncodeBuffer(input_buffer, + &encoded_rotated_buffer); + reference_source_mono_codec->EncodeBuffer(input_buffer, + &encoded_reference_buffer); + // Rotate the test soundfield by |rotation_angle| degrees wrt the given + // |rotation_axis|. + const WorldRotation rotation = WorldRotation( + AngleAxisf(rotation_angle * kRadiansFromDegrees, rotation_axis)); + const bool result = hoa_rotator_->Process(rotation, encoded_rotated_buffer, + &encoded_rotated_buffer); + EXPECT_TRUE(result); + // Check if the sound source in the reference and rotated buffers are in the + // same direction. + // If the buffer size is more than |kSlerpFrameInterval|, due to + // interpolation, we expect that the last |kSlerpFrameInterval| frames have + // undergone the full rotation. + // If the buffer size is less than |kSlerpFrameInterval|, because no + // interpolation is applied, the rotated soundfield should result from + // frame 0. + for (size_t channel = 0; channel < encoded_rotated_buffer.num_channels(); + ++channel) { + const int num_frames = + static_cast<int>(encoded_rotated_buffer[channel].size()); + const int interval = static_cast<int>(kSlerpFrameInterval); + const int frames_to_compare = + (num_frames % interval) ? (num_frames % interval) : interval; + const int start_frame = std::max(0, num_frames - frames_to_compare); + ASSERT_LT(start_frame, num_frames); + for (int frame = start_frame; frame < num_frames; ++frame) { + EXPECT_NEAR(encoded_rotated_buffer[channel][frame], + encoded_reference_buffer[channel][frame], kEpsilonFloat); + } + } + } + + // Higher Order Ambisonic rotator instance. + std::unique_ptr<HoaRotator> hoa_rotator_; +}; + +// Tests that no rotation is aplied if |kRotationQuantizationRad| is not +// exceeded. +TEST_F(HoaRotatorTest, RotationThresholdTest) { + const size_t kFramesPerBuffer = 16; + const std::vector<float> kInputVector( + kFramesPerBuffer * kNumThirdOrderAmbisonicChannels, 1.0f); + AudioBuffer input_buffer(kNumThirdOrderAmbisonicChannels, kFramesPerBuffer); + FillAudioBuffer(kInputVector, kNumThirdOrderAmbisonicChannels, &input_buffer); + AudioBuffer output_buffer(kNumThirdOrderAmbisonicChannels, kFramesPerBuffer); + const WorldRotation kSmallRotation = + WorldRotation(1.0f, 0.001f, 0.001f, 0.001f); + const WorldRotation kLargeRotation = WorldRotation(1.0f, 0.1f, 0.1f, 0.1f); + EXPECT_FALSE( + hoa_rotator_->Process(kSmallRotation, input_buffer, &output_buffer)); + EXPECT_TRUE( + hoa_rotator_->Process(kLargeRotation, input_buffer, &output_buffer)); +} + +typedef tuple<WorldPosition, SphericalAngle> TestParams; +class HoaAxesRotationTest : public HoaRotatorTest, + public ::testing::WithParamInterface<TestParams> {}; + +// Tests third order soundfield rotation against the x, y and z axis using long +// input buffers (>> slerp interval). +TEST_P(HoaAxesRotationTest, CompareWithExpectedAngleLongBuffer) { + const WorldPosition& rotation_axis = ::testing::get<0>(GetParam()); + const SphericalAngle& expected_angle = ::testing::get<1>(GetParam()); + const size_t kLongFramesPerBuffer = 512; + const std::vector<float> kLongInputData(kLongFramesPerBuffer, 1.0f); + CompareRotatedAndReferenceSoundfields(kLongInputData, kAngleDegrees, + rotation_axis, expected_angle); +} + +// Tests third order soundfield rotation against the x, y and z axes using short +// input buffers (< slerp interval). +TEST_P(HoaAxesRotationTest, CompareWithExpectedAngleShortBuffer) { + const WorldPosition& rotation_axis = ::testing::get<0>(GetParam()); + const SphericalAngle& expected_angle = ::testing::get<1>(GetParam()); + const size_t kShortFramesPerBuffer = kSlerpFrameInterval / 2; + const std::vector<float> kShortInputData(kShortFramesPerBuffer, 1.0f); + CompareRotatedAndReferenceSoundfields(kShortInputData, kAngleDegrees, + rotation_axis, expected_angle); +} + +// Tests third order soundfield rotation against the x, y and z axes using +// buffer sizes that are (> slerp interval) and not integer multiples of +// the slerp interval. +TEST_P(HoaAxesRotationTest, CompareWithExpectedAngleOddBufferSizes) { + const WorldPosition& rotation_axis = ::testing::get<0>(GetParam()); + const SphericalAngle& expected_angle = ::testing::get<1>(GetParam()); + const size_t frames_per_buffer = kSlerpFrameInterval + 3U; + const std::vector<float> kShortInputData(frames_per_buffer, 1.0f); + CompareRotatedAndReferenceSoundfields(kShortInputData, kAngleDegrees, + rotation_axis, expected_angle); +} + +INSTANTIATE_TEST_CASE_P( + TestParameters, HoaAxesRotationTest, + Values(TestParams({1.0f, 0.0f, 0.0f}, kXrotatedSourceAngle), + TestParams({0.0f, 1.0f, 0.0f}, kYrotatedSourceAngle), + TestParams({0.0f, 0.0f, 1.0f}, kZrotatedSourceAngle))); + +} // namespace + +} // namespace vraudio diff --git a/src/3rdparty/resonance-audio/resonance_audio/ambisonics/stereo_from_soundfield_converter.cc b/src/3rdparty/resonance-audio/resonance_audio/ambisonics/stereo_from_soundfield_converter.cc new file mode 100644 index 000000000..94462a445 --- /dev/null +++ b/src/3rdparty/resonance-audio/resonance_audio/ambisonics/stereo_from_soundfield_converter.cc @@ -0,0 +1,52 @@ +/* +Copyright 2018 Google Inc. 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 "ambisonics/stereo_from_soundfield_converter.h" + +#include "base/constants_and_types.h" +#include "dsp/gain.h" + +namespace vraudio { + +namespace { + +const float kMidSideChannelGain = 0.5f; + +} // namespace + +void StereoFromSoundfield(const AudioBuffer& soundfield_input, + AudioBuffer* stereo_output) { + DCHECK(stereo_output); + DCHECK_EQ(kNumStereoChannels, stereo_output->num_channels()); + DCHECK_EQ(soundfield_input.num_frames(), stereo_output->num_frames()); + DCHECK_GE(soundfield_input.num_channels(), kNumFirstOrderAmbisonicChannels); + const AudioBuffer::Channel& channel_audio_space_w = soundfield_input[0]; + const AudioBuffer::Channel& channel_audio_space_y = soundfield_input[1]; + AudioBuffer::Channel* left_channel_output = &(*stereo_output)[0]; + AudioBuffer::Channel* right_channel_output = &(*stereo_output)[1]; + // Left = 0.5 * (Mid + Side). + *left_channel_output = channel_audio_space_w; + *left_channel_output += channel_audio_space_y; + ConstantGain(0 /* no offset */, kMidSideChannelGain, *left_channel_output, + left_channel_output, false /* accumulate_output */); + // Right = 0.5 * (Mid - Side). + *right_channel_output = channel_audio_space_w; + *right_channel_output -= channel_audio_space_y; + ConstantGain(0 /* no offset */, kMidSideChannelGain, *right_channel_output, + right_channel_output, false /* accumulate_output */); +} + +} // namespace vraudio diff --git a/src/3rdparty/resonance-audio/resonance_audio/ambisonics/stereo_from_soundfield_converter.h b/src/3rdparty/resonance-audio/resonance_audio/ambisonics/stereo_from_soundfield_converter.h new file mode 100644 index 000000000..44d2144e7 --- /dev/null +++ b/src/3rdparty/resonance-audio/resonance_audio/ambisonics/stereo_from_soundfield_converter.h @@ -0,0 +1,35 @@ +/* +Copyright 2018 Google Inc. 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 RESONANCE_AUDIO_AMBISONICS_STEREO_FROM_SOUNDFIELD_CONVERTER_H_ +#define RESONANCE_AUDIO_AMBISONICS_STEREO_FROM_SOUNDFIELD_CONVERTER_H_ + +#include "base/audio_buffer.h" + +namespace vraudio { + +// Performs Ambisonic to stereo decode by performing Mid-Side (M-S) matrixing. +// Soundfield format assumed is Ambix (ACN/SN3D). Only first order channels +// are used for the stereo decode. +// +// @param soundfield_input Soundfield of arbitrary order. +// @param stereo_output Pointer to stereo output buffer. +void StereoFromSoundfield(const AudioBuffer& soundfield_input, + AudioBuffer* stereo_output); + +} // namespace vraudio + +#endif // RESONANCE_AUDIO_AMBISONICS_STEREO_FROM_SOUNDFIELD_CONVERTER_H_ diff --git a/src/3rdparty/resonance-audio/resonance_audio/ambisonics/stereo_from_soundfield_converter_test.cc b/src/3rdparty/resonance-audio/resonance_audio/ambisonics/stereo_from_soundfield_converter_test.cc new file mode 100644 index 000000000..aa1f5c70a --- /dev/null +++ b/src/3rdparty/resonance-audio/resonance_audio/ambisonics/stereo_from_soundfield_converter_test.cc @@ -0,0 +1,103 @@ +/* +Copyright 2018 Google Inc. 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 "ambisonics/stereo_from_soundfield_converter.h" + +#include <iterator> + +#include "third_party/googletest/googletest/include/gtest/gtest.h" +#include "base/audio_buffer.h" +#include "base/constants_and_types.h" +#include "utils/planar_interleaved_conversion.h" + +namespace vraudio { + +namespace { + +// Number of frames per buffer. +const size_t kFramesPerBuffer = 2; + +// First order ambisonic signal in the AmbiX format (W, Y, Z, X), sound source +// in the front. +const float kFirstOrderSourceFront[] = {1.0f, 0.0f, 0.0f, 1.0f, + 1.0f, 0.0f, 0.0f, 1.0f}; + +// First order ambisonic signal in the AmbiX format (W, Y, Z, X), sound source +// to the left. +const float kFirstOrderSourceLeft[] = {1.0f, 1.0f, 0.0f, 0.0f, + 1.0f, 1.0f, 0.0f, 0.0f}; + +// Tests whether the conversion to stereo from soundfield results in equal +// signal in both L/R output channels when there is a single source in front of +// the soundfield. +TEST(StereoFromSoundfieldConverterTest, StereoFromSoundfieldFrontTest) { + // Initialize the soundfield input buffer. + const std::vector<float> soundfield_data(std::begin(kFirstOrderSourceFront), + std::end(kFirstOrderSourceFront)); + AudioBuffer soundfield_buffer( + kNumFirstOrderAmbisonicChannels, + soundfield_data.size() / kNumFirstOrderAmbisonicChannels); + FillAudioBuffer(soundfield_data, kNumFirstOrderAmbisonicChannels, + &soundfield_buffer); + + // Output buffer is stereo. + AudioBuffer output(kNumStereoChannels, + soundfield_data.size() / kNumFirstOrderAmbisonicChannels); + + StereoFromSoundfield(soundfield_buffer, &output); + + // Test for near equality. + ASSERT_EQ(kNumStereoChannels, output.num_channels()); + const AudioBuffer::Channel& output_channel_left = output[0]; + const AudioBuffer::Channel& output_channel_right = output[1]; + for (size_t frame = 0; frame < kFramesPerBuffer; ++frame) { + EXPECT_NEAR(output_channel_left[frame], output_channel_right[frame], + kEpsilonFloat); + } +} + +// Tests whether the conversion to stereo from soundfield, when the sound source +// in the soundfield is to the left, results in a signal only in the L output +// channel. +TEST(StereoFromSoundfieldConverterTest, StereoFromSoundfieldLeftTest) { + // Initialize the soundfield input buffer. + const std::vector<float> soundfield_data(std::begin(kFirstOrderSourceLeft), + std::end(kFirstOrderSourceLeft)); + AudioBuffer soundfield_buffer( + kNumFirstOrderAmbisonicChannels, + soundfield_data.size() / kNumFirstOrderAmbisonicChannels); + FillAudioBuffer(soundfield_data, kNumFirstOrderAmbisonicChannels, + &soundfield_buffer); + + // Output buffer is stereo. + AudioBuffer output(kNumStereoChannels, + soundfield_data.size() / kNumFirstOrderAmbisonicChannels); + + StereoFromSoundfield(soundfield_buffer, &output); + + // Test for near equality. + ASSERT_EQ(kNumStereoChannels, output.num_channels()); + const AudioBuffer::Channel& output_channel_left = output[0]; + const AudioBuffer::Channel& output_channel_right = output[1]; + for (size_t frame = 0; frame < kFramesPerBuffer; ++frame) { + EXPECT_NEAR(output_channel_left[frame], 1.0f, kEpsilonFloat); + EXPECT_NEAR(output_channel_right[frame], 0.0f, kEpsilonFloat); + } +} + +} // namespace + +} // namespace vraudio diff --git a/src/3rdparty/resonance-audio/resonance_audio/ambisonics/utils.h b/src/3rdparty/resonance-audio/resonance_audio/ambisonics/utils.h new file mode 100644 index 000000000..8c71944d6 --- /dev/null +++ b/src/3rdparty/resonance-audio/resonance_audio/ambisonics/utils.h @@ -0,0 +1,120 @@ +/* +Copyright 2018 Google Inc. 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 RESONANCE_AUDIO_AMBISONICS_UTILS_H_ +#define RESONANCE_AUDIO_AMBISONICS_UTILS_H_ + +#include <cmath> + +#include "base/constants_and_types.h" +#include "base/logging.h" +#include "base/misc_math.h" + +namespace vraudio { + +// Returns ACN channel sequence from a degree and order of a spherical harmonic. +inline int AcnSequence(int degree, int order) { + DCHECK_GE(degree, 0); + DCHECK_LE(-degree, order); + DCHECK_LE(order, degree); + + return degree * degree + degree + order; +} + +// Returns normalization factor for Schmidt semi-normalized spherical harmonics +// used in AmbiX. +inline float Sn3dNormalization(int degree, int order) { + DCHECK_GE(degree, 0); + DCHECK_LE(-degree, order); + DCHECK_LE(order, degree); + return std::sqrt((2.0f - ((order == 0) ? 1.0f : 0.0f)) * + Factorial(degree - std::abs(order)) / + Factorial(degree + std::abs(order))); +} + +// Returns the number of spherical harmonics for a periphonic ambisonic sound +// field of |ambisonic_order| at compile-time. +// We have to use template metaprogramming because MSVC12 doesn't support +// constexpr. +template <size_t AmbisonicOrder> +struct GetNumPeriphonicComponentsStatic { + enum { value = (AmbisonicOrder + 1) * (AmbisonicOrder + 1) }; +}; + +// Returns the number of spherical harmonics for a periphonic ambisonic sound +// field of |ambisonic_order|. +inline size_t GetNumPeriphonicComponents(int ambisonic_order) { + return static_cast<size_t>((ambisonic_order + 1) * (ambisonic_order + 1)); +} + +// Returns the number of periphonic spherical harmonics (SHs) for a particular +// Ambisonic order. E.g. number of 1st, 2nd or 3rd degree SHs in a 3rd order +// sound field. +inline size_t GetNumNthOrderPeriphonicComponents(int ambisonic_order) { + if (ambisonic_order == 0) return 1; + return static_cast<size_t>(GetNumPeriphonicComponents(ambisonic_order) - + GetNumPeriphonicComponents(ambisonic_order - 1)); +} + +// Calculates the ambisonic order of a periphonic sound field with the given +// number of spherical harmonics. +inline int GetPeriphonicAmbisonicOrder(size_t num_components) { + DCHECK_GT(num_components, 0); + const int ambisonic_order = static_cast<int>(std::sqrt(num_components)) - 1; + // Detect when num_components is not square. + DCHECK_EQ((ambisonic_order + 1) * (ambisonic_order + 1), + static_cast<int>(num_components)); + return ambisonic_order; +} + +// Calculates the order of the current spherical harmonic channel as the integer +// part of a square root of the channel number. Please note, that in Ambisonics +// the terms 'order' (usually denoted as 'n') and 'degree' (usually denoted as +// 'm') are used in the opposite meaning as in more traditional maths or physics +// conventions: +// [1] C. Nachbar, F. Zotter, E. Deleflie, A. Sontacchi, "AMBIX - A SUGGESTED +// AMBISONICS FORMAT", Proc. of the 2nd Ambisonics Symposium, June 2-3 2011, +// Lexington, KY, https://goo.gl/jzt4Yy. +inline int GetPeriphonicAmbisonicOrderForChannel(size_t channel) { + return static_cast<int>(sqrtf(static_cast<float>(channel))); +} + +// Calculates the degree of the current spherical harmonic channel. Please note, +// that in Ambisonics the terms 'order' (usually denoted as 'n') and 'degree' +// (usually denoted as 'm') are used in the opposite meaning as in more +// traditional maths or physics conventions: +// [1] C. Nachbar, F. Zotter, E. Deleflie, A. Sontacchi, "AMBIX - A SUGGESTED +// AMBISONICS FORMAT", Proc. of the 2nd Ambisonics Symposium, June 2-3 2011, +// Lexington, KY, https://goo.gl/jzt4Yy. +inline int GetPeriphonicAmbisonicDegreeForChannel(size_t channel) { + const int order = GetPeriphonicAmbisonicOrderForChannel(channel); + return static_cast<int>(channel) - order * (order + 1); +} + +// Returns whether the given |num_channels| corresponds to a valid ambisonic +// order configuration. +inline bool IsValidAmbisonicOrder(size_t num_channels) { + if (num_channels == 0) { + return false; + } + // Number of channels must be a square number for valid ambisonic orders. + const size_t sqrt_num_channels = static_cast<size_t>(std::sqrt(num_channels)); + return num_channels == sqrt_num_channels * sqrt_num_channels; +} + +} // namespace vraudio + +#endif // RESONANCE_AUDIO_AMBISONICS_UTILS_H_ diff --git a/src/3rdparty/resonance-audio/resonance_audio/ambisonics/utils_test.cc b/src/3rdparty/resonance-audio/resonance_audio/ambisonics/utils_test.cc new file mode 100644 index 000000000..51dfa13d5 --- /dev/null +++ b/src/3rdparty/resonance-audio/resonance_audio/ambisonics/utils_test.cc @@ -0,0 +1,65 @@ +/* +Copyright 2018 Google Inc. 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 "ambisonics/utils.h" + +#include "third_party/googletest/googletest/include/gtest/gtest.h" +#include "base/constants_and_types.h" + +namespace vraudio { + +namespace { + +// Ambisonic Channel Numbers (ACNs) used in the tests. +const size_t kAcnChannels[] = {0, 1, 4, 8, 16, 32}; + +// Number of different ACNs to be tested. +const size_t kNumAcnChannels = sizeof(kAcnChannels) / sizeof(kAcnChannels[0]); + +// Tests that the Ambisonic order for a given channel returned by the +// GetPeriphonicAmbisonicOrderForChannel() method is correct. +TEST(UtilsTest, GetPeriphonicAmbisonicOrderForChannelTest) { + const int kExpectedOrders[] = {0, 1, 2, 2, 4, 5}; + for (size_t channel = 0; channel < kNumAcnChannels; ++channel) { + EXPECT_EQ(kExpectedOrders[channel], + GetPeriphonicAmbisonicOrderForChannel(kAcnChannels[channel])); + } +} + +// Tests that the Ambisonic degree for a given channel returned by the +// GetPeriphonicAmbisonicDegreeForChannel() method is correct. +TEST(UtilsTest, GetPeriphonicAmbisonicDegreeForChannelTest) { + const int kExpectedDegrees[] = {0, -1, -2, 2, -4, 2}; + for (size_t channel = 0; channel < kNumAcnChannels; ++channel) { + EXPECT_EQ(kExpectedDegrees[channel], + GetPeriphonicAmbisonicDegreeForChannel(kAcnChannels[channel])); + } +} + +// Tests that the ambisonic order validation returns the expected results for +// arbitrary number of channels. +TEST(UtilsTest, IsValidAmbisonicOrderTest) { + for (size_t valid_channel : {1, 4, 9, 16, 25, 36}) { + EXPECT_TRUE(IsValidAmbisonicOrder(valid_channel)); + } + for (size_t invalid_channel : {2, 3, 5, 8, 50, 99}) { + EXPECT_FALSE(IsValidAmbisonicOrder(invalid_channel)); + } +} + +} // namespace + +} // namespace vraudio |