之前的文章講了使用經典濾波器來抑制噪聲。

裡面提到，“用經典濾波器抑制噪聲，非常簡單。如果噪聲的功率譜PSD和有用訊號功率譜PSD沒有重疊的話，那可以實現非常好的效果。

但是，如果有重疊，去噪的效果就不是特別理想了。因為在復指數訊號空間裡面，沒辦法把有用訊號和噪聲訊號分離啊。”

正是由於“噪聲譜和有用訊號譜可能重疊”，所以發展了維納濾波器。

前面的文章對維納濾波器的設計也講過了。

這篇文章，就是真實地來操作一下，設計一個維納濾波器來抑制噪聲。

因為沒有去錄音，所以噪聲源還是matlab裡的randn產生的高斯過程的資料。

再假設高斯過程並不是直接加入有用資料的，而是經過了一個“通道”，發生了一些變化，比如，AR過程。

然後就是分兩步了。

1. 訓練過程，用已知的訓練訊號和已知的接收訊號，通過解方程，求得濾波器係數；

2. 濾波器係數不變，用這組係數對此後接收到的訊號進行濾波，實現噪聲抑制。

程式碼比較簡單，所以直接出來了。

主函式。

%% 維納去噪，基於訓練的
% 通過解方程求濾波器係數
% 作者：qcy

% 版本說明：補充濾波階段
% 版本：v1.1 
% 時間：2016年10月29日18:18:21

% 版本說明：訓練階段
% 版本：v1.0
% 時間：2016年10月29日16:52:37

clear;
close all;
clc

%% 讀入音訊
filedir=[];                             % 設定路徑
filename='bluesky1.wav';                % 設定檔名
fle=[filedir filename];                 % 構成完整的路徑和檔名
[s, fs] = audioread(fle);               % 讀入資料檔案

s=s-mean(s);                            % 消除直流分量
s=s/max(abs(s));                        % 幅值歸一
N=length(s);                            % 語音長度
time=(0:N-1)/fs;                        % 設定時間刻度

%% 產生噪聲並加入
SNR = 0;                                  % 設定信噪比
r2=randn(size(s));                      % 產生隨機噪聲
b=fir1(31,0.5);                         % 設計FIR濾波器,代替H
r21=filter(b,1,r2);                     % FIR濾波
% [signal_with_noise,noise]=add_noisedata(s,data,fs,fs1,snr);
[r1,r22]=add_noisedata(s,r21,fs,fs,SNR);% 產生帶噪語音，信噪比為SNR　

%% 解方程
% 【訓練階段】
h_length = 100; % 這個東西我目前所知道的，就只有憑感覺去設定了。
desired_signal = s;
observed_signal = r1;
h = my_weiner_filter1( h_length,desired_signal,observed_signal );

% 濾波
% 用維納濾波器，作用在接收到的訓練訊號上。看看效果-->這還是屬於訓練階段
y = filter(h,1,r1);
output = y;

snr1=SNR_singlech(s,r1);                % 計算初始信噪比
snr2=SNR_singlech(s,output);            % 計算濾波後的信噪比
snr=snr2-snr1;

%% 畫圖
figure;
subplot 311; plot(time,s,'k'); ylabel('幅值');
ylim([-1 1 ]); title('原始語音訊號');
subplot 312; plot(time,r1,'k'); ylabel('幅值') 
ylim([-1 1 ]); title('帶噪語音訊號');
subplot 313; plot(time,output,'k'); 
ylim([-1 1 ]); title('濾波輸出語音訊號');
xlabel('時間/s'); ylabel('幅值')

%% 列印SNR
fprintf('[訓練]\n',snr1);
fprintf('濾波前 SNR = %f [dB] \n',snr1);
fprintf('濾波後 SNR = %f [dB] \n',snr2);
fprintf('提升 %f [dB] \n',snr);
fprintf('\n');

%% 聽效果
% sound(s,fs); % 乾淨的語音
% sound(r1,fs); % 含噪的語音
% sound(output,fs); % 濾波後的語音

%% 【降噪階段】

% 假設這是後來錄的一段音，混入了性質類似的噪聲。
% 現在想用剛剛得到的濾波器係數，來去除掉現在這段含噪訊號中的噪聲。

filedir=[];                             % 設定路徑
filename='risesun.wav';                % 設定檔名
fle=[filedir filename];                 % 構成完整的路徑和檔名
[s, fs] = audioread(fle);               % 讀入資料檔案

s=s-mean(s);                            % 消除直流分量
s=s/max(abs(s));                        % 幅值歸一
N=length(s);                            % 語音長度
time=(0:N-1)/fs;                        % 設定時間刻度

