波束形成实例

本文基于webRTC源码剖析其基于delay-sum和T-F masking混合的最大SIRsignal Interference Ratio声源分离方法.该算法在计算的复杂度,性能和鲁棒性之间做了很好的均衡.鲁棒性主要是目标估计错误的情况而提出的.该算法的核心是采用了非线性后置滤波处理.

编译测试文件

gsc@gsc-250:~/webrtc-checkout/src/webrtc/modules/audio_processing$ 目录下的"audio_processing_tests.gypi" 文件的153 行，如下：

gypi编译脚本
该文件是gypi（generated your project included）格式的文件，语法类似json和python。编译的目标是nonlinear_beamformer_test。该文件编译的源文件是

beamformer/nonlinear_beamformer_test.cc

查找该目标的所在位置：

测试文件处理流程:

该文件处理流程比较清晰，首先根据读入的wav文件，获得采样率和通道数（麦克风个数），然后获得麦克的物理位置，

从输出可以看出，输入信号采样率是16KHz，三通道，输出信号的采样率是16KHz，一个通道，将三个通道合成为了一个通道。根据这个可执行程序，就可以一步步打印相关的中间执行结果，这有利于更好的弄懂WebRTC的beamforming算法.

算法的基本思想

每个频点的幅度增益(实部和虚部)根据,

目标信号的空域协方差矩阵.
干扰信号的空域协方差矩阵.
临近频点之间波动情况.

由于语音信号是宽带信号,所以适合将语音信号从时域转换到频域,然后针对每一个频点做处理,这样处理的粒度可以非常的细.对每个频点处理的细则如下,如果说使用delay-sum方法合成的波束变换到频域后,对应频点能量较大,则说明这个频点含有目标信息的量就比较大,应该保留的比例会相对较大,同理,如果是干扰方向计算得到的频点能量较大,则削弱这个频点的力度也要大,如果没有显示给出干扰方向,可以选择第一主瓣做为干扰方向,如果用户的使用场景可以确定目标方向的,则就可以直接使用用户给定的目标方向.

对最终波束形成之后的频域信息,可以统计得到目标是否存在的概率,一个直观个理解是如果目标(声源)正确得大,则对应波束输出频点中应该含有声源的语音信息,不论是基于能量还是统计模型计算结果都可以估计目标是否存在.

该算法的后置滤波以均方误差为准则对每一个频点计算一个增益值,就实际的语音场景,其动态范围一般较宽,所以在绝大多数场景每一个频点的绝大部分能量会来自同一个声源.所以就需要根据声源信息为每一个频点选择合适的增益.

在进入代码之前,还有几个数学公式: 矩阵归一化(范数归一化):∥Rβ∥α=VHRβV,V=argmaxWWHRαW,WWH=1(11.1)(11.1)∥Rβ∥α=VHRβV,V=arg⁡maxWWHRαW,WWH=1

则归一化的R¯¯¯¯αR¯α如下:R¯¯¯¯α=1∥Rα∥αRα(11.2)(11.2)R¯α=1∥Rα∥αRα

设dndn是第n个目标信号,umum是第m干扰信号,则对每一个时频点有一个对应的目标干扰组合(d′n,u′m)(dn′,um′),则增益的表达式可以写成下式:gi=maxn′minm′gin′m′(11.3)(11.3)gi=maxn′minm′gin′m′实际使用时,为了鲁棒性会将11.3调整成下式:

gi=1−∏n′⎡⎣1−∏m′gn′m′⎤⎦(11.4)(11.4)gi=1−∏n′[1−∏m′gn′m′]

