www.pudn.com > aec_回音消除.rar > corr2292.m
% 噪声系列1
%%%%%x = [3,11,7,0,-1,4,2];
%%%%%nx =[-3:3];
load r.txt;
load r1.txt;
load echo.txt;
%pro(1:length(echo)) = round(1000*echo(1:length(echo)));
pro(1:length(echo)) = 1*echo(1:length(echo));
framelength = 160; %%%%%每帧（20ms）共160点
number1 = 123; %%%%%总帧数
number2 = 10; %%%%%时延帧数，或搜索帧数范围(number2,number3-k);number2:最大时延帧数
number3 = 6; %%%%%number3:最小时延帧数
k = 2; %%%%%求相关系数的帧数
per1 = 1; %%%%%从echo中取初始帧的移动的最小单位= 1/per1帧
per2 = 4; %%%%%从r中搜索最大相关系数的最小单位= 1/per2帧
a = 0.8;
b = 0.5;%%%%%b = 0.22;
c = 0.3;%%%%%c = 0.18;
attenuation_powerfactor = (27/17); %%%%%((number2-(number2-k))/k) --> 1/{[(2/3)^2+1+(2/3)^2)]/3]}=27/17；(9/10.5)^2 = 0.74
attenuation_absfactor = 0.5; %(9/7); %%%%% 1/[(2/3+1+2/3)/3] = 9/7
xgxs_1 = zeros(2,2);
xgxs = zeros(number1*per1,number2*per2);
opt = zeros(2,number1*per2);
opt(1,1) = 0;
opt(1,2) = 0;
aver_power_x2 = zeros(1,number2*per2);
aver_power_y2 = zeros(1,number1*per1-number2*per1);
x_pre2f = zeros(1,(number2-number3+k)*framelength);
effectcount = 0;
count = 0;
far_power_1 = 2;
far_power_2 = 4;
percent_factor = 1/2;
ratio_factor = 0.1;
for i=((number2)*per1):(number1*per1)
y = echo((i)*framelength/per1+1:(i)*framelength/per1+k*framelength);
ny = [1:k*framelength];
y1 = abs(y);
aver_power_y1(i) = mean(y1);
y2 = (y.^2);
aver_power_y2(i) = mean(y2);
x_pre2f = r((i)*framelength/per1+1-number2*framelength:(i)*framelength/per1+k*framelength-number3*framelength);
x1_pre2f = abs(x_pre2f);
aver_power_x1(i) = mean(x1_pre2f);
%for (ii=1:(number2-number3+k))
% count(ii) =
% (max(size(find(r((i)*framelength/per1+1-number2*framelength+(ii-1)*framelength:(i)*framelength/per1+1-number2*framelength+ii*framelength) > far_power))));
%end
%effectcount = max(size(find(x1_pre2f > far_power_1)));
x1_pre2f_2 = x1_pre2f(x1_pre2f > far_power_1);
%count = max(size(find((x_pre2f > percent_factor * max(x_pre2f))&(x_pre2f > far_power_1))));
%count = max(size(find(x_pre2f > far_power_2)));
x2_pre2f = (x_pre2f.^2);
aver_power_x2(i) = mean(x2_pre2f);
%if ((aver_power_x1(i) < max(aver_power_x1(i)))&(aver_power_y2(i) < (aver_power_x2(i) * attenuation_powerfactor)))
%if ((max(count) > 0.08 * framelength )&(aver_power_y2(i) < (aver_power_x2(i) * attenuation_powerfactor)))
%if ((aver_power_y2(i) < (aver_power_x2(i) * attenuation_powerfactor)))
%if (count > ratio_factor * effectcount)
%if ((count > ratio_factor * effectcount)&(aver_power_y2(i) < (aver_power_x2(i) * attenuation_powerfactor)))
if (((sum(x1_pre2f_2)/max(size(x1_pre2f_2))) > ratio_factor * max(x1_pre2f_2) )&(aver_power_y2(i) < (aver_power_x2(i) * attenuation_powerfactor)))
pro((i-1)*framelength/per1+1:(i-1)*framelength/per1+k*framelength) = zeros(1,k*framelength);
else
%y_max(i) = max(y);
for j= 1 :number2*per2
x = r((i)*framelength/per1+1-number2*framelength+j*framelength/per2:(i)*framelength/per1-number2*framelength+j*framelength/per2+k*framelength);
nx = [1:k*framelength];
%%%%%if (max(y)>0.3)
%x1 = (x.^2);
%aver_power_x(j-(i - number2)*per2) = mean(x1);
%x_max(j-(i - number2)*per2) = max(x);
if (max(y)>0.3)
if ((max(x)==0)|(max(y)==0))
xgxs(i,j) = 0;
else
xgxs_1 = corrcoef(x,y);
xgxs(i,j) = xgxs_1(2,1);
end
else
xgxs(i,j) = 0;
end
end
opt(1,i) = max(xgxs(i,:));
opt(2,i) = i - rem(min(find(xgxs(i,:)==max(xgxs(i,:)))),number1*per2);
%x1 = (x.^2)';
%y1 = (y.^2)';
%if ((opt(1,i) > 0.08) & (cumsum(y1) < cumsum(x1)))
%if ((opt(1,i) > a)|(aver_power_y(i) < (max(aver_power_x) * attenuation_factor)))
%if ((opt(1,i) > a))
%if ((opt(1,i) > a)|(y_max(i) < (max(x_max) * attenuation_factor)))
if ((opt(1,i) > a))
%pro((i-1)*framelength/per1+1:(i-1)*framelength/per1+k*framelength) = zeros(1,k*framelength);
%pro((i-1)*framelength/per1+1:(i-1)*framelength/per1+framelength) = (pro((i-1)*framelength/per1+1:(i-1)*framelength/per1+framelength)).^(0.3);
%for k =1:framelength
% if (pro((i-1)*framelength/per1+k) > mean(pro((i-1)*framelength/per1+1:(i-1)*framelength/per1+framelength)))
% pro((i-1)*framelength/per1+k) = mean(pro((i-1)*framelength/per1+1:(i-1)*framelength/per1+framelength));
% end
%end
%pro((i-1)*framelength/per1+1:(i-1)*framelength/per1+framelength) = (1000*(pro((i-1)*framelength/per1+1:(i-1)*framelength/per1+framelength))).^(0.3)/100;
%pro((i-1)*framelength/per1+1:(i-1)*framelength/per1+framelength) = zeros(1,framelength);
%%%%%else if ((opt(1,i) > b) & (opt(1,i-1) > b) & (opt(1,i-2) > b) & (i > 3))
else if ((opt(1,i) > b) & (opt(1,i-1) > b))
%pro((i-1)*framelength/per1+1:(i-1)*framelength/per1+framelength) = pro((i-1)*framelength/per1+1:(i-1)*framelength/per1+framelength)/10;
pro((i-1)*framelength/per1+1:(i-1)*framelength/per1+k*framelength) = zeros(1,k*framelength);
else if ((opt(1,i) > c) & (opt(1,i-1) > c) & (opt(1,i-2) > c))
%pro((i-1)*framelength/per1+1:(i-1)*framelength/per1+framelength) = pro((i-1)*framelength/per1+1:(i-1)*framelength/per1+framelength)/10;
%pro((i-1)*framelength/per1+1:(i-1)*framelength/per1+k*framelength) = zeros(1,k*framelength);
end
end
end
end
end
figure(3);
%subplot(411),plot(r),title('r');axis([0, 40000, -10, 10]);
%subplot(412),plot(r1),title('r1');axis([0, 40000, -10, 10]);
%subplot(413),plot(echo),title('echo');axis([0, 40000, -10, 10]);
%subplot(414),plot(pro),title('filter');axis([0, 40000,-10, 10]);
subplot(411),plot(r),title('r');
subplot(412),plot(r1),title('r1');
subplot(413),plot(echo),title('echo');
subplot(414),plot(pro),title('filter');
filter = 1000*pro';
auwrite(filter,8000,16,'linear','filter');
%auwrite(echo,8000,16,'linear','echo');
%format long;
save filter.txt filter -ascii -double;
save pro.txt pro -ascii -double;
%save filter.txt filter -ascii;
%save filter.txt filter;
%for i=1:length(pro)
% filter(2*i-1) = rem(pro(i),256);
% filter(2*i) = fix(pro(i)/256);
%end
%save filter.dat filter -ascii;
%opt(2,i) = fix(find(xgxs(i,:)==opt(1,i)),number1*per);
%opt(2,i) = fix(find(xgxs(i,:)==opt(1,i))/(number1*per))-i;
%i -rem(find(xgxs(i,:)==max(xgxs(i,:))),number1*per);
%opt(2,8) = 8 - rem(find(xgxs(8,:)==opt(1,8)),number1*per);
%opt(1,8) = max(xgxs(8,:));
% opt(2,8) = rem(find(xgxs(8,:)==opt(1,8)),number1*per);
% opt(1,16) = max(xgxs(16,:));
% opt(2,16) = rem(find(xgxs(16,:)==opt(1,16)),number1*per);
% opt(1,25) = max(xgxs(25,:));
% opt(2,25) = rem(find(xgxs(25,:)==opt(1,25)),number1*per);
%%%%%[y,ny] = sigshift(x, nx, 2);
%%%%%w = randn(1, length(y));
%%%%%nw = ny;
%%%%%[y, ny] = sigadd(y, ny, w, nw);