语音增强经典模型论文精读与实现

发表于 2026-06-25 更新于 2026-07-17

FullSubNet

核心思路：将一个纯全频带模型和一个纯子频带模型依次连接起来，并利用实际的联合训练将这两种模型的优点结合起来。

全频带模型：输入全频带带噪语音频谱，输出全频带预测纯净语音的模型。全频带模型可以捕获全局上下文谱和长距离交叉频带依赖，但缺乏信号平稳性建模和关注局部谱模式的能力
子频带模型：模型的输入由一个频率和多个上下文频率组成。输出是对应频率的纯净语音。所有频率都是独立处理。由于噪声较于语音更加平稳，因此其可以通过对局部频带的平稳性建模区分噪声和语音；但对于信噪比极低的子带，效果不佳，因为其无法利用全频带信息。

通俗理解，子频带模型负责学习不同频段输入的分布，以得到不同频段的增强规则；全频带模型则负责学习全频带上频率的分布情况，以补充子频带模型没有利用全频带信息的缺陷

输入

首先需要强调一点，在FullSubNet中，两个模型的建模方向均为时间轴（核心！！！），也就是输出序列的长度为T

全频带模型：输入直接是带噪语音的stft频谱
子频带模型：全频带输出的格式为（B，1，F，T），原始FullSubNet会先将其reshape为（BxF，1，1，T），第一个1代表通道数，在此处没有意义；第二个1代表选择的当前频段，即全频带输出会将每个频段分离，每个频段信号都视作单个样本；此外，FullSubNet还会从原始带噪语音中选择相同频段和周围频段（（BxF，1，N，T），N为全频带当前样本频段的上下N个相邻频段），最后与全频带输出进行拼接，也就是子频带模型的输入为（BxF，1，1+N，T）

模型结构

全频带模型和子频带模型的结构完全一致，但由于子频带模型样本更小，因此其LSTM中隐藏层也相对小了一点；同时，由于子带模型输出的就是当前T时刻，对应频段的mask，因此也没有布置激活函数

原文使用的归一化是直接对每个样本的所有维度求均值再除以该均值

训练目标

模型的输出，也就是子频带模型输出为（B，F，2，T），分别代表实部虚部的mask，这样就可以同时修正幅度和相位：

$\hat M(t,f)=\hat M_r(t,f)+j\hat M_i(t,f)$

而对于带噪语音，其可以表示为：

$Y=S+N → S=MY$

因此理想mask可以表示为：

$M(t,f)= \frac {S}{Y}$

然后直接拟合mask即可

总结

其实FullSubNet的核心就是在全频带模型后接了一个子频带模型，该子频带模型用共享的参数处理所有的子频带特征，其他就没啥了

DCCRN

CRN

在介绍DCCRN前，我们需要先知道它的核心框架，CRN：Convolutional Recurrent Network

在语音领域，CRN同时结合了CNN局部能力强、RNN能学习长程依赖的优点，提出了一个经典的Encoder-Decoder结构：

Encoder

由若干层Conv2D和BN、ELU组成

每层Conv2D的参数通常为：kernel=(5/3,2),stride=(2,1)对应输入形状（B，C，F，T），即Encoder会通过卷积逐步降采样频率信息，而时间长度则保持不变（但时间轴上也会做卷积）

对于离线增强任务，使用普通卷积，而对于在线增强推理任务，由于未来信息不可知，因此会使用掩码卷积（causal conv），由于卷积核是有宽度的，因此其原理会比Transformer中的causal mask使用的三角矩阵更复杂一丢丢：

假设现在数据为x0 x1 x2，卷积核大小为3，步长为1

那对x0的传统卷积窗口则为：pad x0 x1，即pad默认在两边

而在实时场景中t0时刻我们只有x0，没有x1，因此我们的causal conv的窗口即为：pad pad x0即全部在左边pad

再举个例子就是，我们如果对x2进行causal conv，窗口即为x0 x1 x2，在数值上等同于对x1进行普通conv，但是他们代表的意义是不同的（对应的时间步不同）

源码中的实现如下：

1
2
3

def forward(self,inputs):
    if self.padding[1] != 0 and self.causal:
        inputs = F.pad(inputs,[self.kernel_size-1, 0,0,0])

这里顺便说一下F.pad的用法，对于第二个参数，也就是这里的[self.kernel_size-1, 0,0,0]，其代表的是：

在最后一维的左边添加self.kernel_size-1个pad，右边不添加pad，在倒数第二维左右都不添加pad。最重要的就是其顺序是从最后一维开始指定，该参数最小为2维，最多不设上限，也就是这里也可以就是[self.kernel_size-1, 0]

LSTM

Encoder最终会输出（B，C’，F‘，T），CRN中将其reshape为（B，T，C‘ x F’）作为输入

# [2, 256, 4, 200] = [2, 1024, 200] => [2, 200, 1024]
lstm_in = e_5.reshape(batch_size, n_channels * n_f_bins, n_frame_size).permute(0, 2, 1)
lstm_out, _ = self.lstm_layer(lstm_in)  # [2, 200, 1024]
lstm_out = lstm_out.permute(0, 2, 1).reshape(batch_size, n_channels, n_f_bins, n_frame_size)  # [2, 256, 4, 200]

Decoder

Decoder使用ConvTranspose2d将LSTM输出还原为原始输入形状，对于反卷积，其原理就是插零+普通卷积进行上采样。为了防止上采样丢失高频细节，Decoder与Encoder之间还存在对应的skip connection以保留高频细节：

d_1 = self.tran_conv_block_1(torch.cat((lstm_out, e_5), 1))
d_2 = self.tran_conv_block_2(torch.cat((d_1, e_4), 1))
d_3 = self.tran_conv_block_3(torch.cat((d_2, e_3), 1))
d_4 = self.tran_conv_block_4(torch.cat((d_3, e_2), 1))
d_5 = self.tran_conv_block_5(torch.cat((d_4, e_1), 1))

