A spherical harmonic-domain spatial audio signal enhancement method based on minimum variance distortionless response
Huawei Zhang Jihui Zhang Huiyuan Sun Prasanga Samarasinghe








A spherical harmonic-domain spatial audio signal enhancement method based on minimum variance distortionless response 
This work is sponsored by the Australian Research Council (ARC) Discovery Projects funding scheme with project number DE230101567.

    Huawei Zhang^1, Jihui (Aimee) Zhang^2, 1, Huiyuan (June) Sun^1, Prasanga Samarasinghe^1
^1Audio & Acoustic Signal Processing Group, The Australian National University, Canberra, Australia 

^2Institute of Sound and Vibration Research, University of Southampton, Southampton, U.K.

    September 5, 2024
============================================================================================================================================================================================================================================================================================



§ ABSTRACT 
Spatial audio signal enhancement aims to reduce interfering source contributions while preserving the desired sound field with its spatial cues intact. 

Existing methods generally rely on impractical assumptions (e.g. no reverberation or accurate estimations of impractical information) or have limited applicability.   

This paper presents a spherical harmonic (SH)-domain minimum variance distortionless response (MVDR)-based spatial signal enhancer using Relative Harmonic Coefficients (ReHCs) to extract clean SH coefficients from noisy ones in reverberant environments.




A simulation study shows the proposed method achieves lower estimation error, higher speech-distortion-ratio (SDR), and comparable noise reduction (NR) within the sweet area in a reverberant environment, compared to a beamforming-and-projection method as the baseline.




Spatial audio signal enhancement, relative harmonic coefficient, relative transfer function, minimum variance distortionless response (MVDR), beamforming, spherical harmonic.




§ INTRODUCTION

















With the rapid increase of spatial audio applications in the last decade, there has been a growing demand for spatial audio signal enhancement (also known as ambisonic-to-ambisonic separation) <cit.>, which separates a desired sound field from interference source contributions while preserving the desired spatial cues.



This technology acts as a fundamental step in many audio signal processing applications including sound field recording <cit.>, sound field reproduction <cit.>, and spatial active noise control <cit.>.










We categorize existing methods for spatial audio enhancement into three kinds: (i) Beamforming-and-projection methods <cit.>, (ii) Multi-channel Wiener Filter methods <cit.>, and (iii) Learning based methods <cit.>.

Beamforming-and-project methods often include at least two stages: a beamforming stage to capture the desired source signal or the desired sound pressure at a point, and a projection stage to reconstruct the desired sound field using the point-to-region acoustic transfer functions (ATFs).

While these solutions have been achieved over space via multi-point processing <cit.> as well as spherical harmonic (SH) domain processing <cit.>, they make one or more impractical assumptions such as no reverberation or the accurate estimation of ATFs. 

Multi-channel Wiener Filter methods achieve enhancement by preserving the estimated power spectral density (PSD) matrix of the desired SH coefficients and reducing the estimated PSD matrix of interference SH coefficients.

However, the accurate estimation of the former PSD matrix is difficult.

Learning-based methods up to today, have mainly focused on separating speech sources in the SH domain and required a suitable dataset for pre-training.

In short, existing methods are generally either with impractical assumptions or applied in limited conditions.






















Relative Harmonic Coefficients (ReHCs) <cit.>, also named as SH-domain relative transfer functions <cit.>, are denoted as the ratio between SH coefficients and the specific SH coefficient at the 0-th order and the 0-th mode.

These ReHCs have been used in beamformer design <cit.> and spatial audio signal enhancement <cit.>.

These ReHCs can be estimated based on corresponding SH-domain PSD matrices obtained from the received microphone signals <cit.>.

However, the potential of ReHCs has not been fully studied.

In this paper, we develop an SH-domain minimum variance distortionless response (MVDR) spatial audio signal enhancement method with multiple beamformers based on ReHCs.


Motivated by a binaural MVDR beamforming method <cit.> and a multi-output MVDR beamforming method <cit.>, we extend the conventional SH-domain MVDR beamforming method <cit.> to be a multi-output variant.

This proposed method extracts the SH coefficients due to the desired sound source from the mixed SH domain recording.

Different from the conventional single-output MVDR beamforming method (losing the spatial information) <cit.> and the beamforming-and-projection methods (requiring a projection stage to reconstruct spatial information), this proposed method can directly extract the spatial information of the desired sound field.