设ΦΦ是目标声源ξξ的估计值,则对于线性波束形成算法,有:Φ=WHY(11.5)(11.5)Φ=WHY后置滤波操作将上述目标估计值ξ^ξ^乘以11.4中的增益,则有:Φ^=giΦ=giWHY(11.6)(11.6)Φ^=giΦ=giWHYW由声源目标的空间位置和麦克风坐标求得,gigi通过最小均方误差求得:g∗i=argmaxg[Ai[|wHhξξ−gΦ|2]](11.7)(11.7)gi∗=arg⁡maxg[Ai[|wHhξξ−gΦ|2]]将式子11.7展开后得:g∗i=argmaxg[Ai[|wHhξξ−gΦ|2]]=argmaxg[Ai[g2|Φ|2−gwHhξξΦH−ghHξwξHΦ]](11.7)(11.7)gi∗=arg⁡maxg[Ai[|wHhξξ−gΦ|2]]=arg⁡maxg[Ai[g2|Φ|2−gwHhξξΦH−ghξHwξHΦ]]令上式等于零,可以得到最优解:g∗i=E[A_i[|wHhξ∥2∥ξ∥2]]E[Ai[∥Φ∥2]]==A_i[∥ξ∥2]wHRξwAi∥ξ∥2wHRξw+Ai[∥ϕ∥2]wHRξw=Ai[∥ξ∥2]wHRξwwHRMw(11.8)(11.8)gi∗=E[A_i[|wHhξ∥2∥ξ∥2]]E[Ai[∥Φ∥2]]==A_i[∥ξ∥2]wHRξwAi∥ξ∥2wHRξw+Ai[∥ϕ∥2]wHRξw=Ai[∥ξ∥2]wHRξwwHRMw其中RMRM是观测到信号的协方差矩阵,RM=Ai[yyH]RM=Ai[yyH],RξRξ是麦克风接收到的声源部分(目标信号加传递函数的影响)的协方差矩阵,ϕϕ指示麦克风采集到的干扰和噪声部分. 上述式子中RξRξ，RϕRϕ和RMRM是容易求得的，式11.8分母部分可以看成是波束形成器的输出功率,对其进行矩阵归一化(用到11.1和11.2),得:HwR¯¯¯¯w=A_i\[∥ϕ∥2\]∥∥Rϕ∥∥ϕ∥∥R_M∥∥_MwHR¯¯¯¯ϕw+Ai\[∥ξ∥2\]∥∥Rξ∥∥ξ∥∥R_M∥∥_MwHR¯¯¯¯ξw=Ai\[∥ϕ∥2\]∥∥Rϕ∥∥M∥∥R_M∥∥_M⋅∥∥Rϕ∥∥ϕ∥∥Rϕ∥∥M⋅∥∥∥R¯¯¯¯ϕ∥∥_w=Ai\[∥ξ∥2\]∥∥Rξ∥∥M∥∥R_M∥∥_M⋅∥∥Rξ∥∥ξ∥∥Rξ∥∥M⋅∥∥∥R¯¯¯¯ξ∥∥_w=$1−λ$⋅∥∥Rϕ∥∥ϕ∥∥Rϕ∥∥_M⋅∥∥∥R¯¯¯¯ϕ∥∥w+λ⋅∥∥Rξ∥∥ξ∥∥Rξ∥∥M⋅∥∥∥R¯¯¯¯ξ∥∥_w(11.9)(11.9)wHR¯w=A_i\[∥ϕ∥2\]∥∥Rϕ∥∥ϕ∥∥R_M∥∥_MwHR¯ϕw+Ai\[∥ξ∥2\]∥∥Rξ∥∥ξ∥∥R_M∥∥_MwHR¯ξw=Ai\[∥ϕ∥2\]∥∥Rϕ∥∥M∥∥R_M∥∥_M⋅∥∥Rϕ∥∥ϕ∥∥Rϕ∥∥M⋅∥∥∥R¯ϕ∥∥_w=Ai\[∥ξ∥2\]∥∥Rξ∥∥M∥∥R_M∥∥_M⋅∥∥Rξ∥∥ξ∥∥Rξ∥∥M⋅∥∥∥R¯ξ∥∥_w=$1−λ$⋅∥∥Rϕ∥∥ϕ∥∥Rϕ∥∥_M⋅∥∥∥R¯ϕ∥∥w+λ⋅∥∥Rξ∥∥ξ∥∥Rξ∥∥M⋅∥∥∥R¯ξ∥∥_w其中λλ是对Ai[|ξ|2]Ai[|ξ|2]的归一化:λ=Ai[|ξ|2]||Rξ||M||RM||M(11.10)(11.10)λ=Ai[|ξ|2]||Rξ||M||RM||M

从上式可以看出对波束形成的输出贡献最大的部分是期望的目标信号. 根据11.9,可以得到归一化的λλ的解:λ=∥∥R¯¯¯¯M∥∥_M−∥∥Rϕ∥∥ϕ∥∥Rϕ∥∥M∥∥R¯¯¯¯ϕ∥∥w∥∥Rξ∥∥ξ∥∥Rξ∥∥M∥∥R¯¯¯¯ξ∥∥w−∥∥Rϕ∥∥ϕ∥∥Rϕ∥∥M∥∥R¯¯¯ϕ∥∥_w=∥∥R¯¯¯¯M∥∥_w−∥∥R¯¯¯ϕ∥∥w∥∥R¯¯¯ϕ∥∥M∥∥R¯¯¯ξ∥∥w∥∥R¯¯¯ξ∥∥M−∥∥R¯¯¯ϕ∥∥w∥∥R¯¯¯ϕ∥∥_M(11.11)(11.11)λ=∥∥R¯M∥∥_M−∥∥Rϕ∥∥ϕ∥∥Rϕ∥∥M∥∥R¯ϕ∥∥w∥∥Rξ∥∥ξ∥∥Rξ∥∥M∥∥R¯ξ∥∥w−∥∥Rϕ∥∥ϕ∥∥Rϕ∥∥M∥∥R¯ϕ∥∥_w=∥∥R¯M∥∥_w−∥∥R¯ϕ∥∥w∥∥R¯ϕ∥∥M∥∥R¯ξ∥∥w∥∥R¯ξ∥∥M−∥∥R¯ϕ∥∥w∥∥R¯ϕ∥∥_M可以按能量方式去理解11.11.这样的化,式子11.8可以变为:g∗=λ⋅||RM||M||Rξ||M⋅||Rξ||w||RM||w=1−||R¯¯¯ϕ||w||R¯¯¯M||M1||R¯¯¯M||w1−||Rϕ||w||R¯¯¯ϕ||M||R¯¯¯ξ||M||R¯¯¯ξ||w(11.12)(11.12)g∗=λ⋅||RM||M||Rξ||M⋅||Rξ||w||RM||w=1−||R¯ϕ||w||R¯M||M1||R¯M||w1−||Rϕ||w||R¯ϕ||M||R¯ξ||M||R¯ξ||w