%% 產生噪聲並加入
% SNR = 0;                                  % 設定信噪比
r2=randn(size(s));                      % 產生隨機噪聲
b=fir1(31,0.5);                         % 設計FIR濾波器,代替H
r21=filter(b,1,r2);                     % FIR濾波
% [signal_with_noise,noise]=add_noisedata(s,data,fs,fs1,snr);
[r1,r22]=add_noisedata(s,r21,fs,fs,SNR);% 產生帶噪語音，信噪比為SNR　

%% 用上一階段求得的h，降噪

% 濾波
y = filter(h,1,r1);
output = y;

snr1=SNR_singlech(s,r1);                % 計算初始信噪比
snr2=SNR_singlech(s,output);            % 計算濾波後的信噪比
snr=snr2-snr1;

%% 畫圖
figure;
subplot 311; plot(time,s,'k'); ylabel('幅值');
ylim([-1 1 ]); title('原始語音訊號');
subplot 312; plot(time,r1,'k'); ylabel('幅值') 
ylim([-1 1 ]); title('帶噪語音訊號');
subplot 313; plot(time,output,'k'); 
ylim([-1 1 ]); title('濾波輸出語音訊號');
xlabel('時間/s'); ylabel('幅值')

%% 列印SNR
fprintf('去噪 \n',snr);
fprintf('濾波前 SNR = %f [dB] \n',snr1);
fprintf('濾波後 SNR = %f [dB] \n',snr2);
fprintf('提升 %f [dB] \n',snr);

%% 聽效果
% sound(s,fs); % 乾淨的語音
% sound(r1,fs); % 含噪的語音
sound(output,fs); % 濾波後的語音
 

維納濾波器設計的函式。

function h = my_weiner_filter1( h_length,desired_signal,observed_signal )
%function h = my_weiner_filter1( h_length,desired_signal,observed_signal )
%   維納濾波器的實現
%   輸入引數
%       h_length: 返回的FIR濾波器的長度
%       desired_signal: 所期望的訊號,訓練訊號,乾淨訊號
%       observed_signal: 觀測到的訊號
%   返回引數
%       h: FIR濾波器係數
%
%   作者：qcy
%   版本：v1.0
%   時間：2016年10月29日18:50:56


% 0. 定義線性方程組的大小
row_number = h_length;
col_number = row_number;

% 1. Rx --> observed_signal
M = col_number;
% lags = -(N-1):(N-1);
Rx_c_full = xcorr(observed_signal);
[~,k] = max(Rx_c_full);
Rx_c = Rx_c_full(k:k+M-1);
Rx_c = Rx_c.';

% 2. Rdx
Rdx_c_full = xcorr(desired_signal,observed_signal);
Rdx_c = Rdx_c_full(k:k+M-1);

% 3. 求h, Ax = b
% (1) 生成自相關A陣
A = toeplitz(Rx_c,Rx_c);
b = Rdx_c;
h = A\b;

end
 

其中，加入噪聲的程式碼如下。

function [signal,noise]=add_noisedata(s,data,fs,fs1,snr)
s=s(:);                          % 把訊號轉換成列資料
s=s-mean(s);                     % 消除直流分量
sL=length(s);                    % 求出訊號的長度

if fs~=fs1                       % 若純語音訊號的取樣頻率與噪聲的取樣頻率不相等
    x=resample(data,fs,fs1);     % 對噪聲重取樣，使噪聲取樣頻率與純語音訊號的取樣頻率相同
else
    x=data;
end

x=x(:);                          % 把噪聲資料轉換成列資料
x=x-mean(x);                     % 消除直流分量
xL=length(x);                    % 求噪聲資料長度
if xL>=sL                        % 如果噪聲資料長度與訊號資料長度不等，把噪聲資料截斷或補足
    x=x(1:sL);
else
    disp('Warning: 噪聲資料短於訊號資料，以補0來補足！')
    x=[x; zeros(sL-xL,1)];
end

Sr=snr;
Es=sum(x.*x);                    % 求出訊號的能量
Ev=sum(s.*s);                    % 求出噪聲的能量
a=sqrt(Ev/Es/(10^(Sr/10)));      % 計算出噪聲的比例因子
noise=a*x;                       % 調整噪聲的幅值
signal=s+noise;                  % 構成帶噪語音

效果示意圖。

訓練階段。

濾波階段。

求解方程得到的濾波器的頻率響應。

看得出來，因為是語音嘛，所以，設計出來的濾波器的性質可能更多的還是低通…

總之，這是一種在整個頻率範圍內的，均方意義下的最優……

最後，因為維納濾波器的前提是，訊號與噪聲不相關，訊號與噪聲寬平穩，……

反正約束條件還是挺多的，所以，雖然對訊號的噪聲抑制效果，比經典濾波器的要好，

但是，為了得到更好的效果，還需要利用更先進的技術。比如，之後要講的，自適應濾波器（adaptive filter）。

謝謝觀賞。

敬請期待。