Introduction to Wiener Filter

维纳滤波信号领域最基本，最常见的滤波方法之一，陈老师的麦克风阵列处理一书中首先介绍的就是维纳滤波。本篇文章中我简单的介绍一下维纳滤波。

因为是从ASR看到前端来的，因此，对维纳滤波的理解做不到十分深入。在我看来，维纳滤波的典型特征是MMSE准则，由此就带来一个困惑，没有参考信号或者目标信号，如何做均方误差最小化。换句话说，很多时候，如果已经知道参考信号了，那何须继续滤波呢……

后来我突然想到手机上大部分人容易忽略的一个器件，辅助降噪麦克风（就是那个小圆圈），才对这个疑惑有了解释，在使用维纳滤波的系统中，多少是存在办法获取参考信号的，辅助降噪麦克风的功能不出意外的话，就是采集环境噪声，给系统提供参考信号的统计量。有时间我觉得有必要研究一下webrt中的维纳滤波算法，听说是目前效果最好的实现。

针对单通道的denoise任务，时域上MSE准则写成：

$$\mathcal{J}(\mathbf{h}) = E\left[\left(z(k) - \mathbf{h}^\top\mathbf{y}(k)\right)^2\right] \tag{1}$$

其中：

$$\mathbf{y}(k) = [y(k), y(k -1), \cdots, y(k - L + 1)]$$

$\mathbf{h}$表示要估计的滤波器系数（有限冲击响应，FIR），$L$表示$\mathbf{h}$的阶数。

解目标函数$(1)$可以得到解析解：

$$\mathbf{h}_W = \arg \min_{\mathbf{h}} \mathcal{J}(\mathbf{h}) = \mathbf{R}_{yy}^{-1} \mathbf{r}_{yx}$$

其中$\mathbf{R}_{yy}$表示观测信号的相关矩阵：

$$\mathbf{R}_{yy} = E[\mathbf{y}(k)\mathbf{y}(k)^T ]$$

$\mathbf{r}_{yx}$表示期望信号和观测信号的协相关矩阵：

$$\mathbf{r}_{yx} = E[\mathbf{y}(k) x(k)]$$

在噪声和期望信号相互独立的假设下，$\mathbf{r}_{yx} $可以写成：

$$\mathbf{r}_{yx} = \mathbf{r}_{yy} - \mathbf{r}_{vv} = E[\mathbf{y}(k) y(k)] - E[\mathbf{v}(k) v(k)] \tag{2}$$

以降噪麦克风为例，可以获取参考噪声，因此估计$\mathbf{r}_{vv}$十分方便。如果要直接计算$\mathbf{r}_{yx}$，就要得到目标信号参考，这个难度显然要大很多，因此，多用$(2)$式间接计算。在我后面的实验中，取起始的一段噪声信号计算$\mathbf{r}_{vv}$，$\mathbf{r}_{yy}, \mathbf{R}_{yy}$根据观测信号均可直接算出。

下面用一个toy简单验证一下，噪声类型为白噪声

wav = 'wav/egs1.wav';
snr = 1;
filter_size = 512;
warm_up = 51200;

speech = audioread(wav);
M = norm(speech, inf);
% modify matlab's awgn
[noisy, noise] = my_awgn(speech, snr, 'measured');
T = length(noisy);

warm_up = min(warm_up, T);
fprintf("using first %d samples to estimate rvv, totally %d samples...\n", warm_up, T);

padding = zeros(1, filter_size - 1);
Y = toeplitz([noisy' padding], [noisy(1), padding]);
V = toeplitz([noise' padding], [noise(1), padding]);

Ryy = Y' * Y / T;

ryy = mean(Y(1: T, :) .* noisy, 1);
rvv = mean(V(1: warm_up, :) .* noise(1: warm_up), 1);
ryx = ryy - rvv;

Hw = Ryy \ ryx';

denoise = Y * Hw;
denoise = denoise(1: T, :);

subplot(3, 1, 1);
plot(speech);title("Clean Speech");
xlim([1, T]); ylim([-1, 1]);
subplot(3, 1, 2);
plot(noisy / norm(noisy, inf) * M);title("Noisy Speech");
xlim([1, T]); ylim([-1, 1]);
subplot(3, 1, 3);
plot(denoise / norm(denoise, inf) * M);title("Wiener Denoised Speech");
xlim([1, T]); ylim([-1, 1]);

测试的时候其实发现denoise效果一般（如果上面的toy实现没有什么大问题的话）。下面给出一个$\text{SNR}=2$的例子（还是可以看出噪声部分得到了部分削弱）。