11.12的分子和分母都有||R¯¯¯¯ϕ||w||R¯¯¯¯ϕ||M||R¯ϕ||w||R¯ϕ||M,该项在波束指向方向时,值较小,但是其它方向时较大(通常小于1).为了增强鲁棒性,常使用下式:g∗=1−min{α,||R¯¯¯ϕ||w||R¯¯¯ϕ||M1||R¯¯¯M||w}1−min{α,||R¯¯¯ϕ||w||R¯¯¯ϕ||M||R¯¯¯ξ||M||R¯¯¯ξ||w}(11.13)(11.13)g∗=1−min{α,||R¯ϕ||w||R¯ϕ||M1||R¯M||w}1−min{α,||R¯ϕ||w||R¯ϕ||M||R¯ξ||M||R¯ξ||w}其中αα是选定的常数因子,通常取值为0.999. 在使用delay-sum时,由于宽带信号的空域混跌问题,权重ww在200Hz以下,以及混叠频率处需要额外处理,通常取均值做为整个频带的权重因子.

空域混叠对于时域有奈奎斯特定理:fs=1Ts≥2fmaxfs=1Ts≥2fmax,在空间上类似的有:fxa=1d≥2fxmaxfxa=1d≥2fxmax. 根据麦克距离和声源方向可以计算混跌频率：fAlias=sd⋅(1+cos(target))fAlias=sd⋅(1+cos(target))d是两麦克直线距离，

测试文件

nonlinear_beamformer_test.cc

6 NonlinearBeamformer bf(array_geometry, array_geometry.size());
7 bf.Initialize(kChunkSizeMs, in_file.sample_rate());
8 while (in_file.ReadSamples(interleaved.size(),

&
interleaved[0]) == interleaved.size() ) {
9 FloatS16ToFloat(
&
interleaved[0], interleaved.size(), 
&
interleaved[0]);
10 Deinterleave(
&
interleaved[0], buf.num_frames(),
buf.num_channels(), buf.channels());
count++;
11 bf.AnalyzeChunk(buf);
12 bf.PostFilter(
&
buf);
13 Interleave(buf.channels(), buf.num_frames(),
buf.num_channels(), 
&
interleaved[0]);
14 FloatToFloatS16(
&
interleaved[0], interleaved.size(), 
&
interleaved[0]);
out_file.WriteSamples(
&
interleaved[0], interleaved.size());
}

6是波束形成的实例化,需要把麦克风的坐标和数量放在参数里,7行的第一个参数是数据长度(8ms~32ms)之间都可以,第二个参数是采样率.9,10行是对读入的数据进行了重新排列.11行主要是用于计算11.13中的g∗g∗,12行主要是对对其进行T-F相乘. 第6行创建的对象是整个算法的而核心:

// Enhances sound sources coming directly in front of a uniform linear array
// and suppresses sound sources coming from all other directions. Operates on
// multichannel signals and produces single-channel output.
//
// The implemented nonlinear postfilter algorithm taken from "A Robust Nonlinear
// Beamforming Postprocessor" by Bastiaan Kleijn.
class NonlinearBeamformer : public LappedTransform::Callback {
public:
//半波数弧度值,是角分辨率的度量,webrtc默认设置成了20度.
static const float kHalfBeamWidthRadians;
//这里设置了目标的初始方向.
explicit NonlinearBeamformer(
const std::vector
<
Point
>
&
 array_geometry,
size_t num_postfilter_channels = 1u,
SphericalPointf target_direction =
SphericalPointf(static_cast
<
float
>
(M_PI) / 2.f, 0.f, 1.f));
~NonlinearBeamformer() override;

// Sample rate corresponds to the lower band.
// Needs to be called before the NonlinearBeamformer can be used.
virtual void Initialize(int chunk_size_ms, int sample_rate_hz);

//根据初始化的参数,将时域信号分到频域,计算11.3的每个频点的g值
virtual void AnalyzeChunk(const ChannelBuffer
<
float
>
&
 data);

//根据AnalyzeChunk计算到的g值,做后滤波,然后在变换到时域
virtual void PostFilter(ChannelBuffer
<
float
>
* data);
//设置目标的坐标的API(球坐标:依次是方位角,俯仰角,半径)
virtual void AimAt(const SphericalPointf
&
 target_direction);
//根据波束形成后的情况,判断目标是否存在
virtual bool IsInBeam(const SphericalPointf
&
 spherical_point);
virtual bool is_target_present();

protected:
//时频处理的核心函数
void ProcessAudioBlock(const complex
<
float
>
* const* input,
size_t num_input_channels,
size_t num_freq_bins,
size_t num_output_channels,
complex
<
float
>
* const* output) override;

private:
FRIEND_TEST_ALL_PREFIXES(NonlinearBeamformerTest,
InterfAnglesTakeAmbiguityIntoAccount);

typedef Matrix
<
float
>
 MatrixF;
typedef ComplexMatrix
<
float
>
 ComplexMatrixF;
typedef complex
<
float
>
 complex_f;
//低频部分是200!~400,这段不是每个频点使用一个g,而是使用这个频段的平均值.
void InitLowFrequencyCorrectionRanges();
//根据阵列信息计算得到混跌频率,从混跌频率到二分之一采样率范围为高频部分,处理方式和200Hz~400Hz类似
void InitHighFrequencyCorrectionRanges();
//如果没有指定干扰方向的信息,则这里可以按照第一旁瓣接收到的信号做为干扰信号,也可以是目标方向6dB带宽之外的信号做为干扰信号.
void InitInterfAngles();
//目标方向的信号
void InitDelaySumMasks();
//计算目标信号的协方差矩阵
void InitTargetCovMats();
//计算散射协方差矩阵
void InitDiffuseCovMats();
//计算干扰方向的协方差矩阵
void InitInterfCovMats();
//归一化方法11.1
void NormalizeCovMats();
//g值计算核心函数
float CalculatePostfilterMask(const ComplexMatrixF
&
 interf_cov_mat,
float rpsiw,
float ratio_rxiw_rxim,
float rmxi_r);

// Prevents the postfilter masks from degenerating too quickly (a cause of
// musical noise).
void ApplyMaskTimeSmoothing();
void ApplyMaskFrequencySmoothing();

// The postfilter masks are unreliable at low frequencies. Calculates a better
// mask by averaging mid-low frequency values.
void ApplyLowFrequencyCorrection();

// Postfilter masks are also unreliable at high frequencies. Average mid-high
// frequency masks to calculate a single mask per block which can be applied
// in the time-domain. Further, we average these block-masks over a chunk,
// resulting in one postfilter mask per audio chunk. This allows us to skip
// both transforming and blocking the high-frequency signal.
void ApplyHighFrequencyCorrection();

// Compute the means needed for the above frequency correction.
float MaskRangeMean(size_t start_bin, size_t end_bin);

// Applies post-filter mask to |input| and store in |output|.
void ApplyPostFilter(const complex_f* input, complex_f* output);

void EstimateTargetPresence();

static const size_t kFftSize = 256;
static const size_t kNumFreqBins = kFftSize / 2 + 1;

// Deals with the fft transform and blocking.
size_t chunk_length_;
std::unique_ptr
<
LappedTransform
>
 process_transform_;
std::unique_ptr
<
PostFilterTransform
>
 postfilter_transform_;
float window_[kFftSize];

// Parameters exposed to the user.
const size_t num_input_channels_;
const size_t num_postfilter_channels_;
int sample_rate_hz_;

const std::vector
<
Point
>
 array_geometry_;
// The normal direction of the array if it has one and it is in the xy-plane.
const rtc::Optional
<
Point
>
 array_normal_;

// Minimum spacing between microphone pairs.
const float min_mic_spacing_;

// Calculated based on user-input and constants in the .cc file.
size_t low_mean_start_bin_;
size_t low_mean_end_bin_;
size_t high_mean_start_bin_;
size_t high_mean_end_bin_;

// Quickly varying mask updated every block.
float new_mask_[kNumFreqBins];
// Time smoothed mask.
float time_smooth_mask_[kNumFreqBins];
// Time and frequency smoothed mask.
float final_mask_[kNumFreqBins];

float target_angle_radians_;
// Angles of the interferer scenarios.
std::vector
<
float
>
 interf_angles_radians_;
// The angle between the target and the interferer scenarios.
const float away_radians_;

// Array of length |kNumFreqBins|, Matrix of size |1| x |num_channels_|.
ComplexMatrixF delay_sum_masks_[kNumFreqBins];

// Arrays of length |kNumFreqBins|, Matrix of size |num_input_channels_| x
// |num_input_channels_|.
ComplexMatrixF target_cov_mats_[kNumFreqBins];
ComplexMatrixF uniform_cov_mat_[kNumFreqBins];
// Array of length |kNumFreqBins|, Matrix of size |num_input_channels_| x
// |num_input_channels_|. The vector has a size equal to the number of
// interferer scenarios.
std::vector
<
std::unique_ptr
<
ComplexMatrixF
>
>
 interf_cov_mats_[kNumFreqBins];

// Of length |kNumFreqBins|.
float wave_numbers_[kNumFreqBins];

// Preallocated for ProcessAudioBlock()
// Of length |kNumFreqBins|.
float rxiws_[kNumFreqBins];
// The vector has a size equal to the number of interferer scenarios.
std::vector
<
float
>
 rpsiws_[kNumFreqBins];

// The microphone normalization factor.
ComplexMatrixF eig_m_;

// For processing the high-frequency input signal.
float high_pass_postfilter_mask_;
float old_high_pass_mask_;

// True when the target signal is present.
bool is_target_present_;
// Number of blocks after which the data is considered interference if the
// mask does not pass |kMaskSignalThreshold|.
size_t hold_target_blocks_;
// Number of blocks since the last mask that passed |kMaskSignalThreshold|.
size_t interference_blocks_count_;
};

权重计算函数

