【语言信号增强算法研究-1】维纳滤波（Wiener Filter）

1 语音增强方法分类

2 维纳滤波的局限性

对于非线性和非高斯噪声的处理效果不佳；对于信号和噪声的统计特性要求比较高，需要准确地了解信号和噪声的分布规律（说明自适应很差）；在处理复杂信号时，需要进行多次计算和优化，计算量较大。

3 维纳滤波的原理

4 实现源码与程序

4-1 维纳滤波器wiener_filter函数的实现代码

%% ============================================================ %%
%%                      维纳滤波器函数实现
%% ============================================================ %%

%% ------------------------------------------------------------ %%
%% param 
%  noisy：含噪信号  
%  clean：纯净信号(无噪信号)  
%  noise：纯噪音信号 
%  para ：(集总了滤波器各参数的)结构体

%% retval
%  H：滤波器系数
%  enhenc：增强后的(时域形式)信号

%% (集总了滤波器各参数的)结构体para成员介绍：
%  n_fft      : FFT
%  hop_length :
%  win_length :
%  alpha      :
%  beta       :
%% ------------------------------------------------------------ %%
function [H, enhenc] = wiener_filter(noisy, clean, noise, para) % 使用维纳滤波函数
    %% 从结构体获取到相关参数
    n_fft      = para.n_fft;
    hop_length = para.hop_length;
    win_length = para.win_length;
    alpha      = para.alpha;
    beta       = para.beta;

    %% 使用短时傅里叶变换（STFT）后取绝对值，并转换为相应的频谱表示
    S_noisy = abs(stft(noisy, n_fft, hop_length, win_length)); % 得到含噪信号频谱
    S_noise = abs(stft(noise, n_fft, hop_length, win_length)); % 得到纯噪音信号频谱
    S_clean = abs(stft(clean, n_fft, hop_length, win_length)); % 得到无噪信号频谱

    %% 计算(无噪/纯噪音)信号的频谱平方后的均值
    Pxx = mean(S_clean.^2, 1);
    Pnn = mean(S_noise.^2, 1);

    %% 根据维纳滤波器的公式计算滤波器系数H
    H = (Pxx ./ (Pxx + alpha * Pnn)) .^ beta;   

    %% 将噪声信号频谱与滤波器系数相乘，得到增强后的频谱S_enhec
    S_enhec = S_noisy .* H;                     

    %% 使用逆短时傅里叶变换（ISTFT）将增强后的频谱S_enhec转换回时域信号，得到增强后的信号enhenc。
    enhenc = istft(S_enhec, hop_length, win_length);
end

4-2 滤波效果仿真代码

援引CSDN上一篇blog的代码改写如下

因为维纳滤波器的参数求解过程中必须要用到参考信号

仿真信号 $S g n$ 为两个正弦信号的叠加： $Sgn=A*sin(2πf_1t)+B*sin(2πf_2t)$

close all;
clear;
clc;

%% 信号产生，对原始信号进行采样
A=10;
B=15;
f1=1000;
f2=2000;
fs=10^5;
t=0:999;                            % 采样点t = [0,999],长度1000
M = length(t);                      % 信号长度
Sgn = A*sin(2*pi*f1*t/fs) + B*sin(2*pi*f2*t/fs);

%% 给仿真信号加入高斯白噪声
SNR = 3;                            % 初始信噪比
Obs_Sgn = awgn(Sgn,SNR,'measured'); % 给仿真信号加入信噪比为-3dB的高斯白噪声
WGN = Obs_Sgn - Sgn;                % 观测信号Obs_Sgn = 仿真信号Sgn + 高斯白噪声WGN(误差信号，每次运行都不一样)

%% 第一种情况：(期望)感兴趣信号SOI为的原信号Sgn
SOI = Sgn; 

%% 维纳滤波
N = floor(length(Obs_Sgn)*0.1);     % 滤波器的阶数，向下取整
Rxx=xcorr(Obs_Sgn, N-1, 'biased');  % 自相关函数1*(2N-1)维度，返回一个延迟范围在[-N，N]的互相关函数序列

%% 变成矩阵(N*N维度)
for i=1:N
    for j=1:N
        mRxx(i,j)=Rxx(N-i+j); % N*N维度
    end
end

%% 产生维纳滤波中x方向上观测信号与期望信号d的互相关矩阵
Rxd=xcorr(Obs_Sgn, SOI, N-1, 'biased');        % 互相关函数1*(2N-1)维度

%% 变成矩阵1*N维度
for i=1:N
    mRxd(i)=Rxd(N-1+i); % 1*N维度
end

h = inv(mRxx)*mRxd';                % 由wiener-Hopf方程得到滤波器最优解, h是N*1维度