最终，Decoder输出增强的Mask

CRN和DCCRN都是参考的U-Net的skip connection结构，即使用cat，而不是ResNet系的叠加，这样可以更好保留encoder和decoder各自的特征，但计算消耗会变大

CRN以及传统AE任务的缺陷

CRN预测的仅为幅值谱，在最终增强中使用的是noisy的相位，在低信噪比下性能有限。而对于预测复数谱的模型来说，输出的实部和虚部仅仅是输出的两个channel，模型并没有被显式告知这两个channel满足复数乘法规律，并且这些模型内部计算也全部都是实数表示，也就是说仅仅通过拟合来使模型最终输出倾向于一个实部值一个虚部值，并没有还原复数本身的结构

Complex Conv

因此，DCCRN直接提出了一个可以直接接收复数输入的卷积模块：

对于输入：

$X_r+jX_i$

卷积核对应的也为：

$W_r+jW_i$

按照复数乘法规律，则卷积输出为：（DCCRN核心公式）

$Y_r=W_r∗X_r−W_i∗X_i$ $Y_i=W_r∗X_i+W_i∗X_r$

也就是说，一层Complex Conv就进行了4次卷积，更符合STFT的物理意义

其具体代码就是将Wr和Wi分别定义为两个卷积层：

if self.complex_axis == 0:#已将实部虚部cat在同一维
    real = self.real_conv(inputs)
    imag = self.imag_conv(inputs)
    real2real,imag2real = torch.chunk(real,2, self.complex_axis)
    real2imag,imag2imag = torch.chunk(imag,2, self.complex_axis)
    
else:
    if isinstance(inputs, torch.Tensor):
        real,imag = torch.chunk(inputs, 2, self.complex_axis)

    real2real = self.real_conv(real,)
    imag2imag = self.imag_conv(imag,)

    real2imag = self.imag_conv(real)
    imag2real = self.real_conv(imag)

real = real2real - imag2imag
imag = real2imag + imag2real
out = torch.cat([real, imag], self.complex_axis)

对应的，Batchnorm也加入了协方差来使复数整体归一化，激活函数也改为了复数输入形式

complex LSTM

格式与complex conv基本一致，也是将两层LSTM结构化复数格式，不再赘述

输出

DCCRN最终输出的是CRM（复数比例mask），即：

$(\hat S_r+ j\hat S_i)=(M_r+jM_i)(Y_r+jY_i)$

损失函数为复数谱的MSE和SI-SNR

LibriSpeech实战

模型参数配置

# 创建模型 (默认配置)
model = DCCRN(
  rnn_layers=2,        # LSTM 层数
  rnn_units=128,       # LSTM 隐藏单元数
  win_len=400,         # STFT 窗长，16K对应25ms
  win_inc=100,         # STFT 帧移
  fft_len=512,         # FFT 点数
  win_type='hann',     # 窗函数类型
  masking_mode='E',    # 掩蔽模式: 'E', 'C', 'R'
  use_clstm=False,     # 是否使用复数LSTM
  use_cbn=False,       # 是否使用复数BatchNorm
  kernel_size=5,       # 卷积核大小
  kernel_num=[16,32,64,128,256,256]  # 各层通道数
)

对于masking_mode，对应网络最终输出的mask格式：

E：估计幅度掩蔽 + 相位修正 (tanh 限制幅度)
C：复数乘法掩蔽 (直接估计复数 mask)
R：实数掩蔽 (实部和虚部分别乘 mask)

use_clstm和use_cbn则是选择是否要在LSTM和BN中也使用复数计算模型，当使用复数LSTM时，官方代码的卷积核对应也增大为了kernel_num=[32, 64, 128, 256, 256, 256]

在论文中，如下参数是性能最好的配置：

model = DCCRN(
    rnn_layers=2,
    rnn_units=128,
    win_len=400,
    win_inc=100,
    fft_len=512,
    win_type='hann',
    masking_mode='E',
    use_clstm=True,
    use_cbn=True,
    kernel_size=5,
    kernel_num=[32, 64, 128, 256, 256, 256],
)

前向传播

