📄 cp3m2.m
字号:
clear all;
[y1,fs,bits]=wavread('D:\语音降噪\voise\3.wav');
%y1=y1/max(abs(y1));%语音信号归一化
wavwrite(y1,8000,8,'D:\语音降噪\wav\4.wav');
figure(1);
plot(y1);
[noise,fs1,bits1]=wavread('D:\语音降噪\voise\3_noise.wav');
y=mixsig(y1,noise,0);% 混合
%y=y/max(abs(y));%归一化
wavwrite(y,8000,8,'D:\语音降噪\wav\mymasking_s&w.wav');%0db带噪信号
figure(2);
plot(y);
frame = 256; % Defining frame size
shift=64;
win=hamming(256);
for j1 = 1:length(y),
signal(j1) = y(j1);
end;
ps_noise=zeros(length(signal)/shift-3,frame);
frame_temp = zeros(length(signal)/shift-3,frame);
framenoise_temp = zeros(length(signal)/shift-3,frame);
hh = 0;
for k = 1 : 5,
for l = 1 : frame,
b(l) = signal(hh+l);
end;
hh = hh + frame;
frame_temp(k,1:frame) = abs(fft(b)); %fft for the first 50 frames
ps_noise(k,1:frame) = (frame_temp(k,1:frame).*conj(frame_temp(k,1:frame)))/frame;
end;%语音前五真早声能量
ps_noise2(1,1:frame)= (sum(ps_noise(1:k,1:frame))/5);
head = 0;
mm=1;
nn=1;
% START OF THE NOISE ELIMINATION THROUGH SPECTRAL SUBTRACTION BASED ON THE THRESHOLD SET
% for k = 1 : length(signal)/shift-3,
% for m = 1 : frame,
% abc1(m) = signal(head+m);
% end;
% abc1=abc1.*win';
%
% head = head +shift;
% frame_temp(k,1:frame) = abs(fft(abc1));% FFT OF THE SIGNAL + NOISE FRAME BY FRAME
%
% frame_angle(k,1:frame) = angle(fft(abc1));% ANGLE OF FFT OF THE SIGNAL + NOISE FRAME BY FRAME
%
% ps_signal(k,1:frame) = (frame_temp(k,1:frame).*conj(frame_temp(k,1:frame)))./frame;
%
% %ps_temp=zeros(1,frame);
% ps_temp(k,1:frame)=ps_signal(k,1:frame);
%
% if k==1
% ps_noise(k,1:frame)=ps_noise(1,1:frame);
% % ps_signal(k,1:frame)=0.98*ps_signal(k,1:frame)+0.02*ps_signal(k,1:frame);
% else
% ps_noise(k,1:frame)=ps_noise(1,1:frame);%0.99*ps_noise(1,1:frame)+0.01*ps_signal(k,1:frame);
% ps_signal(k,1:frame)=ps_signal(k,1:frame);%0.99*ps_signal(k,1:frame)+0.01*ps_noise(1,1:frame);%+0.01*ps_signal(k-1,1:frame);
% end
% ps_noise(k,1:frame)=ps_noise(1,1:frame);
% frame_ps(1,k) = sum(ps_signal(k,1:frame));%一阵内信号能量之和
% frame_pn(1,k)=sum(ps_noise(k,1:frame));
% ps_final(1,k) = frame_ps(1,k)- 0.8*frame_pn(1,k); % 0.8*
%
% aa=1;bb=1;
% h(1,k)=power(ps_final(1,k)/(ps_final(1,k)+aa*frame_pn(1,k)),bb);
% frame1(k,1:frame) = h(1,k).*(frame_temp(k,1:frame));%为那滤波
%
% ps_signal(k,1:frame) = (frame1(k,1:frame).*conj(frame1(k,1:frame)))./frame;
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%
% T=zeros(1,256);
% b=zeros(1,22);c=zeros(1,22);o=zeros(1,22);
% sf=zeros(22,22);
%
%
%
% for i1=1:3,
% b(1)=b(1)+ps_signal(k,i1);
% end;
% for i1=4:6,
% b(2)=b(2)+ps_signal(k,i1);
% end;
% for i1=7:10,
% b(3)=b(3)+ps_signal(k,i1);
% end;
% for i1=11:13,
% b(4)=b(4)+ps_signal(k,i1);
% end;
% for i1=14:16,
% b(5)=b(5)+ps_signal(k,i1);
% end;
% for i1=17:20,
% b(6)=b(6)+ps_signal(k,i1);
% end;
% for i1=21:25,
% b(7)=b(7)+ps_signal(k,i1);
% end;
% for i1=26:29,
% b(8)=b(8)+ps_signal(k,i1);
% end;
% for i1=30:35,
% b(9)=b(9)+ps_signal(k,i1);
% end;
% for i1=36:41,
% b(10)=b(10)+ps_signal(k,i1);
% end;
% for i1=42:47,
% b(11)=b(11)+ps_signal(k,i1);
% end;
% for i1=48:55,
% b(12)=b(12)+ps_signal(k,i1);
% end;
% for i1=56:64,
% b(13)=b(13)+ps_signal(k,i1);
% end;
% for i1=65:74,
% b(14)=b(14)+ps_signal(k,i1);
% end;
% for i1=75:86,
% b(15)=b(15)+ps_signal(k,i1);
% end;
% for i1=87:101,
% b(16)=b(16)+ps_signal(k,i1);
% end;
% for i1=102:118,
% b(17)=b(17)+ps_signal(k,i1);
% end;
% for i1=119:141,
% b(18)=b(18)+ps_signal(k,i1);
% end;
% for i1=142:170,
% b(19)=b(19)+ps_signal(k,i1);
% end;
% for i1=171:205,
% b(20)=b(20)+ps_signal(k,i1);
% end;
% for i1=206:246,
% b(21)=b(21)+ps_signal(k,i1);
% end;
% for i1=247:256,
% b(22)=b(22)+ps_signal(k,i1);
% end;
%
%
% % for i1=1:22,
% % for j1=1:22,
% % sf(i1,j1)=15.81+7.5*((i1-j1)+0.474)-17.5*sqrt(1+((i1-j1)+0.474)*((i1-j1)+0.474)); %$ (18)
% % sf(i1,j1)=power(10,sf(i1,j1)/(20)); %$ transform from dB
% % end;
% % end;
%
% for i1=1:22,
% sf(i1)=15.81+7.5*(i1+0.474)-17.5*sqrt(1+(i1+0.474)*(i1+0.474));
% end;
%
% % for j1=1:22,
% % c(j1)=0;
% % for i1=1:22,
% % c(j1)=c(j1)+b(i1)*sf(i1,j1); %$ (19)
% % end;
% % end;
%
% cc_temp=conv2(b,sf);
% for i1=1:22
% for j1=1:22
% c(i1)=sum(cc_temp(i1,j1));
% end;
% end;
%
% temp_value=0.0;
% for i1=1:22,
% temp_value=temp_value+b(i1);
% end;
% ua=temp_value/256.0;
%
% temp_value=0;
% for i1=1:256,
% temp_value=temp_value+log10(ps_signal(k,i1));
% end;
% temp_value=temp_value/256.0;
% uj=power(10,temp_value);
%
% sfm=-10*log10(uj/ua); %$ (20)
%
% u=min(sfm/(-60),1); %$ (21)
%
% for i1=1:22,
% O(i1)=u*(14.5+i)+(1-u)*5.5;
% T(i1)=power(10,log10(c(i1))-O(i1)/10);
% c(i1)=T(i1);%c[i]暂存T[i]
% end;
% for i1=1:3,
% T(i1)=c(1);
% end;
% for i1=4:6,
% T(i1)=c(2);
% end;
% for i1=7:10,
% T(i1)=c(3);
% end;
% for i1=11:13,
% T(i1)=c(4);
% end;
% for i1=14:16,
% T(i1)=c(5);
% end;
% for i1=17:20,
% T(i1)=c(6);
% end;
% for i1=21:25,
% T(i1)=c(7);
% end;
% for i1=26:29,
% T(i1)=c(8);
% end;
% for i1=30:35,
% T(i1)=c(9);
% end;
% for i1=36:41,
% T(i1)=c(10);
% end;
% for i1=42:47,
% T(i1)=c(11);
% end;
% for i1=48:55,
% T(i1)=c(12);
% end;
% for i1=56:64,
% T(i1)=c(13);
% end;
% for i1=65:74,
% T(i1)=c(14);
% end;
% for i1=75:86,
% T(i1)=c(15);
% end;
% for i1=87:1010,
% T(i1)=c(16);
% end;
% for i1=102:118,
% T(i1)=c(17);
% end;
% for i1=119:141,
% T(i1)=c(18);
% end;
% for i1=142:170,
% T(i1)=c(19);
% end;
% for i1=171:205,
% T(i1)=c(20);
% end;
% for i1=206:246
% T(i1)=c(21);
% end;
% for i1=247:256,
% T(i1)=c(22);
% end;
%
%
% %////计算绝对听阈////
% mm=0.0;
% for i1=1:256,
% f(i1)=mm;
% mm=mm+8/256;
% end;
%
% f(1)=f(2);
% for i1=1:256,
% %f(i1)=(3.64*power(f(i1),-0.8)-6.5*exp(-0.6*power((f(i1)-3.3),2))+0.001*power(f(i1),4));
% f(i1)=3.64*(f(i1).^(-0.8))- 6.5*exp(-0.6*((f(i1)-3.3).^2))+0.001*(f(i1).^4);
% end;
%
% for i1=1:256,
% T(i1)=max(T(i1),f(i1));
% end;
% tmax=0;
% tmin=0;
% for i1=1:256,
% if(T(i1)>tmax)
% tmax=T(i1);
% end;
% if(T(i1)<tmin)
% tmin=T(i1);
% end;
% end;
%
%
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%
% for k1=1:256,
%
% alfa(k1)=(tmax*6-T(k1)*6+1*T(k1)-1*tmin)/(tmax-tmin);
% beita(k1)=(tmax*0.02-T(k1)*0.02+0*T(k1)-0*tmin)/(tmax-tmin);
%
%
% temppow=abs(ps_noise(k,k1))/abs(ps_temp(k,k1));
% if(power(temppow,2)<(1.0/(alfa(k1)+beita(k1))))
% %if(1-alfa[k]*pow(temppow,2)>0)
% frame1( k,k1)=power((1-alfa(k1)*power(temppow,2)),0.5)*ps_temp(k,k1);
% else
% frame1(k,k1)=power((beita(k1)*power(temppow,2)),0.5)*ps_temp(k,k1);
% end;
% end;
% end;
% hh = 0;
% for k = 1 : 5,
% ps_noise(k,1:frame) = (frame_temp(k,1:frame).*conj(frame_temp(k,1:frame)))/frame;
% end;
% %语音前五真早声能量
% ps_noise2(1,1:frame)= (sum(ps_noise(1:5,1:frame))/5);
% % frame1_initial(1:frame)=sum(frame1(1:5,1:frame))/5;
for k = 1 : length(signal)/shift-3,
for m = 1 : frame,
abc1(m) = signal(head+m);
end;
abc1=abc1.*win';
head = head +shift;
frame_temp(k,1:frame) = abs(fft(abc1));% FFT OF THE SIGNAL + NOISE FRAME BY FRAME
frame_angle(k,1:frame) = angle(fft(abc1));% ANGLE OF FFT OF THE SIGNAL + NOISE FRAME BY FRAME
ps_noise(k,1:frame)=ps_noise2(1,1:frame);
ps_signal(k,1:frame) = (frame_temp(k,1:frame).*conj(frame_temp(k,1:frame)))./frame;
frame_ps(1,k) = (sum(ps_signal(k,1:frame)));%一阵内信号能量之和
frame_pn(1,k)=sum(ps_noise(k,1:frame));
ps_final(1,k) = frame_ps(1,k)- 0.8*frame_pn(1,k);
if ps_final(1,k)<0
ps_final(1,k)=0;
end
aa=1;bb=1;
h(1,k)=power(ps_final(1,k)/(ps_final(1,k)+aa*frame_pn(1,k)),bb);
frame1(k,1:frame) = h(1,k).*(frame_temp(k,1:frame));%为那滤波
frame2(k,1:frame)=ifft(frame1(k,1:frame));
%signal(1,(((k-1)*frame)+1):(k*frame)) = frame2(k,1:frame); % Retriving back the signal(after spectral subtraction)
if k==1
signal(1,1:shift)=frame2(k,1:shift);
else if k==2
signal(1,(shift+1):(2*shift))=(frame2(k,1:shift)+frame2(k-1,(shift+1):(shift*2)))/2;
else if k==3
signal(1,(shift*2+1):(3*shift))=(frame2(k,1:shift)+frame2(k-1,(shift+1):(shift*2))+frame2(k-2,(shift*2+1):(shift*3)))/3;
%else if k==length(signal)/shift-3
%signal(1,(k*shift+1):((k+1)*shift))=frame2(k,(3*shift+1):frame);
else signal(1,(((k-1)*shift+1):(k*shift)))=(frame2(k,1:shift)+frame2(k-1,(shift+1):(shift*2))+frame2(k-2,(shift*2+1):(shift*3))+frame2(k-3,(shift*3+1):frame))/4;
%end;
end;
end;
end;
% frame1(k,1:frame) = frame1(k,1:frame).*(exp(i*frame_angle(k,1:frame)));
% frame2(k,1:frame)=ifft(frame1(k,1:frame));
% signal(1,(((k-1)*frame)+1):(k*frame)) = frame2(k,1:frame); % Retriving back the signal(after spectral subtraction)
end
signal((length(y)-1000):(length(y)))=[]; %give up 4 frames in the end
y1((length(y1)-1000):(length(y1)))=[];
figure(3);
plot(1:length(frame_ps),frame_ps,1:length(ps_final),ps_final);
figure(4);
signal=signal';
overall_snr = 10*log10(sum(abs(y1).^2)/sum((abs(y1-signal)).^2))
signal=signal/max(abs(signal));%归一化
plot(1:length(signal),signal);
%以下画语谱图
%map=(log10(1:0.1:10))';
%map=[map map map ];
%subplot(2,2,1)
%spgrambw( noise,8000);title('noise Specgram');colormap(1-map);
%subplot(2,2,2)
%spgrambw(y1,8000);title('sig Specgram');colormap(1-map);
%subplot(2,2,3)
%spgrambw(y,8000);title('mixsig Specgram');colormap(1-map);
%subplot(2,2,4)
%spgrambw(signal,8000);title('Proposed Algorithm Specgram');colormap(1-map);
%spgrambw(sig,8000);title('pure Specgram');colormap(1-map);
%[overall_snr2,seg_nr2]=snr(y1,signal)
%overall_snr = 10*log10(sum(abs(y1).^2)/sum((abs(y1-signal)).^2))
wavwrite(signal,8000,8,'c\wav\mymasking21(0db).wav');⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -