%% 维纳去噪，基于训练的
% 通过解方程求滤波器系数
% 作者：qcy% 版本说明：补充滤波阶段
% 版本：v1.1
% 时间：2016年10月29日18:18:21% 版本说明：训练阶段
% 版本：v1.0
% 时间：2016年10月29日16:52:37clear;
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);     % 对噪声重采样，使噪声采样频率与纯语音信号的采样频率相同
elsex=data;
endx=x(:);                          % 把噪声数据转换成列数据
x=x-mean(x);                     % 消除直流分量
xL=length(x);                    % 求噪声数据长度
if xL>=sL                        % 如果噪声数据长度与信号数据长度不等，把噪声数据截断或补足x=x(1:sL);
elsedisp('Warning: 噪声数据短于信号数据，以补0来补足！')x=[x; zeros(sL-xL,1)];
endSr=snr;
Es=sum(x.*x);                    % 求出信号的能量
Ev=sum(s.*s);                    % 求出噪声的能量
a=sqrt(Ev/Es/(10^(Sr/10)));      % 计算出噪声的比例因子
noise=a*x;                       % 调整噪声的幅值
signal=s+noise;                  % 构成带噪语音