Moreover, the proposed beamforming method is developed with a more practical foreknowledge of ReHCs, which is used to construct a set of spatial constraints to preserve desired SH coefficients. 








































§ PROBLEM FORMULATION








As shown in Fig. <ref>, consider a scenario with a desired sound source, an interference sound source, with a spherical sweet area with the radius r_s and a Q-microphone spherical array of the radius r_a placed concentrically at the origin O in a reverberated room.

The two sound sources are placed at fixed positions outside the spherical area, causing a desired sound field and an interference sound field within the sweet area, respectively.

We assume the desired sound source is uncorrelated with the interference sound source, and the interference sound source signal is stable.













In this paper, we focus on far-field scenarios.

























    
    

    








Each sound field pressure within the sweet area is denoted as x(t, k, r, θ, ϕ) in the time-frequency (TF) domain, where t is the time frame, k = 2 π f / c is the wavenumber, f is the frequency, c is the speed of sound, r, θ, and ϕ are the radius, the elevation angle, and the azimuth angle, respectively. 

The sound field pressure x(t, k, r, θ, ϕ) can be decomposed as


    x(t, k, r, θ, ϕ) = d(t, k, r, θ, ϕ) + v(t, k, r, θ, ϕ),
 

where d(t, k, r, θ, ϕ) and v(t, k, r, θ, ϕ) are the sound pressure of the desired field and the interference field, respectively.

The received microphone signal at the q-th microphone is denoted as

    x_q(t, k)    = x(t, k, r_q, θ_q, ϕ_q) + u(t, k) 
       = d_q(t, k) + v_q(t, k) + u(t, k),

where d_q(t, k), v_q(t, k), u(t, k) are the desired microphone signal, the interference microphone signal, and the random sensor noise, respectively, r_q, θ_q, and ϕ_q are the radius, the elevation angle, and the azimuth angle, respectively, at the q-th microphone. 








    
    

    









The microphone signal x_q(t, k) can be decomposed into SH coefficients to represent the mixed field as <cit.>


    x_q (t, k) = ∑_n=0^N_k∑_m=-n^n x̃_nm(t, k) j_n(kr_q) Y_n m(θ_q, ϕ_q),
 

where x̃_nm(t, k) is the mixed SH coefficient representing the mixed field at the n-th order and the m-th mode, j_n(·) is the Spherical Bessel function of the first kind, Y_nm(·) is the SH function, and N_k is the corresponding maximum order for the wavenumber k defined as <cit.> 


    N_k = ceil(kr_a) = ceil(2 π f/cr_a),
 

where ceil(  ) is the ceiling function.










For brevity, the time frame t and the wave number k will be omitted in the rest of this paper.




Obtained L = (N+1)^2 mixed SH coefficients 𝐱̃ = [x̃_00, x̃_1-1, ⋯, x̃_nm, ⋯, x̃_N N]^T can be divided as


    𝐱̃   = SHT(𝐱) = SHT(𝐝 + 𝐯 + 𝐮) =  𝐝̃ + 𝐯̃ + 𝐮̃,
 



where the SHT is the Spherical Harmonic Transform based on Eq. (<ref>) <cit.>, 𝐱 = [x_1, ⋯, x_q, ⋯, x_Q]^T are received microphone signals at the array, 𝐝, 𝐯 and 𝐮 are received signals from the desired field, the interference field, and sensor noises, respectively, 𝐝̃, 𝐯̃ and 𝐮̃ are SH coefficients of the desired field, the interference field, and sensor noises, respectively.




The task of this work is to extract desired SH coefficients 𝐝̃ from mixed SH coefficients 𝐱̃ obtained from the mixture recording 𝐱.



































  







§ PROPOSED MULTI-OUTPUT MVDR METHOD






In this section, we use ReHCs to develop an SH-domain multi-output MVDR beamforming method to estimate desired SH coefficients.




































The proposed method is a multi-output extension of a SH-domain MVDR beamforming method, consisting of L beamformers 𝐰̃ = [𝐰̃_00^T, 𝐰̃_1-1^T, ⋯, 𝐰̃_nm^T, ⋯, 𝐰̃_NN^T]^T as a L^2 × 1 vector <cit.>.