void
 NonlinearBeamformer::ProcessAudioBlock(
const
 complex_f* 
const
* input,

size_t
 num_input_channels,

size_t
 num_freq_bins,

size_t
 num_output_channels,
complex_f* 
const
* output) {
RTC_CHECK_EQ(kNumFreqBins, num_freq_bins);
RTC_CHECK_EQ(num_input_channels_, num_input_channels);
RTC_CHECK_EQ(
0
, num_output_channels);

#
if
 0
for
(
size_t
 ch =
0
; ch 
<
num_input_channels_; ++ch ){

for
(
size_t
 f_ix=
0
; f_ix 
<
 kNumFreqBins; ++ f_ix)

std
::
cout
<
<
"fft out["
<
<
ch 
<
<
"]["
<
<
 f_ix 
<
<
"]"
<
<
 input[ch][f_ix] 
<
<
std
::
endl
;
}


#
endif
std
::
cout
<
<
"NonlinearBeamformer::ProcessAudioBlock"
<
<
std
::
endl
;

// Calculating the post-filter masks. Note that we need two for each
// frequency bin to account for the positive and negative interferer
// angle.
//对每个频点进行的操作
for
 (
size_t
 i = low_mean_start_bin_; i 
<
= high_mean_end_bin_; ++i) {

//计算输入信号的归一化矩阵,实际上是噪声+目标+干扰的结果

eig_m_.CopyFromColumn(input, i, num_input_channels_);

float
 eig_m_norm_factor = 
std
::
sqrt
(SumSquares(eig_m_));

if
 (eig_m_norm_factor != 
0.f
) {
eig_m_.Scale(
1.f
 / eig_m_norm_factor);
}

//11.1中的范数归一化,由此计算得到(目标+噪声+干扰)的能量.
float
 rxim = Norm(target_cov_mats_[i], eig_m_);

float
 ratio_rxiw_rxim = 
0.f
;

if
 (rxim 
>
0.f
) {
ratio_rxiw_rxim = rxiws_[i] / rxim;
}

//计算得到目标的能量

complex_f rmw = 
abs
(ConjugateDotProduct(delay_sum_masks_[i], eig_m_));
rmw *= rmw;

float
 rmw_r = rmw.real();

//根据干扰和目标计算得到权重值

new_mask_[i] = CalculatePostfilterMask(*interf_cov_mats_[i][
0
],
rpsiws_[i][
0
],
ratio_rxiw_rxim,
rmw_r);

//对于面阵,或者有明确干扰源的情况下,还要进行计算
for
 (
size_t
 j = 
1
; j 
<
 interf_angles_radians_.size(); ++j) {

float
 tmp_mask = CalculatePostfilterMask(*interf_cov_mats_[i][j],
rpsiws_[i][j],
ratio_rxiw_rxim,
rmw_r);

//保守的权值跟新策略
if
 (tmp_mask 
<
 new_mask_[i]) {
new_mask_[i] = tmp_mask;
}
}
}

//在时域上对权值进行平滑

ApplyMaskTimeSmoothing();

//低频情况下使用平均代替单个频点

ApplyLowFrequencyCorrection();

//混叠频率的处理和低频情况下一致

ApplyHighFrequencyCorrection();

//频域上的平滑处理

ApplyMaskFrequencySmoothing();
}

权重相乘操作

