【生物医学信号处理及其MATLAB应用】Chapter 6 维纳滤波

设有一个线性系统，它的单位脉冲响应是$h(n)$，当输入一个观测到的随机信号$x(n)$，简称观测值，且该信号包含噪声$w(n)$和有用信号$s(n)$，简称信号，也即

​ $x(n)=s(n)+w(n)$

​ 则输出$y(n)$为：

​ $y(n)=x(n)\ast h(n)={\sum }_{m=-\infty }^{+\infty }h(m)x(n-m)$

​ 我们希望输出得到的$y(n)$与有用信号$s(n)$尽量接近，因此称$y(n)$为$s(n)$的估计值，用$\widehat{s(n)}$来表示$y(n)$，就有了维纳滤波器的系统框图，这个系统的单位脉冲响应也称为对于$s(n)$的一种估计器。

​ 如果该系统是因果系统，上式中的$m=0,1,2,...$，则输出的$\hat{s}(n)$可以看成是由当前时刻的观测值$x(n)$和过去时刻的观测值$x(n-1)、x(n-2)、x(n-3)......$的估计值。用当前的和过去的观测值来估计当前的信号$y(n)=\hat{s}(n)$，称为滤波；用过去的观测值来估计当前的或将来的信号$y(n)=\hat{s}(n+N)(N\ge 0)$，称为预测；用过去的观测值来估计过去的信号$y(n)=\hat{s}(n-N)(N\ge 1)$，称为平滑或者内插。

​ 从系统中估计得到的$\hat{s}(n)$信号和我们期望得到的有用信号$s(n)$不可能完全相同，这里用$e(n)$来表示真值和估计值之间的误差：

​ $e(n)=s(n)-\hat{s}(n)$

​ 显然，$e(n)$是随机变量，维纳滤波的误差准则就是最小均方误差准则：

​ $E[e^2(n)]=E[(s(n)-\hat{s}(n){)}^2]$

​ 维纳滤波是解决线性滤波和预测问题的方法，且是以均方误差最小为准则的，只适用于平稳随机过程。

6.1 维纳滤波器的时域解

​ 设计维纳滤波器的过程就是寻求在最小均方误差下滤波器的单位脉冲响应$h(n)$或传递函数$H(z)$的表达式，其实质就是解维纳-霍夫（$Wiener-Hopf$）方程。

6.1.1 因果的维纳滤波器

​ 设$h(n)$是物理可实现的，也即是因果序列：

​ $h(n)=0,n<0$

​ 因此，

​ $y(n)=\hat{s}(n)={\sum }_{m=0}^{+\infty }h(m)x(n-m)$

​ $E[e^2(n)]=E[(s(n)-{\sum }_{m=0}^{+\infty }h(m)x(n-m){)}^2]$

​ 要使得均方误差最小，则将上式对各$h(m)(m=0,1,...)$求偏导，并且等于零，得：

​ $2E[(s(n)-{\sum }_{m=0}^{+\infty }{h}_{opt}(m)x(n-m))x(n-j)]=0j=0,1,2,...$

​ 即：

​ $E[s(n)x(n-j)]={\sum }_{m=0}^{+\infty }{h}_{opt}(m)E[x(n-m)x(n-j)]j\ge 0$

​ 用相关函数$R$来表达上式，则得到维纳-霍夫方程的离散形式：

​ $R_{xs}(j)={\sum }_{m=0}^{+\infty }{h}_{opt}(m)R_{xx}(j-m)j\ge 0$

​ 从维纳-霍夫方程中解出的$h$就是最小均方误差下的最佳$h$：$h_{opt}(n)$。

6.1.2 有限脉冲响应法求解维纳-霍夫方程

​ 设$h(n)$是一个因果序列且可以用有限长（$N$点长）的序列去逼近，则

​ $y(n)=\hat{s}(n)={\sum }_{m=0}^{N-1}h(m)x(n-m)$

​ $E[e^2(n)]=E[(s(n)-{\sum }_{m=0}^{N-1}h(m)x(n-m){)}^2]$

​ $2E[(s(n)-{\sum }_{m=0}^{N-1}{h}_{opt}(m)x(n-m))x(n-j)]=0j=0,1,2,...,N-1$

​ $E[s(n)x(n-j)]={\sum }_{m=0}^{N-1}{h}_{opt}(m)E[x(n-m)x(n-j)]j=0,1,...,N-1$

​ $R_{xs}(j)={\sum }_{m=0}^{N-1}{h}_{opt}(m)R_{xx}(j-m)j=0,1,2,...,N-1$