The estimated desired SH coefficient d̂̃̂_nm at the n-th order and the m-th mode is obtained with the corresponding beamformer 𝐰̃_nm as


    d̂̃̂_nm = 𝐰̃_nm^H 𝐱̃.
 





The beamformer 𝐰̃_nm is required to preserve the desired SH coefficient and reduce SH coefficients of interference and sensor noises.

As it is difficult to obtain accurate ATFs, similar to <cit.>, we use ReHCs to develop a cost function with a spatial constraint as


    minimize_𝐰̃_nm     𝐰̃_nm^H 𝐑̃_𝐯+𝐮𝐰̃_nm
    subject  to     𝐰̃_nm^H 𝐡̃ = h̃_nm,


where 𝐑̃_𝐯+𝐮 = 𝐑̃_𝐯 + 𝐑̃_𝐮 is the summaration of the SH-domain PSD matrices of interference 𝐑̃_𝐯 = 𝔼 [𝐯̃𝐯̃^H] and that of sensor noises 𝐑̃_𝐮 = 𝔼 [𝐮̃𝐮̃^H], 𝔼[  ] is math expectation, 𝐡̃ = [h̃_00, h̃_1-1, ⋯, h̃_nm, ⋯, h̃_NN]^T are ReHCs for the desired field, and h̃_nm is the corresponding ReHC at the n-th order and the m-th mode.

These ReHCs 𝐡̃ can be obtained as <cit.>


    𝐡̃
        = 𝐝̃/d̃_00
        =𝐑̃_𝐝𝐞̃_1/𝐞_1^H 𝐑̃_𝐝𝐞̃_1≈𝐑̂̃̂_𝐝+𝐮𝐞̃_1/𝐞̃_1^H 𝐑̂̃̂_𝐝+𝐮𝐞̃_1,


where d̃_00 is the desired SH coefficient at the 0-th order and the 0-th mode, 𝐑̃_𝐝 = 𝔼 [𝐝̃𝐝̃^H] is the SH-domain PSD matrix of desired SH coefficients, 𝐞̃_1 = [1   0_1 × (L-1)]^T, 𝐑̂̃̂_𝐝+𝐮 is the estimated SH-domain PSD matrix obtained from received microphone signals when only the desired sound source is active.














Combining cost functions as Eq. (<ref>) in all orders and all modes, a new cost function is derived <cit.> to design 𝐰̃ as 


    minimize_𝐰̃     𝐰̃^H ℛ̃_𝐯+𝐮𝐰̃
    subject  to     𝐂̃^H 𝐰̃ = 𝐛̃,


where 𝐛̃ = conj(𝐡̃) is a L × 1 column vector, conj(  ) is conjugate operation, ℛ̃_𝐯+𝐮 and 𝐂̃ are two large matrices constructed by 𝐑̃_𝐯+𝐮 and 𝐡̃, respectively, as


    ℛ̃_𝐯+𝐮 = 
      [ 𝐑̃_𝐯+𝐮                     ;        𝐑̃_𝐯+𝐮              ;                    ⋯       ;                      𝐑̃_𝐯+𝐮 ]_L^2 × L^2,



    𝐂̃ = 
      [ 𝐡̃         ;    𝐡̃      ;        ⋯   ;          𝐡̃ ]_L^2 × L.


The optimal solution of Eq. (<ref>) can be obtained as <cit.>


    𝐰̃ = ℛ̃_𝐯+𝐮^-1𝐂̃ [𝐂̃^H ℛ̃_𝐯+𝐮^-1𝐂̃]^-1𝐛̃.


Each beamformer 𝐰̃_nm can be obtained from 𝐰̃ directly.

After that, we can obtain an estimation 𝐝̂̃̂ of desired SH coefficients for all orders and all modes with Eq. (<ref>).


The estimated sound field pressures with spatial cues over the sweet area can be further obtained based on the Inverse Spherical Harmonic Transform (ISHT) <cit.> from the estimated desired SH coefficients 𝐝̂̃̂.

However, the spatial cues of the residual interference sound field may not be preserved, similar to <cit.>.






















































§ SIMULATION AND RESULTS


In this section, we evaluate the performance of the proposed method compared with a beamforming-and-projection method <cit.> as the baseline in a reverberant room.