%% 检验wiener滤波效果
y   = conv(Obs_Sgn, h);             % 滤波后的输出,长度为M+N-1，要截取前M个。
y   = y(1:M);                       % yy = filter(h,1,x);  % 用卷积或者直接用filter都可以
Py  = sum(power(y,2))/M;            % 滤波后信号y的功率
ERR = SOI - y;                      % 输出减去期望等于滤波误差
Pn  = sum(power(ERR,2))/M;          % 滤波后噪声功率
Delta_P = Py-Pn                     % 功率增量
SNR_Filtered = 10*log10(Delta_P/Pn) % 滤波后信噪比(单位db) 
imv = 10*log10(Delta_P/Pn/power(10,SNR/10)) % 滤波后较滤波前信噪比提高了imv dB。

figure(1), subplot(311);
plot(t,Sgn);
title('Signal原信号');

subplot(312);
plot(t,WGN);
title('WGN高斯白噪声');

subplot(313)
plot(t,Obs_Sgn)
title('观测波形')

在这里插入图片描述

%% 期望和滤波后的信号对比
figure(2), subplot(211);
plot(t, SOI, 'r:', t, y, 'b-','LineWidth',1);
legend('期望信号','滤波后结果'); 
title('期望信号与滤波结果对比');
xlabel('观测点数'), ylabel('信号幅度');
axis([0 1000 -50 50]);

subplot(212), plot(t, ERR);
title('输出误差');
xlabel('观测点数'), ylabel('误差幅度');
axis([0 1000 -50 50]);


%% 滤波前后对比
figure(3), subplot(211);
plot(t, Obs_Sgn);
title('维纳滤波前');
xlabel('观测点数'), ylabel('信号幅度');
axis([0 1000 -50 50])

subplot(212), plot(t, y);
title('维纳滤波后');
xlabel('观测点数'), ylabel('误差幅度');
axis([0 1000 -50 50]);

在这里插入图片描述

%% 情况二：(期望)感兴趣信号SOI为的加性高斯白噪声WGN
SOI = WGN;

...
%% 期望信号与滤波结果对比
figure(2)
plot(t, SOI, 'r:', t, y, 'b-','LineWidth',1);
legend('期望信号','滤波后结果'); 
title('期望信号与滤波结果对比');
xlabel('观测点数'), ylabel('信号幅度');
axis([0 1000 -50 50])


%% Signal原始信号与噪声抵消后结果对比
figure(3)
plot(t, Sgn, 'r:', t, Obs_Sgn - y, 'b-','LineWidth',1);
legend('Signal原始信号','噪声抵消后结果');
title('Signal原始信号与噪声抵消后结果对比');
xlabel('观测点数'), ylabel('信号幅度');
axis([0 1000 -50 50])

在这里插入图片描述

%% 情况三：(期望)感兴趣信号SOI为的观测信号
%% 观测信号与滤波结果对比
figure(2), subplot(211);
plot(t, Obs_Sgn, 'r:', t, y, 'b-','LineWidth',1);
legend('观测信号','滤波后结果'); 
title('观测信号与滤波结果对比');
xlabel('观测点数'), ylabel('信号幅度');
axis([0 1000 -50 50])

subplot(212);
plot(t, ERR);
title('输出误差');
xlabel('观测点数'), ylabel('误差幅度');
axis([0 1000 -50 50]);


%% 滤波前后对比
figure(3), subplot(211);
plot(t, Obs_Sgn);
title('维纳滤波前');
xlabel('观测点数'), ylabel('信号幅度');
axis([0 1000 -50 50])

subplot(212), plot(t, y);
title('维纳滤波后');
xlabel('观测点数'), ylabel('误差幅度');
axis([0 1000 -50 50]);

在这里插入图片描述

4-3 仿真实验结论

（期望）感兴趣信号设置为原始信号 $S g n$ 时，效果最好

5 详细说明

原始信号 $S g n$ 经过噪声信道(加入白噪声 $W GN$ )，变成含噪信号 $O b s S g n$ （ $O b s S g n = S g n + W GN$ ）。
$O b s S g n$ 经过维纳滤波器(Wiener Filter)后得到 $S g n F i lt ere d$ ，由于期望的感兴趣信号为 $SO I$ ，则误差信号 $ERR = SO I - S g n F i lt ere d$ ，滤波器系数为H
按照均方误差最小原则，从 $O b s S g n$ 信号中剥离得到 $S g n$ ，叫维纳滤波理论
不能使用维纳滤波法对未知信号进行滤波！！！因为设计滤波器系数时需要使用到期望信号 $SO I$