%*****参考文献：EMD-Based Filtering (EMDF) of Low-Frequency Noise for Speech
clear all;
close all;
clc;
fs=16000;
Frame_Len=512;
Shift_Len=128;
fid = fopen ('SI923.WAV', 'r');% Open file for reading
if fid ==- 1
errordlg ( 'Specify _valid file name!', 'ERROR: Invalid file name');
else
fseek (fid, 1024, -1);
Input_data0 = fread (fid, inf, 'int16');
fclose (fid);
end;
figure;
plot(Input_data0,'y');
title('input Data ') ;
Input_data=Input_data0-mean(Input_data0);
inputnoise=100*(rand(length(Input_data0),1)-0.5);
Input_data=Input_data+ inputnoise;
% figure;
% plot(inputnoise.^2);
figure;
plot(Input_data);
title('Data of Adding Noise') ;
% d0=round(Input_data);
% d_t=int16(d0);
% wavplay(d_t,fs);
% wavwrite(d_t,fs,16,'includenoise.wav');
%
%********************传统信噪比 ***************
% SNR=10*log10((sum(Input_data0.^2)/sum(Input_data.^2)))
%***********************************************
Input_data=Input_data';
DataLen=length(Input_data); %输入信号长度
Frame_Num=fix(DataLen/Shift_Len); %计算需要处理的帧数
Win_hm=hamming(Frame_Len)'; %汉明窗
%*******************P_speech_absence****************************
%*********参考文献：Noise Spectrum Estimation in Adverse Environments:Improved Minima
%Controlled Recursive Averaging********
%定义部分
% guji_noise=[];
Noise_power(1:Frame_Len)=0; % 噪声估计
Sf_first=zeros(1,Frame_Len); % Sf
Smin_first=zeros(1,Frame_Len); %第一次更新最小功率频谱值Smin
S_first=zeros(1,Frame_Len); %第一次回归平滑 S
Smin_second=zeros(1,Frame_Len); %第二次更新最小功率谱 Smin^
S_second=zeros(1,Frame_Len); %第二次回归平滑 S^
Sf_second=zeros(1,Frame_Len); %Sf^
I=zeros(1,Frame_Len); %初步估计语音存在概率 I(k,l)
D_window_first=1e20*ones(120, Frame_Len); % 求S最小值搜索窗 D
D_window_second=1e20*ones(120, Frame_Len);% 求S^最小值搜索窗 D
Frame_Data(1:Frame_Len)=Input_data(1:Frame_Len).*Win_hm; %第一帧
Fre_Data(1:Frame_Len)=fft(Frame_Data,Frame_Len); %变换到频域
Data_Amp=abs(Fre_Data);
Noise_power=Data_Amp.^2; %初始化Noise_power
%******************用前10帧估算噪声，计算输入分段信噪比*************
Noise_power0=0;
SEGSNR_in=0;
SEGSNR_out=0;
for j=1:5
Frame_Data(1:Frame_Len)=Input_data((j-1)*Shift_Len+1:(j-1)*Shift_Len+Frame_Len);
Noise_power0=Noise_power0+ sum(Frame_Data.^2);
end
Noise_power0= Noise_power0/5;
%******************************************************************
% guji_noise=[guji_noise,Noise_power];
% figure;
% plot( Noise_power);
Sf_first=Data_Amp(1)^2; %计算Sf(k,0)
for k=2:Frame_Len-1
Sf_first(k)=0.25*Data_Amp(k-1)^2+0.5*Data_Amp(k)^2+0.25*Data_Amp(k+1)^2;
end
Sf_first(Frame_Len)=Data_Amp(Frame_Len)^2;
%初始化
Smin_first=Sf_first; % Smin(k,0) = Sf(k,0)
S_first=Sf_first; % S(k,0)=Sf(k,0)
Smin_second=Sf_first; %Smin(k,0)^=Sf(k,0)
S_second=Sf_first; %S(k,0)^=Sf(k,0)
Sf_second=Sf_first; %Sf(k,0)^=Sf(k,0)
D_window_first(1,:)=Sf_first; %搜索窗内第一个值 为Sf(k,0)
D_window_second(1,:)=Sf_first;
D_count_first=1; %搜索窗内帧数统计
D_count_second=1;
Out_data(1:Frame_Num*Shift_Len)=0; %增强后语音赋初值
% Prio_snr=0.92; %附先验信噪比初值
Input_snr=0; %输入信噪比
Output_snr=0; %输出信噪比
Post_snr_iter=1.0;
Gh1=1.0;
dlmwrite('guji.txt',Noise_power,'precision','%10.0f');
for n=1:1:Frame_Num-3
Frame_Data(1:Frame_Len)=Input_data((n-1)*Shift_Len+1:(n-1)*Shift_Len+Frame_Len).*Win_hm; %进行分帧加窗
Fre_Data(1:Frame_Len)=fft(Frame_Data,Frame_Len); %变换到频域
Data_Amp=abs(Fre_Data); %计算幅值
Data_Ang=angle(Fre_Data); %计算相位
Frame_Data1(1:Frame_Len) = Input_data0((n-1)*Shift_Len+1:(n-1)*Shift_Len+Frame_Len)'.*Win_hm; %纯净语音分帧，计算分段信噪比
Sf_first=Data_Amp(1)^2;
for k=2:Frame_Len-1
Sf_first(k)=0.25*Data_Amp(k-1)^2+0.5*Data_Amp(k)^2+0.25*Data_Amp(k+1)^2; %计算Sf b=[0.25,0.5,0.25]
end
Sf_first(Frame_Len)=Data_Amp(Frame_Len)^2;
if n>=2
S_first=0.9*S_first+0.1*Sf_first; %S(k,l)平滑迭代
D_count_first=D_count_first+1;
if D_count_first>120 %对搜索窗内的数值更新
if rem(D_count_first,120)==0
D_window_first(120,:) = S_first;
else
D_window_first(rem(D_count_first,120),:) = S_first;
end
else
D_window_first(n,:)=S_first;
end
for k=1:Frame_Len %求第k频点 此帧前面D帧的最小值
Smin_first(k)=min( D_window_first(:,k));
end
end
Rmin_first=(Data_Amp.^2)./(1.66.*Smin_first);
Bmin_first=S_first./(1.66.*Smin_first);
for k=1:Frame_Len %初步估计语音存在的概率
if (Rmin_first(k)<4.6)&(Bmin_first(k)<1.67)
I(k)=1;
else
I(k)=0;
end
end
if n>=2
for k=2:Frame_Len-1 %计算Sf^
if (I(k-1)+I(k)+I(k+1))~=0
Sf_second(k)=(0.25*(Data_Amp(k-1)^2)*I(k-1)+0.5*Data_Amp(k)^2*I(k)+0.25*(Data_Amp(k+1)^2)*I(k+1))/(0.25*I(k-1)+0.5*I(k)+0.25*I(k+1));
end
end
S_second=0.9*S_second+0.1*Sf_second; %迭代更新S^
D_count_second=D_count_second+1;
if D_count_second>120
if rem(D_count_second,120)==0
D_window_second(120,:) = S_second;
else
D_window_second(rem(D_count_second,120),:) = S_second;
end
else
D_window_second(n,:)=S_second;
end
for k=1:Frame_Len %计算Smin^
Smin_second(k)=min(D_window_second(:,k));
end
end
Rmin_second=(Data_Amp.^2)./(1.66*Smin_second);
Bmin_second=S_first./(1.66*Smin_second);
P_speech_absence=zeros(1,Frame_Len); %求语音不存在的概率 q(k,l)
for k=1:Frame_Len
if ( Rmin_second(k)<=1)&(Bmin_second(k)<1.67)
P_speech_absence(k)=1;
elseif (1<Rmin_second(k))&(Rmin_second(k)<3)&(Bmin_second(k)<1.67)
P_speech_absence(k)=(3-Rmin_second(k))/2;
else
P_speech_absence(k)=0;
end
end
%*************************************
Post_snr=Data_Amp.^2./Noise_power; %计算后验信噪比
Prio_snr=0.92*(Gh1.^2).* Post_snr_iter +0.08*max(Post_snr-1,0);
Post_snr_iter= Post_snr;
% Prio_snr_iter=(0.95.*Prio_snr+0.05.*max(Post_snr-1,0))*1.0;%对先验信噪比进行迭代更新 E{X^2}/D
% Prio_snr_H1=1.25*Prio_snr_iter; % E{X^2|H1}/D q=0.2
% V_para=Prio_snr_H1.*Post_snr./(1+Prio_snr_H1); %计算MMSE幅度谱估计的V参数
V_para=Prio_snr.*Post_snr./(1+Prio_snr);
% G=( 1+( (P_speech_absence) ./(1-P_speech_absence) ).*( (1+Prio_snr_H1).*exp(- V_para) )).^(-1); %求语音存在的概率 p(k,l)
G=( 1+( (P_speech_absence) ./(1-P_speech_absence) ).*( (1+Prio_snr).*exp(- V_para) )).^(-1); %求语音存在的概率 p(k,l)
Gexp=zeros(1,Frame_Len);
for k=1:Frame_Len
d=0.01;
t=V_para(k):d:15;
Gexp(k) = sum((exp(-t)/t)*d);
end
% Gh1=Prio_snr_H1./(1+Prio_snr_H1).*exp(0.5*Gexp); %求语音存在时的 增益函数 Gh1
Gh1=Prio_snr./(1+Prio_snr).*exp(0.5*Gexp); %求语音存在时的 增益函数 Gh1
Gain=(Gh1.^G).*((0.01).^(1-G)); % Gain增益函数 Gmin=-25
% Gain=(Gh1.^G);
M=Data_Amp.*Gain;
Mn=M.*exp(i.*Data_Ang);
% Prio_snr=(M.^2)./Noise_power; %先验信噪比更新
Noise_power=Noise_power.*G+(0.85.*Noise_power+0.15.*(Data_Amp.^2)).*(1-G); %噪声估计更新 平滑参数=0.85
for k=1:Frame_Len
if Prio_snr(k)==0
Noise_power(k)=1.21*Noise_power(k); %偏差补偿取 1.47
end
end
% guji_noise=[guji_noise,Noise_power];
% guji=ifft( Noise_power,Frame_Len);
% guji_noise=real(guji).