def forward(self, inputs, lens=None):
    specs = self.stft(inputs)
    real = specs[:,:self.fft_len//2+1]
    imag = specs[:,self.fft_len//2+1:]
    spec_mags = torch.sqrt(real**2+imag**2+1e-8)
    spec_mags = spec_mags
    spec_phase = torch.atan2(imag, real)
    spec_phase = spec_phase
    cspecs = torch.stack([real,imag],1)
    cspecs = cspecs[:,:,1:]
    '''
    means = torch.mean(cspecs, [1,2,3], keepdim=True)
    std = torch.std(cspecs, [1,2,3], keepdim=True )
    normed_cspecs = (cspecs-means)/(std+1e-8)
    out = normed_cspecs
    ''' 

    out = cspecs
    encoder_out = []
    
    for idx, layer in enumerate(self.encoder):
        out = layer(out)
    #    print('encoder', out.size())
        encoder_out.append(out)
    
    batch_size, channels, dims, lengths = out.size()
    out = out.permute(3, 0, 1, 2)
    if self.use_clstm:
        r_rnn_in = out[:,:,:channels//2]
        i_rnn_in = out[:,:,channels//2:]
        r_rnn_in = torch.reshape(r_rnn_in, [lengths, batch_size, channels//2*dims])
        i_rnn_in = torch.reshape(i_rnn_in, [lengths, batch_size, channels//2*dims])
    
        r_rnn_in, i_rnn_in = self.enhance([r_rnn_in, i_rnn_in])

        r_rnn_in = torch.reshape(r_rnn_in, [lengths, batch_size, channels//2, dims])
        i_rnn_in = torch.reshape(i_rnn_in, [lengths, batch_size, channels//2, dims]) 
        out = torch.cat([r_rnn_in, i_rnn_in],2)
    
    else:
        # to [L, B, C, D]
        out = torch.reshape(out, [lengths, batch_size, channels*dims])
        out, _ = self.enhance(out)
        out = self.tranform(out)
        out = torch.reshape(out, [lengths, batch_size, channels, dims])
   
    out = out.permute(1, 2, 3, 0)
    
    for idx in range(len(self.decoder)):
        out = complex_cat([out,encoder_out[-1 - idx]],1)
        out = self.decoder[idx](out)
        out = out[...,1:]
    #    print('decoder', out.size())
    mask_real = out[:,0]
    mask_imag = out[:,1] 
    mask_real = F.pad(mask_real, [0,0,1,0])
    mask_imag = F.pad(mask_imag, [0,0,1,0])
    
    if self.masking_mode == 'E' :
        mask_mags = (mask_real**2+mask_imag**2)**0.5
        real_phase = mask_real/(mask_mags+1e-8)
        imag_phase = mask_imag/(mask_mags+1e-8)
        mask_phase = torch.atan2(
                        imag_phase,
                        real_phase
                    ) 

        #mask_mags = torch.clamp_(mask_mags,0,100) 
        mask_mags = torch.tanh(mask_mags)
        est_mags = mask_mags*spec_mags 
        est_phase = spec_phase + mask_phase
        real = est_mags*torch.cos(est_phase)
        imag = est_mags*torch.sin(est_phase) 
    elif self.masking_mode == 'C':
        real,imag = real*mask_real-imag*mask_imag, real*mask_imag+imag*mask_real
    elif self.masking_mode == 'R':
        real, imag = real*mask_real, imag*mask_imag 
    
    out_spec = torch.cat([real, imag], 1) 
    out_wav = self.istft(out_spec)
     
    out_wav = torch.squeeze(out_wav, 1)
    #out_wav = torch.tanh(out_wav)
    out_wav = torch.clamp_(out_wav,-1,1)
    return out_spec,  out_wav

ConvSTFT

在DCCRN中，由于2020年torch还没有推出支持自动求导的STFT，因此其使用Conv1d函数来实现了一个可微的STFT，需要注意的是，此处使用Conv1d正常来说仅是为了实现STFT，因此其权重是固定为傅里叶基的

由代码可得知，模型的输入直接就是（b，采样点），因此直接输入音频即可

另外还有一点要说明一下，在语音领域一个常用的STFT配置是：

fs=16000
win_length=400
hop=100
n_fft=512

此处窗长小于了傅里叶变换次数，但对于离散傅里叶变换（可见信号处理基础算法 | 小董的BLOG），我们要求时间采样点（也就是STFT的窗长）要等于傅里叶变换次数，在几乎所有库函数中，会自动补0，使样本长度=512。像这样配置的好处是可以用更少的时域数据获得更平滑的频率轴，同时也可以增大输入维度，但这实际上不能真正增加频域分辨率

加载数据集

代码比较繁长，就不放了，注意以下几点即可

噪声信噪比范围问题

语音使用librispeech，噪声使用选煤厂噪声，当时测下来噪声平均声压级为100 db左右，如果假设大声说话声压级为80 db左右，因此最大信噪比会来到<-20 db，实际训练信噪比取[-20,10]

音频长度截取

DCCRN和大部分离线训练的AE模型很多都是取4s为一个样本。对于每个音频，<4s的进行重复拼接，>4s的则随机裁剪4s的长度

归一化问题

如果是直接加载的wav，那其默认数值范围就是[-1,1]。由于信号混合时都是直接通过信噪比控制数值，而通常做法就是根据语音的能量（RMS）控制噪声的增益来实现信噪比的控制，所以要先将语音控制在[-1,1]，通常有两种做法来实现归一化：

由于wav加载的数值都是[-1,1]（如果不是就直接归一化），因此直接对带噪信号进行峰值归一化：

noisy = clean + noise
# 峰值归一化: 混合后可能超出 [-1, 1], 同步缩放 noisy 和 clean
peak = np.max(np.abs(noisy)) + 1e-10
if peak > 1.0:
    noisy = noisy / peak
    clean = clean / peak

即仅对较强信号进行一个峰值缩放，注意语音信号也要进行相同缩放，这样才能使语音本身数值一致

直接控制clean和noise本身的能量RMS，先对齐到同一较低级别，这样即使相加也不会超过[-1,1]的范围，且可以确保整个数据集中所有样本的能量分布符合真实的声学统计规律，而不是被强制拉伸到统一的峰值。

另外，为了防止模型输出音频也超过范围，通常会在模型最后添加一个Tanh激活函数（直接生成波形）或者直接映射到sigmoid（mask方法）

损失函数

源论文是直接使用的仅SI-SDR，也可以使用MSE等常规loss

def si_snr(s1, s2, eps=1e-8):
    #s1 = remove_dc(s1)
    #s2 = remove_dc(s2)
    s1_s2_norm = l2_norm(s1, s2)
    s2_s2_norm = l2_norm(s2, s2)
    s_target =  s1_s2_norm/(s2_s2_norm+eps)*s2
    e_nosie = s1 - s_target #分母
    target_norm = l2_norm(s_target, s_target)#分子
    noise_norm = l2_norm(e_nosie, e_nosie)
    snr = 10*torch.log10((target_norm)/(noise_norm+eps)+eps)
    return torch.mean(snr)

SI-SDR值越大越好，因此作为loss的话要取一个负号

离线训练转在线流式推理

流式缓存维护

LSTM缓存维护

由于在流式推理中，是把很多chunk拼在一起，而网络是一个chunk一个chunk的处理数据，为了保证处理当前chunk时保留之前chunk的记忆，需要对模型中LSTM进行缓存维护：

def forward(self, inputs, states=None):
    if isinstance(inputs,list):
        real, imag = inputs
    elif isinstance(inputs, torch.Tensor):
        real, imag = torch.chunk(inputs,-1)

    # states format:
    #   ( (h_r, c_r), (h_i, c_i) )  for real_lstm
    #   ( (h_r, c_r), (h_i, c_i) )  for imag_lstm
    if states is not None:
        (h0_rl_r, c0_rl_r), (h0_rl_i, c0_rl_i) = states[0]  # real_lstm
        (h0_il_r, c0_il_r), (h0_il_i, c0_il_i) = states[1]  # imag_lstm
    else:
        h0_rl_r = c0_rl_r = h0_rl_i = c0_rl_i = None
        h0_il_r = c0_il_r = h0_il_i = c0_il_i = None

    # real_lstm processes 'real' and 'imag' independently
    if h0_rl_r is not None:
        r2r_out, (h_rl_r, c_rl_r) = self.real_lstm(real, (h0_rl_r, c0_rl_r))
        i2r_out, (h_rl_i, c_rl_i) = self.real_lstm(imag, (h0_rl_i, c0_rl_i))
    else:
        r2r_out, (h_rl_r, c_rl_r) = self.real_lstm(real)
        i2r_out, (h_rl_i, c_rl_i) = self.real_lstm(imag)

    # imag_lstm processes 'real' and 'imag' independently
    if h0_il_r is not None:
        r2i_out, (h_il_r, c_il_r) = self.imag_lstm(real, (h0_il_r, c0_il_r))
        i2i_out, (h_il_i, c_il_i) = self.imag_lstm(imag, (h0_il_i, c0_il_i))
    else:
        r2i_out, (h_il_r, c_il_r) = self.imag_lstm(real)
        i2i_out, (h_il_i, c_il_i) = self.imag_lstm(imag)

    real_out = r2r_out - i2i_out
    imag_out = i2r_out + r2i_out
    if self.projection_dim is not None:
        real_out = self.r_trans(real_out)
        imag_out = self.i_trans(imag_out)

    new_states = (
        ((h_rl_r, c_rl_r), (h_rl_i, c_rl_i)),   # real_lstm
        ((h_il_r, c_il_r), (h_il_i, c_il_i)),   # imag_lstm
    )
    return [real_out, imag_out], new_states

源码中，共有2个LSTM，这2个LSTM会分别都进行一次实数和虚数计算，而这些计算都应该分开维护，因此最终我们需要维护的缓存共有4套，即对于同一LSTM，实数计算和虚数计算的隐藏层和细胞层都应该独立维护

nn.LSTM的输入可以直接在第二个参数输入隐藏层和细胞层参数（即上一chunk的最后一个时间步的隐藏状态），其输出固定为[output,(h, c)]

卷积层维护

当模型逐个chunk处理时，为了实现causal conv，每个chunk在左侧都会进行1个pad（卷积核大小为2的情况），所以理论上每个chunk在第每一层encoder，decoder中都有一个开头0pad，那我们在后续chunk的建模中，需不需要把这个0pad替换为上一chunk的最后一个时间步呢？具体答案会在下一节流式数据处理中讲到

流式数据处理

对于流式数据的处理，有两种方法，它们的核心就是是否使用额外overlap

首先，无论是哪种方法，chunk之间都一定会存在overlap，它们是为了补充以下内容：

STFT的padding。除非STFT使用center=false，所有的STFT会在左、右侧添加padding（DCCRN自定义的STFT只填充左边），防止对首尾加窗时超出界限。对于第一个之后的chunk，我们希望这个padding由上一段音频替代，这样可以获得更稳定的首帧数据，而我们如果将上一chunk尾部的部分数据拼接到当前chunk，就可以实现。如果overlap的长度=这个stft的padding，那么每个chunk之间的STFT频谱边界就是完全连续的

其实最好的做法就是使用center=false，这样甚至都不用overlap了，但这样在训练中就会麻烦一点点。
causal Conv的padding。前面我们提到过causal conv的实现也是通过padding，我们同样不希望除第一个chunk之外的chunk出现padding，这样每一个chunk在卷积层都会被视作是“第一个chunk”。因此对于每一个causal conv，我们都需要上一个chunk的最后一个时间帧的特征。

举个例子，DCRNN由6层encoder和decoder组成，每层中都有一个causal conv，其会在输入特征第一个时间帧之前做一帧padding，因此对于第一层causal conv，我们需要上一个chunk对应的第一层输入的最后一帧来代替padding，卷积会将这层padding帧和第一个特征帧合为一帧，因此在第二层causal conv，我们就又需要上一个chunk的第二层输入的最后一帧，注意层与层之间的缓存是一一对应的。

通过如上分析，我们想实现conv层的缓存就有两种方法：

维护每一层conv输入的最后一帧缓存供下一个chunk使用，这样就不需要用额外overlap了，这样效率最高，缺点就是每一层都要维护，稍微麻烦一点
加入额外的overlap，让STFT重合的时间帧完全覆盖causal conv的padding范围。对于第一层conv输入，其需要上一chunk的最后一帧，对于第二层，则需要上一chunk的原本倒数第二帧（因为倒一已经被上一层conv融合了），这样算下来，标准DCCRN则需要大概13帧STFT时间帧的额外overlap（大概就是100ms），转换到时域上，这段overlap并不算短，并且这就要求当前chunk要大于100ms，这会引入额外的延迟。优点就是这样就不用专门去维护conv内存了，比较简单暴力

再额外提一嘴，LSTM的缓存是必须维护的，因为其包括了之前所有时间帧记忆

根据上面两种方法，就可以得到两种流式数据处理的方法：

小chunk法，仅考虑STFT的overlap，这样我们每个chunk就可以设的很小（极端地，可以直接把每个时间帧都当做一个chunk，也就是说直接对每个窗长的时域数据做FFT），STFT的overlap长度大概需要12.5ms（400*0.5/16000），我们就可以把chunk定为100甚至50ms以内，可以明显降低延迟

这样操作需要自定义一个STFT，因为padding需要自己替换。当使用torch的stft，center=false时，就不需要考虑overlap，直接拼接各chunk即可
大chunk法，通过大量overlap避免直接维护conv缓存，这样我们的overlap大概会来到100ms以上，因此chunk大小也会跟着变大，从而使延迟增加

这样看下来，小chunk法确实更好，如果我们有效chunk定为50ms，那除了第一个chunk，其他chunk实际长度则为50+12.5=62.5ms（或者说第一个chunk的12.5ms全为0padding），多出来的12.5ms会在转stft中消失，最终模型输出仍为50ms

自适应波束形成

发表于 2026-06-23 更新于 2026-06-24

自适应波束形成

波束形成的原理是调整相位阵列的基本单元参数，使得某些角度的信号获得相长干涉，而另一些角度的信号获得相消干涉。对各个麦克风信号加权求和、滤波，最终得到期望方向的语音信号，相当于形成一个“波束”。

传统波束形成是使用固定的权值来形成某个方向的波束，而自适应波束形成则是根据环境的变化，动态调整阵列麦克风的权重和相位，从而最大限度地抑制干扰噪声并增强目标信号。

基本数学模型

假设一个由 M 个阵元组成的任意构型的接收阵列，空间中存在一个来自方向 θ_0的期望信号，以及来自其他方向的干扰和环境噪声。则各麦克风接收到的信号（Mx1维）可表示为：

$X(t) = A(\theta_0)S(t) + \sum_{i=1}^{K} A(\theta_i)I_i(t) + N(t)$

其中A为导向矢量，0为目标方向，其他为干扰声源，N为环境噪声。

所有波束形成器的输出都可表示为:

$Y(t) = W^H X(t)$

其中权重W则根据各波束形成算法而不同

MVDR最小方差无失真响应

Minimum Variance Distortionless Response，最小方差-无失真响应MVDR的核心思想就是：在确保当前目标方向的信号增益为1（无失真）的约束条件（可以直接理解为使用导向矢量进行约束）下，最小化阵列输出的总功率（即最小化方差），可以表示为：

$\min_{W} \quad E[|Y(t)|^2] = W^H R_{xx} W \quad \text{s.t.} \quad W^H A(\theta_0) = 1$

拉格朗日乘子法求解最优问题

现在我们是要寻找功率函数在无失真约束下的最小值，则可以使用拉格朗日乘子法。

拉格朗日乘子法（Lagrange multipliers）是一种寻找多元函数在一组约束下的极值的方法。

对于求f(x)在g(x)约束下的极值，可以得到如下结论，在最优点x0处，f(x)和g(x)的梯度向量的方向必相同或相反，即存在一个拉格朗日乘子，使得

$f'(x_0)+\lambda g'(x_0)=0$

将该等式与约束条件联立即可求得乘子λ，将λ再带回上式即可获得x_0

在MVDR中，我们的拉格朗日目标函数可表示为：

$L(W, \lambda) = W^H R_{xx} W + \lambda (W^H A(\theta_0) - 1)$

根据上述结论，对W^H求偏导，有：

$\frac{\partial L}{\partial W^H} = R_{xx} W + \lambda A(\theta_0) = 0 \implies W = -\lambda R_{xx}^{-1} A(\theta_0)$

将该结果带入约束条件，最终可获得λ的表达：

$\lambda =- \frac{1}{A(\theta_0)^H R_{xx}^{-1} A(\theta_0)}$

带入W表达式，最终可获得MVDR最佳权重向量：

$W_{\text{MVDR}} = \frac{R_{xx}^{-1} A(\theta_0)}{A(\theta_0)^H R_{xx}^{-1} A(\theta_0)}$

这里就可以看出，W与协方差矩阵直接相关（分子是向量决定方向，分母是标量决定尺度，不能约）。最终功率谱则可表示为：

$P_{\text{MVDR}}(\theta) = \frac{1}{A(\theta)^H R_{xx}^{-1} A(\theta)}$

py代码

for f_idx, freq in enumerate(tqdm(freqs, desc="MVDRing")):
    X_f = data_valid[:,:, f_idx]	#提取当前频率分量
    R = (X_f.T @ X_f.conj())  #协方差矩阵
    eps = 1e-3 * np.trace(R)  #视情况选择大小，越大mvdr性能越强，鲁棒性越差
    R = R + eps * np.eye(self.Nchan) #对角加载
    R_inv = np.linalg.pinv(R)	#求逆
    a = np.exp(-1j * 2 * np.pi * freq * tdoa) #根据toda构建导向矢量
    Ra = R_inv @ a
    denominator = -np.sum(a.conj() * Ra, axis=0)  # a^H R^-1 a
    power_at_freq = 1.0 / np.real(denominator/ X_f.shape[0] + 1e-12) #功率谱

对同一段信号在相同频段进行mvdr和das功率扫描，可以看出MVDR在低频段上表现明显优于DAS，这是由于MVDR对干扰源的抑制能力极强，在低频段主瓣都很大的情况下也能区分声源位置

关于对角加载

MVDR中比较关键的一步就是需要求协方差矩阵R的逆，而从代码中我们可以看到，在求逆之前，我们进行了对角加载，那么对角加载有什么用呢？

核心痛点

在理想状态下，MVDR 依赖真实的协方差矩阵（数学期望）。但在实际工程中，我们只能通过有限个采样快照来估计。这就会导致一个问题：当快拍数较少（甚至小于阵元数），协方差矩阵会发生秩亏缺或条件数极高。此时对求逆会带来巨大的数值计算误差，导致波束形成器的权重剧烈抖动，旁瓣自适应抬高。

数学原理

因此，在估计的协方差矩阵的对角线上，我们加入了一个底噪声功率：（I为单位矩阵）

$R_{\text{DL}} = \hat{R}_{xx} + \sigma_{\text{DL}}^2 I$

加入这个底噪后，原本快拍数不足导致的最小特征值接近于0的问题（会使求逆输出极大值）会因为引入对角加载而被强行拉高特征值，令矩阵的条件数（最大特征值/最小特征值）大幅降低。矩阵求逆运算变得极度稳定，抑制了权重的随机剧烈抖动。

加载量的选择

对于添加的底噪值，其越小则代表自适应能力越强（太强的自适应会导致鲁棒性变差，对导向矢量准确性要求极高），当其接近于无限大时，此时MVDR则退化为DAS，因此需要选择一个合适的值

基于噪声基底

直接将加载量设定为系统背景噪声功率的 10到 100 倍

基于矩阵迹（trace）

$\sigma_{\text{DL}}^2 = \alpha \cdot \text{trace}(\hat{R}_{xx})/M,\alpha =10^{-3} \sim 10^{-1}$

根据对角线元素总和，即阵列所有阵元接收到的总功率，也就是将加载量也设置为一个根据采集信号动态调整的值

在我实际使用中，对角加载对mvdr结果影响很大，因此不能随便选

MVDR的优缺点

优点

高分辨率： 空间分辨能力显著优于传统的延迟求和波束形成。
强干扰抑制（抑制其他声源对目标声源干扰）： 能够自动在强干扰源方向形成零陷，大幅提升强干扰环境下的输出信噪比。
无需先验干扰信息： 只需要知道目标方向，不需要知道干扰的数量和方向。

缺点与工程挑战

对模型误差极其敏感（鲁棒性差，对导向矢量要求高）： 如果目标方向存在偏差（Pointing Error）或者阵元位置存在校准误差，MVDR会误将目标信号视为“干扰”并进行自适应压制，导致严重的自消散（Signal Canceling）现象。
相干源失效（快衰落）：当目标信号与干扰信号高度相关或相干（如同频多径反射）时，协方差矩阵 $R_{xx}$ 会趋于退化，导致算法失效。通常需要结合空间平滑（Spatial Smoothing）技术来解相干。
**计算量大：** 每次计算都需要对 $M \times M$ 的矩阵进行求逆操作，在阵元数 $M$ 较大或需要实时更新时，计算复杂度较高（$O(M^3)$）。
**样本快照数限制：** 在实际应用中，$R_{xx}$ 是通过有限的样本快照估计得到的（即 $\hat{R}_{xx} = \frac{1}{N}\sum X X^H$）。如果快照数 $N$ 小于阵元数 $M$，$\hat{R}_{xx}$ 将不可逆。

LCMV线性约束最小方差

线性约束最小方差(Linearly Constrained Minimum Variance, LCMV)波束形成器，是MVDR的扩展形式，其可表示为：

$\min_{W} \quad W^H R_{xx} W \quad \text{s.t.} \quad W^H C = F^H$

其实就是把原来的单个约束替换为了一组线性约束，C为MxL维的约束矩阵，每一列代表一个约束条件；F为Lx1维的增益控制向量，代表对应约束想达到的复增益

同样使用拉格朗日乘子法可以得到LCMV的最佳权重表达：

$W_{\text{LCMV}} = R_{xx}^{-1} C (C^H R_{xx}^{-1} C)^{-1} F$

当约束条件为单个导向矢量时，退化为MVDR，上式也完全等价于MVDR形式

常见约束矩阵设计

多方向无失真约束，即对两个方向同时进行无失真约束：
$C = [A(\theta_1), A(\theta_2)] ,F = [1, 1]^T$
强迫干扰零陷约束，即对已知干扰源直接约束其增益为0：
$C = [A(\theta_0), A(\theta_i)],F = [1, 0]^T$
导数约束，约束导向矢量对角度的一阶导数（甚至二阶导数）为 0，来使主瓣平坦化。
$C = \left[ A(\theta_0), \left. \frac{\partial A(\theta)}{\partial \theta} \right|_{\theta_0} \right],F = [1, 0]^T$
此时主瓣会增大，具备更强稳定性

GSC广义旁瓣相消器

GSC可以看作是LCMV的一种无约束等效实现，其示意图如下：

其可分为两部分，即上方的主路和下方的辅路，其中主路就是标准DAS波束形成，不再细说。

对于辅路，其任务是尽可能的还原信号中的噪声和干扰，并从DAS权重中减去，因此GSC的最终权重可表示为：

$\mathbf{w} = \mathbf{w}_q - \mathbf{B}^T \mathbf{w}_a$

与MVDR不同，其自适应权重用于计算干扰和噪声，而不是直接用于增强方向信号

阻塞矩阵B

其功能是各通道中的期望信号完全滤除（阻塞），仅允许噪声和干扰信号通过。因此，其与约束矩阵（导向矢量）是一定正交的，因此，当我们已经有比较准确的导向矢量a时，将各信号对齐后即可直接令B让相邻阵元相减来消除同相的目标信号：

$\mathbf{B} = \begin{bmatrix} 1 & -1 & 0 & 0 \\ 0 & 1 & -1 & 0 \\ 0 & 0 & 1 & -1 \end{bmatrix}$

这样得到的干扰和噪声是基于两两相邻阵元信号的，也可以通过均值计算：

$\mathbf{B} = \begin{bmatrix} 3/4 & -1/4 & -1/4 & -1/4 \\ -1/4 & 3/4 & -1/4 & -1/4 \\ -1/4 & -1/4 & 3/4 & -1/4 \end{bmatrix}$

这样每个阵元减去的都是所有阵元和的均值，理论上会更稳定一点点

自适应噪声相消器

信号通过阻塞矩阵后，只剩下干扰和噪声。自适应滤波器（如 LMS、RLS 或直接求逆算法）通过调整权重 w_a，使自适应支路的输出尽可能逼近静态支路中的干扰和噪声成分。

理论上在经过阻塞矩阵B后，信号中已不包含期望信号成分，因此自适应权重w_a的优化完全不会影响期望信号。基于这个前提，自适应部分的优化目标也就等同于使经过B后的信号输出总功率（方差）最小：

$\min_{\mathbf{w}_a} E\left[ |y(t)|^2 \right] = \min_{\mathbf{w}_a} (\mathbf{w}_q - \mathbf{B}^T \mathbf{w}_a)^H \mathbf{R}_{xx} (\mathbf{w}_q - \mathbf{B}^T \mathbf{w}_a)$

也就是最小化噪声和干扰，直接对目标函数求导找零点，即可得到无约束最优解：

$\mathbf{w}_a = (\mathbf{B} \mathbf{R}_{xx} \mathbf{B}^T)^{-1} \mathbf{B} \mathbf{R}_{xx} \mathbf{w}_q$

优缺点

主要优点就是把约束优化转化为了无约束的滤波问题，降低了计算复杂度，并且其物理和数学结构清晰，完全解耦；

缺点就是需要精准的导向矢量，否正阻塞矩阵B会产生泄露，自适应滤波器会将泄露的期望信号视作干扰，导致最后期望信号也受到抑制

为了防止自适应计算过强，也可以使用对角加载限制过大变化

与MVDR的关系

MVDR 与 GSC 等价，本质上是因为 GSC只是MVDR约束优化问题的一种结构化实现方式。

首先，对于GSC的权重表达，将其共轭转置后都右乘导向矢量后，即可得到MVDR的约束；如果将w表达式带入MVDR的目标函数（无约束格式），也可以直接得到GSC最终w_a的解

那GSC存在的意义是什么呢，大概是因为较早时期矩阵求逆运算量大，所以GSC实时性更好吧

GEV广义特征向量波束形成器

更多出现在语音增强中，GEV波束形成的核心思想是通过最大化输出信号的信噪比（SNR）或信号对干扰加噪声比（SINR）来推导最优导向矢量或空间滤波器系数，相比与MVDR要保证目标方向的无失真，GEV更加激进，选择直接最大化信噪比

对于带噪的信号，GEV直接将其视为：

$\mathbf{x}(f, t) = \mathbf{s}(f, t) + \mathbf{n}(f, t)$

所有噪声和干扰均被视为了一整个函数，GEV的目的就是让s和n最终信噪比尽可能大

构造目标与噪声协方差矩阵

根据上面分析，GEV的核心之一就是如何判断信号中哪些是目标信号，哪些是噪声。为了量化信号与噪声的空间分布，引入两个关键的空间相关矩阵：

$\mathbf{R}_{ss/nn}(f) = \frac {\sum_tM_{s/n}{(f,t)}x(f,t)x^H(f,t)}{\sum_tM_{s/n}{(f,t)}}$

其中，M为筛选协方差矩阵的mask，其形状与每个通道的stft结果x(f,t)相同，用于决定当前时-频点加权为噪声还是语音。也就是说，mask的计算是这一步的核心。

在GEV中，mask是一个[0,1]的概率，目前最主流的方法是通过神经网络进行预测，因此，我们在训练中需要同时拥有：

干净语音s
噪声n
混合信号x

而mask的训练目标则为s和n的一些算法结果，例如：

$m=\frac{|s|^2}{|s|^2+|n|^2}$

然后将预测的/hat{m}与mask进行loss计算

GEV优化目标

在获得了两个协方差矩阵后，GEV的优化目标如下：

$\max_{\mathbf{w}} \text{SNR}(\mathbf{w}) = \max_{\mathbf{w}} \frac{\mathbf{w}^{\mathsf{H}} \mathbf{R}_{ss} \mathbf{w}}{\mathbf{w}^{\mathsf{H}} \mathbf{R}_{nn} \mathbf{w}}$

该目标可以直接转化为广义特征值问题：

$\mathbf{R}_{ss} \mathbf{w} = \lambda \mathbf{R}_{nn} \mathbf{w}$

则权重w为最大特征值所对应的特征向量

后置滤波

GEV算出的特征向量具有任意的缩放因子（Scaling），这会导致输出信号产生频率失真（Speech Distortion）。为了解决这一盲源分离中的固有问题，通常引入 BAN （Blind Analytic Normalization）进行增益归一化：

$\mathbf{w}_{\text{BAN}}(f) = \frac{\mathbf{w}_{\text{GEV}}(f)}{\sqrt{\frac{1}{M} \mathbf{w}_{\text{GEV}}^{\mathsf{H}}(f) \mathbf{R}_{nn}(f) \mathbf{w}_{\text{GEV}}(f)}}$

优缺点

首先，我们可以很明显的看出，GEV在计算上完全不依赖导向矢量，也就是说其对阵列的几何结构没有要求，依靠数据驱动，而MVDR的质量几乎是由导向矢量决定，因此GEV理论上鲁棒性更高，由于其优化原理，信噪比也会更高

但是其也需要大量数据进行拟合，并且如果不使用BAN，会出现较大失真

麦克风阵列的一些理论补充

发表于 2026-06-15 更新于 2026-08-28

瑞利极限与远场判据

之前我们提到过判断声场模型的判据：

$R_r = \frac{2D^2}{\lambda}$

当声源与阵列的绝对距离大于这个值时即可认为是远场模型，此时不同麦克风接收信号的幅度差异较小，因此把不同麦克风采集的语音信号的幅值认为都是一样的，只需对各麦克风接收信号的相位差异进行处理即可。声波视为平面波。

并且从公式中我们可以得知波长越小（频率越高），这个判据就会越远，即我们需要更远的距离才能将当前频率下的声场视作远场处理，因此若想使用更简单的远场模型，则要求我们的频段选择偏小

但是这又带来了一个问题，即瑞利极限带来的空间分辨率问题。

$Δ\theta_R=1.22\frac {\lambda}{D}$

瑞利极限最先在光学中提出，在声学成像中，我们可以理解为当两个点源小于这个距离时，系统则无法分辨他们的准确位置，在我个人的使用中，会出现两个声源在中点合成为一个源峰的现象（主要是DAS和SPR-phat）。

也就是说，当阵列尺寸不变时，我们如果想获得较高的空间分辨率，则要求波长偏小，也就是频段选择偏高，这与我们在远场判据中的要求是冲突的

将上两式联立我们可以得到：

$R_r=\frac{2.44D}{Δ\theta_R}$

也就是说阵列分辨率越高，则远场要求更远

因此，在实际成像中，我们应该优先满足高分辨率，则频段选择应该略高一些；即使采集距离较近，也不该使用较低频段强行拉低远场判据，因为这必然会导致分辨率急剧下降，此时就算可以用远场模型，也是没意义的

综上，若想获得质量较高的成像，在工程中首选是提高孔径D与取到合适的频段

主瓣宽度问题

刚刚我们提到，当空间分辨率不足时，也就是波长较长，频率较低时，会出现声源相互合成的现象（或者其他无法准确显示声源位置的现象）。其中很大一部分原因就是声源的主瓣宽度与空间分辨率成强正比关系，也就是说合成的峰值其实就是两个声源的主瓣叠加后超过了各自的主瓣中心强度

也就是说，在不改变孔径的条件下，频率越低，主瓣越宽。当然，主瓣宽度还与孔径大小、阵列数量和使用的定位算法有关，这一节我们主要就讲频率的问题

其实很好理解，在各类利用相位差的算法中（SRP-PHAT，DAS），我们在频域对相位的表达是：

$Δ\phi=2\pi f \tau$

也就是说，f越小，相位差越不明显，最终导致功率图峰值附近都很平滑，变化很慢，就会呈现出宽主瓣的现象

MUSIC为什么可以突破瑞利极限

其实基于上一节主瓣宽度的分析，再结合music的原理，就很好理解了。对于DAS，其定位图像是利用的波束功率，因此非常受限于波长和孔径；而对于SRP-PHAT，由于其舍弃了幅值，相比于DAS会稍微好一点，但是其本质还是基于相位差的互相关功率，也会被上述的主瓣问题影响。

而对于MUSIC，其是直接利用的导向矢量和噪声子空间的正交关系，与相位差无关，因此其空间分辨率理论上可以自己决定其扫描步长，而不是受限于瑞利极限，从而实现超分辨率成像

虽然但是，实际还是会与频率有影响，因为导向矢量中也包含了相位信息，但是没有去专门计算差值了；但实际中其分辨率还是会达到瑞利极限的数倍。

关于TDOA，导向矢量和波束权重的定义习惯

TDOA

本博客中使用的TODA均为下述函数

def get_tdoa(self, phi, theta):
    # phi=0为+x，=180为-x
    # theta=0为+z =180为-z

    a = np.array([  np.sin(theta) * np.cos(phi) ,
                    np.sin(theta) * np.sin(phi),
                    np.cos(theta)])

    tdoa = np.sum(self.mic_struct * a[None, :], keepdims=True, axis=1) / self.sound_speed
    return -tdoa

其中，a代表了从0点指向声源的方向（正负问题的核心，这点ai都容易搞错，需要注意），但实际信号是从声源指向阵列，因此这里得到的tdoa，正值则代表了信号提前于参考点到达，因此，为了和后续标准公式中的τ接轨，该函数最终应返回负值，即：

$\tau = -tdoa$

τ为正代表信号晚到

导向矢量

再说完了上述tdoa的定义后，我们来说导向矢量。现在大部分的导向矢量都定义为：

$a=e^{-j2\pi f\tau}$

而这里的τ则是代表信号晚到了多少（即晚到信号τ为正），现在大部分文献都是用的这种格式（虽然本人觉得tdoa就带值然后导向矢量不带负号更合理），因此本博客中所有导向矢量也采用这种形式

波束权重

最阴的来了，现在对于波束形成，最常见的公式是：

$Y(t) = W^H X(t)$ $\mathbf{R}_{xx} = E[\mathbf{X}(t)\mathbf{X}^H(t)] \approx \frac{1}{N} \sum_{t=1}^{N} \mathbf{X}(t)\mathbf{X}^H(t)$ $P = W^H R_{xx} W$

我其实一直就在想W的这个共轭转置是何意味，按照我的理解来说W就应该直接代表了波束形成的最终权重。但现在大部分的定义就是W形式与导向矢量相同，即

$w=e^{-j2\pi f\tau}$

因此对于DAS来说，若要对齐相位则需要对W取共轭来抵消导向矢量中的相位，这真的太反人类了，也就是说在代码里W和导向矢量a甚至是相同的值。而在各DOA，波束形成算法中，使用的大都也都是上述的导向矢量和权重格式，即默认：

导向矢量指数默认带负号
W^H才代表波束形成最终权重

同时还需要注意的就是，对于协方差矩阵也要配套为X乘X^H，千万不要写反了