The baseline is implemented with TF-domain MVDR beamforming <cit.>, requiring the summation of the TF-domain PSD matrix of interference and that of sensor noises.











 §.§ Setup









A speech signal <cit.> is used as the desired source at (4.60, 4.05,
 1.70) m.

A recorded washer-dryer noise signal <cit.> is used as the interference source at (1.60, 1.05, 1.20) m. 



We set a 32-mic open spherical array at (1.60, 4.05, 1.70) m, with the radius r_a = 0.042 m.

The array has the same size and microphone positions as the em32 Eigenmike spherical microphone array <cit.>. 

White Gaussian noises are added to the microphone array as sensor noises in terms of 35 dB signal-sensor noise-ratio (SSNR), where the ‘signal’ in SSNR refers to the desired source signal.

The room size is 5×6×4 m and the reverberation time is T_60 = 0.2 s.

The sampling frequency is 16000 Hz.

The radius of the sweet area is r_s = r_a = 0.042 m.



The speed of sound is 343 m/s.



All room impulse responses (RIRs) are simulated with a toolbox using the image source method <cit.>.















We first apply the Short-Time Fourier Transform with the frame size as 16384 and 75% overlap to the time-domain received microphone signals. 

Then we apply the SHT with the corresponding maximal order N_k to each frequency bin of obtained TF-domain microphone signals.



We simulate accurate PSD matrices of interference signals and sensor noises as prior knowledge in the TF domain and the SH domain, facilitating the baseline method and the proposed method, respectively.

The PSD matrix of interference signals is calculated based on the simulated interference source signal power and corresponding ATFs from simulated RIRs.

The PSD matrix of sensor noises is calculated by the averaged estimated PSD matrix of simulated sensor noises.





For the baseline method <cit.>, ATFs for the desired sound source are estimated during a mixture recording, requiring the accurate direction of arrival (DoA) as prior knowledge.

In <cit.>, frequency smooth operation has been performed with 9 frequency bins around each target frequency bin.

The following simulation results in Sec. <ref> also show such frequency smooth is robust to different acoustic environments. 


As the output of the baseline method is the estimated desired microphone signals, for a fair comparison with the proposed method, the SHT is applied to the output of the baseline method to estimate desired SH coefficients.





We can then reconstruct the corresponding estimated desired field by applying the ISHT from estimated desired SH coefficients of both the proposed method and the baseline method.


We evaluate the performance of both the proposed method and the baseline method among the common telephone bandwidth (from 300 Hz to 3400 Hz) <cit.> due to: (1) The key information of speech signals exist within the bandwidth; (2) The Spherical Bessel function of the first kind j_n(kr) will achieve values close to 0 among some higher frequency bins, resulting in the Bessel Zero problem <cit.>.

Therefore, at most 3-order SH coefficients are required based on Eq. (<ref>) in the following simulations.

The Bessel Zero problem can be relieved by using the spherical-shell microphone array <cit.>.



























 §.§ Performance Analysis Based on Sound Field Estimation Error








We compare the estimated desired sound fields and evaluate the corresponding estimation errors within the sweet area when the signal-noise-ratio (SNR) is 0 dB, where the ‘signal’
in SNR refers to the desired source signal and the ‘noise’ in SNR refers to the interference source signal.

The normalized square error ϵ(t, k) as the estimation error of the estimated desired sound field (sound pressures) within the sweet area is defined as:


    ϵ(t, k) = 10 log_10‖𝐝_s(t, k) - 𝐝̂_s(t, k) ‖^2/‖𝐝_s(t, k) ‖^2,


where 𝐝_s(t, k) and 𝐝̂_s(t, k) are the true desired sound field pressures and the estimated desired sound field pressures at each observation point respectively, and ‖  ‖^2 denotes the 2-norm operation.




















Fig. <ref> compares the estimated desired field by the proposed method, the estimated desired field by the baseline method, and the true desired field at one time frame at 1500 Hz on the x-y plane. 





Each sound field with corresponding estimation error defined in Eq. (<ref>) is presented for 441 observation points evenly distributed over the x-y plane.



By comparing Fig. <ref>(b), Fig. <ref>(c), and Fig. <ref>(e), we observe that both the proposed method and the baseline method reconstructed a desired field similar to the true desired field within the sweet area.