void
 PostFilterTransform::ProcessAudioBlock(
const
complex
<
float
>
* 
const
* input,

size_t
 num_input_channels,

size_t
 num_freq_bins,

size_t
 num_output_channels,

complex
<
float
>
* 
const
* output) {
RTC_DCHECK_EQ(num_freq_bins_, num_freq_bins);
RTC_DCHECK_EQ(num_input_channels, num_output_channels);

printf
(
"PostFilterTransform::ProcessAudioBlock, \n"
);

for
 (
size_t
 ch = 
0
; ch 
<
 num_input_channels; ++ch) {

for
 (
size_t
 f_ix = 
0
; f_ix 
<
 num_freq_bins_;                  ++f_ix) {
     output[ch][f_ix] =kCompensationGain *    final_mask_[f_ix] * input[ch][f_ix];

   }
}

本章小节

本章主要根据开源的波束算法进行剖析额算法的原理,算法的公式,和相关核心公式的代码实现,要了解算法的精髓还是需要源码一行行的看的.该算法既有delay-sum经典的算法又有T-Fmasking的思想在里面了,弥补了各自的不足,dleay-sum是线性处理,而T-F masking是非线性处理,实际上在有些场景下线性是适用的,在有些场景下非线性是适用的(T-F masking只是一种特例),所以经常需要考虑将多个算法恰到好处的综合成一个,这个难度还是有的,一种新的处理方法正能实现这种综合功能而具有旺盛的生命力,就是基于深度学习的方法,下一章将浅析这一算法.

波束形成进阶

Last chapter a beamforming example based on T-F masking method already analysised.

上一章分析了t-F masking掩码的实例,并给了代码下载的github地址.

That method is based on the geometray of microphone array and DOA(usually base on SRP-PHATH algorithm) estimate. Several questions need to be step furtuer.

那个方法基于麦克风阵列坐标以及先验的声源风向(通常基于SRP-PHAT计算获得).这种方法是大多数基于波束形成的声源分离方法.有些地方还是值得推敲的.

In masse production, if the array corrdinate error is about ± 1 m m ± 1 m m , what's the influnce.
对于大型生产场景,如果麦克风坐标误差对波束形成的影响如何?
DOA不准的时候,波束形成的结果是怎样的?(在其它算法中被称为鲁棒性,如MVDR中会使用对角阵的方式增加鲁棒性)

如何将上面两个对波束形成的因数去除掉?基于盲源分离(PCA,主成份分析)是一种方法,但是这里要分析的依然还是基于波束形成的,但是波束指向的方向向量不再依赖于DOA结果,而是根据观测到的信号自主决定,这就是chime3挑战赛中的方法CGMM.如果细致跟过前面章节中VAD检测的的ML方法,那么这里的方法理解起来也不难.

高斯混合模型

一维高斯分布

p(x)=∑i=1KϕiN(x|μi,σi),∑i=1Kϕi=1(12.1)(12.1)p(x)=∑i=1KϕiN(x|μi,σi),∑i=1Kϕi=1N(x|μi,σi)=1σi2π√exp(−(x−μi)22σ2i)(12.2)(12.2)N(x|μi,σi)=1σi2πexp(−(x−μi)22σi2)

多维高斯分布

CGMM方法

CGMM(complex gaussian mixture model)是基于混合高斯模型的方法,和之前的波束形成方法相比,差异在于根据高斯混合模型求解指向目标方向的向量.这个方法的处理框图如下:

CGMM波束形成架构

该波束形成方法包括波束形成,指向向量以及时频掩码三个模块,这三个模型的最终输出就是增强后的信号.

阵列接收模型假设

设k∈1,2,⋯,Kk∈1,2,⋯,K是声源目标的索引,m∈1,2,⋯,Mm∈1,2,⋯,M是麦克风的索引.在时域每个麦克风接收到的信号可以写做如下形式:ym(t)=∑k∑τh(k)m(τ)s(k)(t−τ)+nm(t)(12.5)(12.5)ym(t)=∑k∑τhm(k)(τ)s(k)(t−τ)+nm(t)上式中,s(k)(t)s(k)(t)和nm(t)nm(t)分别表示第k个声源和以及第m个麦克采集到的噪声,h(k)i(τ)hi(k)(τ)是第k个声源对第m个麦克的脉冲相应函数.对12.5使用STFT变换后可得:ym(f,t)=∑kh(k)m(f)s(k)(f,t)+nm(f,t)(12.6)(12.6)ym(f,t)=∑khm(k)(f)s(k)(f,t)+nm(f,t)向量表示法如下:⎡⎣⎢⎢⎢⎢y1y2⋮ym⎤⎦⎥⎥⎥⎥=⎡⎣⎢⎢⎢⎢⎢h11y12⋮y1m⋯⋯⋯hm1hm2ymm⎤⎦⎥⎥⎥⎥⎥⎡⎣⎢⎢⎢⎢s1s2⋮sm⎤⎦⎥⎥⎥⎥+⎡⎣⎢⎢⎢⎢n1n2⋮nm⎤⎦⎥⎥⎥⎥=y(f,t)=∑kr(k)(f)s(k)(f,t)+n(f,t)(12.7)(12.7)[y1y2⋮ym]=[h11⋯h1my21⋯h2m⋮ym1⋯ymm][s1s2⋮sm]+[n1n2⋮nm]=y(f,t)=∑kr(k)(f)s(k)(f,t)+n(f,t)

波束形成

波束形成方法采用了MVDR,该方法对频域信号乘以权重获得到原始信号的估计:s^(k)f,t=w(k)fHyf,t(12.8)(12.8)s^f,t(k)=wf(k)Hyf,t

s^(k)f,ts^f,t(k)是第k个目标信号的估计值,HH是共轭转置(频域信号), 最小均方误差准则的限制条件是w(k)Hfr(k)f=1wf(k)Hrf(k)=1,则可以推导出MVDR波束形成的权重计算公式是:w(k)f=R(y)−1frfr(k)HfR(y)−1fr(k)f(12.9)(12.9)wf(k)=Rf(y)−1rfrf(k)HRf(y)−1rf(k)其中R(y)fRf(y)是根据观测信号得到的协方差矩阵:R(y)f=1T∑tyf,tyHf,t(12.10)(12.10)Rf(y)=1T∑tyf,tyf,tH当然也可以使用其它波束形成算法,如多通道维纳滤波.

指向向量

这个波束形成算法的关键在于指向向量的计算,该方法直接根据麦克风阵列的协方差矩阵求得,使用协方差矩阵的主特征向量做为波束指向向量,这非常类似盲源分离的PCA思想,协方差矩阵可以使用T-F masking获得. 设λ(k+n)f,tλf,t(k+n)和λ(n)f,tλf,t(n)是时频点(f,t)(f,t)第k阶带噪信号和只有噪声信号的掩码(概率). 则带噪语音信号和只有噪声信号的协方差矩阵可以表示如下:R(k+n)f=1Σtλ(n)f,t∑tλ(k+n)f,tyf,tyHf,t(12.11)(12.11)Rf(k+n)=1Σtλf,t(n)∑tλf,t(k+n)yf,tyf,tH

R(n)f=1Σtλ(n)f,t∑tλ(n)f,tyf,tyHf,t(12.12)(12.12)Rf(n)=1Σtλf,t(n)∑tλf,t(n)yf,tyf,tH

假设噪声是不相关的,则估计的纯净语音信号由下式计算:R(k)f=R(k+n)f−R(n)f(12.13)(12.13)Rf(k)=Rf(k+n)−Rf(n)

则根据12.13可知,第k个声源的指向向量可以通过对R(k)fRf(k)使用特征向量分解法求得,即最大特征值对应的特征向量做为指向向量.

时频掩码估计

复高斯混合模型

时频掩码基于复高斯混合模型,

图12.2 高斯混合模型图

基于语音信号在时频域的稀疏性,多通道的语音生成模型可以使用K+1K+1个的高斯分布来表示,每一个对应于一个带噪语音类别.(这里的分类和VAD里的不太一致的地方,在于,其将带噪语音信号类别的总数被聚类到K+1个),混合高斯模型的参数估计使用EM算法. 式16.716.7可以写为:yf,t=r(v)fs(v)f,t,df,t=v(12.14)(12.14)yf,t=rf(v)sf,t(v),df,t=v其中,df,tdf,t是时频点(f,t)(f,t)所属类别索引,vv可以取k+n(k=1⋯K)或者nk+n(k=1⋯K)或者n,分别表示第kk个带噪语音信号和噪声信号.s(k+n)f,tsf,t(k+n)表示的是第k个类表示的带噪语音信号.s(n)f,tsf,t(n)是时频点(f,t)(f,t)的噪声信号表示.这里假设带噪语音信号和噪声信号的指向向量是时不变的. 基于上述观测模型,可以定义生成模型,假设s(v)f,tsf,t(v)复合复高斯分布:s(v)f,t∼Nc(0,ϕ(v)f,t)(12.15)(12.15)sf,t(v)∼Nc(0,ϕf,t(v))其中,ϕ(v)f,tϕf,t(v)是方差,Nc(x|μ,Σ)=1πΣexp−|x−μ|2ΣNc(x|μ,Σ)=1πΣexp−|x−μ|2Σ.从式12.15和式12.14可得多通道高斯分布如下:yf,t|df,t=v∼Nc(0,ϕ(v)f,tRf(v))(12.16)(12.16)yf,t|df,t=v∼Nc(0,ϕf,t(v)Rf(v)),R(v)fRf(v)是r(v)fr(v)Hfrf(v)rf(v)H,式子12.16中的Nc(x|μ,Σ)=1|πΣ|exp(−(x−μ)HΣ−1(x−μ))Nc(x|μ,Σ)=1|πΣ|exp(−(x−μ)HΣ−1(x−μ)),使用混合权重可以将式子12.16转为下式:yf,t∼∑vα(v)fNc(0,ϕ(v)f,tR(v)f),Σvα(v)f=1(12.17)(12.17)yf,t∼∑vαf(v)Nc(0,ϕf,t(v)Rf(v)),Σvαf(v)=1对r(v)frf(v)的估计通过估计R(v)fRf(v)获得.

基于EM算法的参数估计

CGMM参数,如α(v)f,ϕ(v)f,t,R(v)fαf(v),ϕf,t(v),Rf(v)使用最大似然法估计(maximum likelihood,ML).对于上一小节提出的复高斯混合模型,其对数目标函数是:logp(y|Θ)=∑f,tlog∑vα(v)fNc(yf,t|0,ϕ(v)f,tR(v)f),∑vα(v)f=1(12.18)(12.18)log⁡p(y|Θ)=∑f,tlog⁡∑vαf(v)Nc(yf,t|0,ϕf,t(v)Rf(v)),∑vαf(v)=1这里ΘΘ是CGMM的参数集.在每一次EM计算中的Q函数定义如下:Q(Θ|Θ′)=E[logp(y|Θ,d)]d=∑f,t∑vλ(v)f,tlogα(v)fNc(yf,t|0,ϕ(v)f,tR(v)f)(12.19)(12.19)Q(Θ|Θ′)=E[log⁡p(y|Θ,d)]d=∑f,t∑vλf,t(v)log⁡αf(v)Nc(yf,t|0,ϕf,t(v)Rf(v))E[⋅]xE[⋅]x表示的是xx的期望,Θ′Θ′是前一次的参数估计集合,λ(v)f,t(满足Σvλ(v)f,t=1)λf,t(v)(满足Σvλf,t(v)=1)是df,tdf,t属于vv的后验概率.在E步骤中后验概率的计算公式是:λ(v)f,t←α(v)fp(yf,t|df,t=v,Θ′)∑vα(v)fp(yf,t|df,t=v,Θ′)(12.19)(12.19)λf,t(v)←αf(v)p(yf,t|df,t=v,Θ′)∑vαf(v)p(yf,t|df,t=v,Θ′)其中p(yf,t|df,t=v,Θ′)=Nc(yf,t|0,ϕ(v)f,tR(v)f)p(yf,t|df,t=v,Θ′)=Nc(yf,t|0,ϕf,t(v)Rf(v)).在M步骤,CGMM参数估计如下:ϕ(v)f,t←1Mtr(yf,tyHf,tR(v)−1f)(12.20)(12.20)ϕf,t(v)←1Mtr(yf,tyf,tHRf(v)−1)R(v)f←1∑tλ(v)f,t∑vλ(v)f,t1ϕ(v)f,tyf,tyHf,t(12.21)(12.21)Rf(v)←1∑tλf,t(v)∑vλf,t(v)1ϕf,t(v)yf,tyf,tHα(v)f←1T∑tλ(v)f,t(12.22)(12.22)αf(v)←1T∑tλf,t(v)在EM步骤收敛后,就可以认为T-F掩码就是对应的λnf,tλf,tn值. 对于单声源目标情况,E-M算法的初始设置,可以使用单位阵做为噪声的空域互相关矩阵,使用观测的信号所得的互相关矩阵做为带噪信号的互相关矩阵.多声源情况,可以对观测到的信号的互相关矩阵加满秩随机矩阵做为带噪信号的互相关矩阵. 对于matlab代码以及所实现的效果,可以参考CGMM 波束形成git网址

排列问题

EM算法收敛后,需要对每个频点的带噪语音和噪声进行聚类,这被称为排列问题.这在盲源分离时也会遇到.

基于CGMM方法的在线波束形成

假设语音信号根据,以batch形式获得,每一个batch 按照时间序列依次被采集到,batch的索引l∈1,⋯,Ll∈1,⋯,L,BlBl是第l个batch所包含的数据,对于第ll个batch,λ(v)f,tλf,t(v)和ϕ(v)f,tϕf,t(v)根据式12.19和12.20求得,其中的R(v)f,lRf,l(v)则根据第l−1l−1的R(v)f,l−1Rf,l−1(v)迭代获得,将式12.21修改为下式:R(v)f,l←Λ(v)f,l−1Λ(v)f,l−1+∑t∈Blλ(v)f,tR(v)f,l−1+1Λ(v)f,l−1+∑t∈Blλ(v)f,t∑t∈Blλ(v)f,t1ϕ(v)f,tyf,tyHf,t(12.23)(12.23)Rf,l(v)←Λf,l−1(v)Λf,l−1(v)+∑t∈Blλf,t(v)Rf,l−1(v)+1Λf,l−1(v)+∑t∈Blλf,t(v)∑t∈Blλf,t(v)1ϕf,t(v)yf,tyf,tH这里Λ(v)f,lΛf,l(v)是多有帧λ(v)f,tλf,t(v)的和,其递归跟新方法是:Λ(v)f,l←Λ(v)f,l−1+∑t∈Blλ(v)f,t(12.24)(12.24)Λf,l(v)←Λf,l−1(v)+∑t∈Blλf,t(v)使用序列方法估计的时频掩码λ(v)f,tλf,t(v),则实时波束形成如下:R(v)f,l←Λ(v)f,l−1Λ(v)f,l−1+∑t∈Blλ(v)f,tR(v)f,l−1+1Λ(v)f,l−1+∑t∈Blλ(v)f,t∑t∈Blλ(v)f,tyf,tyHf,t(12.25)(12.25)Rf,l(v)←Λf,l−1(v)Λf,l−1(v)+∑t∈Blλf,t(v)Rf,l−1(v)+1Λf,l−1(v)+∑t∈Blλf,t(v)∑t∈Blλf,t(v)yf,tyf,tH波束指向向量rr依然是根据RR计算获得.然后BlBl内的带噪语音使用MVDR方法获得增强后的信号.带噪语音的协方差矩阵的跟新方式如下:R(y)f,t←Tl−1TlR(y)f,l−1+1Tl∑t∈Blyf,tyHf,t(12.26)(12.26)Rf,t(y)←Tl−1TlRf,l−1(y)+1Tl∑t∈Blyf,tyf,tH其中TlTl是到第ll个batch为止观测到的帧总数.在实时计算算法流程是:

初始化R(y)f,0←0,R(k+n)f,0←0,R(n)f,0←0,Λ(k+n)f,0←0,Λ(n)f,0←0Rf,0(y)←0,Rf,0(k+n)←0,Rf,0(n)←0,Λf,0(k+n)←0,Λf,0(n)←0使用前文提到的方法或者训练数据训练得到:Rk+nf,0Rf,0k+n和R(n)f,0Rf,0(n)for l=1 to L in all k and all f do

根据12.19和12.20对 t ∈ B l t ∈ B l 计算 λ ( k + n ) f , t , λ ( n ) f , t , ϕ ( k + n ) f , t , ϕ ( n ) f , t λ f , t ( k + n ) , λ f , t ( n ) , ϕ f , t ( k + n ) , ϕ f , t ( n )

根据12.23跟新 R ( k + n ) f , l R f , l ( k + n ) 和 R ( n ) f , l R f , l ( n )

根据12.25跟新 R ( k + n ) f , l R f , l ( k + n ) 和 R ( n ) f , l R f , l ( n )

根据12.13跟新 R ( k ) f , l R f , l ( k )

提取 R ( y ) f , l R f , l ( y ) 的特征向量做为波束形成的指向向量 r ( k ) f , l r f , l ( k )

根据12.26更新 R ( y ) f , l R f , l ( y )

根据12.9更新 w ( k ) f , l w f , l ( k )

根据12.8做波束形成,得到语音的估计 s ^ ( k ) f , t , t ∈ B l s ^ f , t ( k ) , t ∈ B l

根据12.24更新 Λ ( k + n ) f , l Λ f , l ( k + n ) 和 Λ ( n ) f , l Λ f , l ( n ) end for

波束成形（Beam-forming）