​ 于是得到$N$个线性方程：
$$\{\begin{array}{l}j=0,R_{xs}(0)=h(0)R_{xx}(0)+h(1)R_{xx}(1)+...+h(N-1)R_{xx}(N-1)\\ j=1,R_{xs}(1)=h(0)R_{xx}(1)+h(1)R_{xx}(0)+...+h(N-1)R_{xx}(N-2)\\ ...\\ j=N-1,R_{xs}(N-1)=h(0)R_{xx}(N-1)+h(1)R_{xx}(N-2)+...+h(N-1)R_{xx}(0)\end{array}$$
​ 写成矩阵形式有：

​ $ \begin{bmatrix} R_{xx}(0) & R_{xx}(1) & … & R_{xx}(N-1) \ R_{xx}(1) & R_{xx}(0) & … & R_{xx}(N-2) \ \vdots & \vdots & \cdots & \vdots \ R_{xx}(N-1) & R_{xx}(N-2) & \cdots & R_{xx}(0) \end{bmatrix} \begin{bmatrix} h(0) \ h(1) \ \vdots \ h(N-1) \end{bmatrix}=\begin{bmatrix} R_{xs}(0) \ R_{xs}(1) \ \vdots \ R_{xs}(N-1) \end{bmatrix} $

​ 简化形式：

​ $R_{xx}H=R_{xs}$

​ 式中，$H=[h(0),h(1),...,h(N-1){]}^{\prime }$，是待求的单位脉冲响应；

​ $R_{xs}=[R_{xs}(0),R_{xs}(1),...,R_{xs}(N-1){]}^{\prime }$，是互相关序列；

​ $R_{xx}$是自相关矩阵。

​ 只要$R_{xx}$是非奇异的，就可以求到$H$：$H={R_{xx}}^{-1}{R}_{xs}$

​ 求到$h_{opt}(n)$后，这时的均方误差为最小。

6.1.3 预白化法求解维纳-霍夫方程

​ 随机信号都可以看成是由一白色噪声$w_1(n)$激励一个物理可实现的系统或模型的响应，用白化法求解维纳-霍夫方程的步骤如下：

​ ①对观测信号$x(n)$的自相关函数$R_{xx}(m)$求$Z$变换得到$R_{xs}(z)$；

​ ②利用等式$R_{xx}(z)={\sigma }_{w_1}^{2}B(z)B(z^{-1})$，找到最小相位系统$B(z)$；

​ ③利用均方误差最小原则求解因果的$G(z)$；

​ ④$H(z)=\frac{G(z)}{B(z)}$，即得到维纳-霍夫方程的系统函数解

6.2 MATLAB应用

6.2.1 非因果维纳滤波器的实现

function [yhat, H] = wienerFilter(desired, observation, graphicsFlagOn, Fs)

% 检测输入信号
if nargin < 4
    Fs = 256;
end
Fs = round(Fs);
if nargin < 3
    graphicsFlagOn = 1;
end
% 估计噪声信号，观察信号=期望信号+噪声信号
noise = observation - desired;

% 输入和噪声信号功率谱估计
Sf = pwelch(desired, [], [], [], [], 'twosided');
Nf = pwelch(noise, [], [], [], [], 'twosided');

% 估计维纳滤波器的幅频响应
H = Sf./(Sf + Nf);

% 计算滤波器的单位脉冲响应
h = ifft(H);

% 用单位脉冲响应和fftfilt函数进行比较
yhat = fftfilt(h, observation);

% 绘制图像
if graphicsFlagOn = 1
    plotWiener(desired, noise, observation, yhat, Fs);
end

end

function plotWiener(model, noise, observation, utilda, Fs, nwindow)

if nargin < 5
    error('对plotWiener的函数调用不当');
end
if nargin < 6
    nwindow = 512;
end

[Psig, Fsig] = pwelch(observation, nwindow, [], [], Fs);
[Pm, Fm] = pwelch(model, nwindow, [], [], Fs);
[Pn, Fn] = pwelch(noise, nwindow, [], [], Fs);
[Pw, Fw] = pwelch(utilda, nwindow, [], [], Fs);

% 绘制时域图
figure;
t=[1/Fs:1/Fs:length(observation)/Fs];
subplot(3, 1, 1); 
plot(t, observation); 
ylabel('观察信号');title('Wiener Filter');
subplot(3, 1, 2);
plot(t, model);
ylabel('真实信号');
subplot(3, 1, 3);
plot(t, utilda);
ylabel('输出信号');xlabel('时间/s');
if 0
    figure; 
    plot(t, observation);hold on;
    plot(t, utilda, 'r');
    ylabel('观察信号与维纳滤波信号的比较');
    xlabel('时间/s');
end

6.2.2 维纳滤波器FIR方法的实现

function [hopt, error] = wein(x, s, N)
% FIR方法实现维纳滤波器

rxx = xcorr(x, 'biased');
rss = xcorr(s, 'biased');

rx = rxx(length(x):end);
Rxx = toeplitz(rx(1:N));

rs = rss(length(x):end);
Rss = toeplitz(rs(1:N));

hopt = inv(Rxx)*rs(1:N);
error = rs(1) - rs(1:N)*hopt;

end

function [hopt, error] = wein(x, s)
% 寻找最优阶数
[h, error(1)] = wein(x, x, 1);
for i = 2:100
    [h, error(i)] = wein(x, x, i);
    if abs(error(i) - error(i-1)) < 0.0001
        N = i - 1;
        break
    end
end
[hopt, error] = wein(x, s, N);

6.2.3 利用维纳滤波算法对语音信号去噪

clear;clc;close all
% 读入语音信号
[data, fs] = audioread('clouds.wav');
data = data([10*fs:100:20*fs], 1);
fs = fs/100;
N = length(data);
t = (0:N-1)/fs;
% 画语音信号时域图
figure;
plot(t, data);
xlabel('时间/s');ylabel('幅度/uV');
title('语音信号时域图');

% 添加随机噪声
xn = randn(1, N)';
x = data + xn;

% 利用非因果维纳滤波器降噪
[yhat, H] = wienerFilter(data, xn, 1, fs);
% 利用因果维纳滤波器降噪
[y, e] = causal_wiener_filter(x, data, 8);
% 绘制输入信号、期望输出、噪声和维纳滤波器的输出
figure;
subplot(3, 1, 1);plot(x);
title('输入信号');
subplot(3, 1, 2);plot(data);
title('期望输出');
subplot(3, 1, 3);plot(y);
title('因果维纳滤波器的输出');

% 调整信噪比观察结果变化
figure;
subplot(3, 2, 1);plot(data);
title('未加噪声的原始信号');
snr = [10 5 0 -5 -10];
w = sqrt(var(data)./(var(x).*(10.^(snr./10))));
for i = 1:length(snr)
    x1(:, i) = data + w(i)*xn;
    % 利用因果维纳滤波器降噪
    [y, e] = causal_wiener_filter(x1(:, i), data, 8);
    str = ['信噪比', num2str(snr(i)), 'dB'];
    subplot(3, 2, i+1);
    plot(y);
    title(str);
end
for j = 1:length(snr)
    x1(:, j) = data + w(j)*xn;
    % 利用非因果维纳滤波器降噪
    [yhat, H] = wienerFilter(data, w(j)*xn, 1, fs);
end
% 以因果维纳滤波器为例，观察信噪比与相关系数的关系
for k = 1:length(snr)
    x1(:, k) = data + w(k)*xn;
    corr_sx(k) = corr(data, x1(:, k));
end
figure;
plot(snr, corr_sx);
xlabel('信噪比');ylabel('相关系数');
title('相关系数与信噪比的关系图');

function [y, e] = causal_wiener_filter(x, d, order)
% x: 输入信号
% d: 期望输出信号
% order: 滤波器阶数

% 该函数接受输入信号x、期望输出信号d和滤波器阶数order作为输入，并返回滤波器的输出信号y和误差信号e。函数首先计算输入信号的自相关矩阵和交叉相关向量，然后使用矩阵求逆运算计算滤波器系数。最后，函数应用滤波器并计算输出信号和误差信号。

% 初始化变量
N = length(x);
y = zeros(N, 1);
e = zeros(N, 1);
w = zeros(order+1, 1);

% 计算自相关矩阵和交叉相关向量
R = zeros(order+1, order+1);
p = zeros(order+1, 1);
for n = order+1:N
    x_vec = x(n:-1:n-order);
    R = R + x_vec * x_vec';
    p = p + x_vec * d(n);
end

% 计算滤波器系数
w = R \ p;

% 应用滤波器
for n = order+1:N
    x_vec = x(n:-1:n-order);
    y(n) = w' * x_vec;
    e(n) = d(n) - y(n);
end

运行结果：

​ （暂未编写维纳预测器的Matlab代码，维纳预测器相关内容请参考 饶妮妮，李凌. 生物医学信号处理[M]，成都，电子科技大学出版社，2005.7）