Fig. <ref>(d) and Fig. <ref>(f) show that the proposed method achieved a lower estimation error of less than about -15 dB, compared with the baseline method. 






 §.§ Performance Analysis over Frequency





In Fig. <ref>, we further evaluate the average of the sound field estimation error, the speech-distortion-ratio (SDR), and the noise reduction (NR) within the sweet area across 15 time frames against different signal frequencies. 

Here, SDR and NR are defined as <cit.>


    SDR(t, k) = 10 log_10‖𝐝_s(t, k) ‖^2/‖𝐝^res_s(t, k) - 𝐝_s(t, k) ‖^2,




    NR(t, k) = 10 log_10‖𝐯_s(t, k) ‖^2/‖𝐯^res_s(t, k) + 𝐮^res_s(t, k) ‖^2,


where 𝐯_s(t, k) is true interference field pressures, 𝐝^res_s(t, k) and 𝐯^res_s(t, k) are residual desired sound field pressures and residual interference sound field pressures after processing, respectively, and 𝐮^res_s(t, k) is representing the influence of residual sensor noises, at these observation points within the spherical sweet area.


The evaluation is based on corresponding sound pressures at 107 observation points evenly distributed over the spherical sweet area.

As shown in Fig. <ref>, the proposed method with estimated ReHCs achieved lower estimation error, higher SDR, and comparable NR (more than 25 dB NR is effective enough) than the baseline method ranging from 300 Hz to 3400 Hz.











In detail, the proposed with estimated ReHCs achieved a lower than -15 dB estimation error and a higher than 15 dB SDR at the majority of chosen frequency bins, which outperforms the baseline.




Fig. <ref> also shows the influence of the accuracy of ReHCs. 

As a comparison, we play a segment of 0 dB white Gaussian noise instead of the original desired source signal at the same source position, and then use the same method in Sec. <ref> to obtain another set of ReHCs as accurate ReHCs.

These ReHCs are more accurate as white Gaussian noise is more stable than speech signals at different frequency bins.

In Fig. <ref>, we can find the proposed method with accurate ReHCs outperformed that with estimated ReHCs in terms of estimation error and SDR, and achieved comparable NR performance, which implies the performance of the proposed method is influenced by the accuracy of estimation ReHCs.


















 §.§ Performance Analysis over Varying Reverberation Times and SNR Levels




Here, we evaluate the performance of both the proposed method and the baseline method for varying reverberation times T_60 and varying SNRs, averaging within the spherical sweet area at the bandwidth for 15 time frames.

The evaluation is based on corresponding sound pressures at 107 observation points as Sec. <ref>.

The estimated ReHCs are used in the proposed method as Sec. <ref>.


Firstly, we evaluate the performance of these two methods over varying reverberation times T_60 as shown in TABLE <ref>.

Here, we vary the T_60 from 0 s to 0.4 s, while the SNR is set to be 0 dB.


With the reverberation time T_60 increased, both methods achieved higher estimation error, lower SDR, and lower NR.

In addition, the proposed method always achieved lower estimation error, higher SDR, and comparable NR than the baseline method for varying T_60 (more than 25 dB NR is effective enough). 


Secondly, we evaluate the performance of these two methods over varying SNR levels as shown in TABLE <ref>.

Here, we vary the SNR level from 5 dB to -5 dB, and remain the T_60 as 0.2 s.

TABLE <ref> shows that the proposed method maintained stable performance for varying SNR levels in terms of estimation error (about -19 dB) and SDR (about 22 dB), always outperforming the baseline method.

By contrast, with the SNR level decreased, both estimation error and SDR achieved by the baseline method declined.

In addition, both the proposed method and the baseline method achieved comparable and stable NR (more than about 28 dB) for varying SNR levels.







§ CONCLUSION


In this paper, we propose a spherical harmonic (SH)-domain minimum variance distortionless response (MVDR) method to estimate the desired sound field from the mixture recording at a spherical microphone array.

We use a cost function with a set of spatial constraints to extract desired SH coefficients and suppress SH coefficients of interference and sensor noises.

The field due to the desired sound source within the sweet area can be hence reconstructed.

Simulation results show that the proposed method can outperform the baseline methods within the sweet area in a reverberant room.


In the future, we plan to further examine the proposed method in real-world scenarios.









IEEEtran




