ITUT-T recommendations G.168 标准回声模型

相关标准定义见《ITUT-T recommendations G.168 标准》附件D，代码中K值、ERL值需要按照标准文件中的值调整一下

回声模型1

clc
clear all
close all


f=0:1:4095;
i=1:8;
m1=[-436 -829 -2797 -4208 -17968 -11215 46150 34480 -10427 9049 -1309 -6320 ...
    390 -8191 -1751 -6051 -3796 -4055 -3948 -2557 -3372 -1808 -2259 -1300 ...
    -1098 -618 -340 -61 328 419 745 716 946 880 1014 976 1033 1091 1053 1042 ...
    794 831 899 716 390 313 304 304 73 -119 -109 -176 -359 -407 -512 -580 -704 ...
    -618 -685 -791 -772 -820 -839 -724];
figure
plot(m1);
title('m1');


L=length(m1);
M1=zeros(1,length(f));
size(M1);


for loop1=0:4095    
    for k=1:L
        M1(loop1+1)=M1(loop1+1)+m1(k)*exp(-2*1j*pi*loop1*k/8192);
    end
end


figure
plot(M1);
title('M1');


% K1=1/(max(abs(M1)));
K1= 1.39*10^-5;


delay=4000;
gk1=zeros(1,L);
for k=1+delay:L+delay
    gk1(k)=(10^(-6/20)*K1)*m1(k-delay);
end


figure
plot(gk1);
title('gk1');


% t=0:0.001:10;
% s=sin(2*pi*50*t);
[s1,fs]=audioread('handel.wav');
s=s1';
Es=conv(gk1,s);


d=zeros(1,length(Es));
for loop1=1:length(s)
    d(loop1)=s(loop1); 
end
d=d+Es;


figure
subplot(3,1,1);
plot(s)
subplot(3,1,2)
plot(Es)
subplot(3,1,3)
plot(d)


sound(s,8000)
pause(10)
sound(d,8000)
           

回声模型1的示意图如下：

代码中让模型产生了一个延时4000点的回声，因为handle.wav的采样率为8000，所以在时间轴上产生的延时是0.5秒，这样在播放音乐的时候可以明显的听到回声的效果。

回声模型2

clc
clear all
close all

f=0:1:4095;
i=1:8;
m2=[-381 658 1730 -51 -3511 -1418 7660 8861 -8106 -21370 -5307 23064 24020 ...
    1020 -12374 -16296 -19524 -7480 13509 17115 13952 13952 97 -9326 -9046 ...
    -15208 -9853 -3858 -1979 6029 5616 7214 6820 3935 3919 921 1316 -693 ...
    -759 -1517 -2176 -2028 -2654 -1814 -2077 -1468 -1221 -842 -463 -298 -68 ...
    64 493 723 789 954 756 839 872 1020 789 822 558 658 476 377 377 262 97 ...
    -68 -183 -232 -331 -347 -430 -314 -430 -463 -463 -414 -381 -479 -479 -512 ...
    -479 -397 -430 -397 -298 -265 -249 -216 -249 -265 -166 -232];
figure
plot(m2);
title('m2');

L=length(m2);
M2=zeros(1,length(f));
size(M2);

for loop1=0:4095    
    for k=1:L
        M2(loop1+1)=M2(loop1+1)+m2(k)*exp(-2*1j*pi*loop1*k/8192);
    end
end

figure
plot(M2);
title('M2');

% K2=1/(max(abs(M2)));
K2= 1.44*10^-5;

delay=4000;
gk2=zeros(1,L);
for k=1+delay:L+delay
    gk2(k)=(10^(-6.55/20)*K2)*m2(k-delay);
end

figure
plot(gk2(delay:delay+L));
title('gk2');

% t=0:0.001:10;
% s=sin(2*pi*50*t);
[s1,fs]=audioread('handel.wav');
s=s1';
Es=conv(gk2,s);

d=zeros(1,length(Es));
for loop1=1:length(s)
    d(loop1)=s(loop1); 
end
d=d+Es;

figure
subplot(3,1,1);
plot(s)
title('语音原声');
subplot(3,1,2)
plot(Es)
title('model2产生的回声');
subplot(3,1,3)
plot(d)
title('听到的声音');

sound(s,8000)
pause(10)
sound(d,8000)


           

回声模型2的冲击响应

回声模型3

clc
clear all
close all

f=0:1:4095;
i=1:8;
m3=[-448 -436 2230 2448 -4178 -7050 5846 18581 2322 -26261 -16249 21637 25649 ...
    -2267 -10311 -4693 -12690 -7428 14164 13467 4438 8627 456 -11879 -6352 -5104 ...
    -7496 3271 6566 4277 11131 7562 1475 3728 -3525 -7301 -3101 -9269 -6146 -2553 ...
    -6272 811 124 788 5147 2172 5387 4598 3535 4004 2311 2150 1017 330 -139 -573 -1100 ...
    -1008 -1077 -1088 -917 -917 -963 -814 -871 -734 -642 -562 -356 -379 -345 -230 ...
    -233 -333 -356 -390 -310 -265 -368 -310 -310 -390 -482 -459 -482 -551 -573];
figure
plot(m3);
title('m3');

L=length(m3);
M3=zeros(1,length(f));
size(M3);

for loop1=0:4095    
    for k=1:L
        M3(loop1+1)=M3(loop3)+m3(k)*exp(-2*1j*pi*loop1*k/8192);
    end
end

figure
plot(M3);
title('M3');

% K3=1/(max(abs(M3)));
K3= 1.52*10^-5;

delay=4000;
gk3=zeros(1,L);
for k=1+delay:L+delay
    gk3(k)=(10^(-6/20)*K3)*m3(k-delay);
end

figure
plot(gk3(delay:delay+L));
title('gk3');

% t=0:0.001:10;
% s=sin(2*pi*50*t);
[s1,fs]=audioread('handel.wav');
s=s1';
Es=conv(gk3,s);

d=zeros(1,length(Es));
for loop1=1:length(s)
    d(loop1)=s(loop1); 
end
d=d+Es;

figure
subplot(3,1,1);
plot(s)
title('语音原声');
subplot(3,1,2)
plot(Es)
title('model3产生的回声');
subplot(3,1,3)
plot(d)
title('听到的声音');

sound(s,8000)
pause(10)
sound(d,8000)


           

回声模型3的冲击响应

回声模型4

clc
clear all
close all

f=0:1:4095;
i=1:8;
m4=[160 312 -241 -415 897 908 -1326 -1499 2405 3347 -3624 -7733 4041 14484 ...
    -1477 -21739 -4470 25356 11458 -19696 -11800 5766 789 6633 14624 -6975 ...
    -17156 -187 149 1515 14907 4345 -7128 -2757 -10185 -7083 6850 3944 6969 ...
    8694 -4068 -3852 -5793 -9371 453 1060 3965 9463 2393 2784 -892 -7366 -3376 ...
    -5847 -2399 3011 1537 6623 4205 1602 1592 -4752 -3646 -5207 -5577 -501 -1174 ...
    4041 5647 4628 7252 2123 2654 -881 -4113 -3244 -7289 -3830 -4600 -2508 431 ...
    -144 4184 2372 4617 3576 2382 2839 -404 539 -1803 -1401 -1705 -2269 -783 -1608 ...
    -220 -306 257 615 225 561 8 344 127 -57 182 41 203 -111 95 -79 30 84 -13 -68 ...
    -241 -68 -24 19 -57 -24 30 -68 84 -155 -68 19];
figure
plot(m4);
title('m4');

L=length(m4);
M4=zeros(1,length(f));
size(M4);

for loop1=0:4095    
    for k=1:L
        M4(loop1+1)=M4(loop4)+m4(k)*exp(-2*1j*pi*loop1*k/8192);
    end
end

figure
plot(M4);
title('M4');

% K4=1/(max(abs(M4)));
K4= 1.77*10^-5;

delay=4000;
gk4=zeros(1,L);
for k=1+delay:L+delay
    gk4(k)=(10^(-6/20)*K4)*m4(k-delay);
end

figure
plot(gk4(delay:delay+L));
title('gk4');

% t=0:0.001:10;
% s=sin(2*pi*50*t);
[s1,fs]=audioread('handel.wav');
s=s1';
Es=conv(gk4,s);

d=zeros(1,length(Es));
for loop1=1:length(s)
    d(loop1)=s(loop1); 
end
d=d+Es;

figure
subplot(3,1,1);
plot(s)
title('语音原声');
subplot(3,1,2)
plot(Es)
title('model4产生的回声');
subplot(3,1,3)
plot(d)
title('听到的声音');

sound(s,8000)
pause(10)
sound(d,8000)



           

模型4的冲击响应

回声模型5

clc
clear all
close all

f=0:1:4095;
i=1:8;
m5=[293 268 475 460 517 704 581 879 573 896 604 787 561 538 440 97 265 -385 ...
    20 -983 -523 -1438 -1134 -1887 -1727 -1698 -4266 -22548 -43424 2743 25897 ...
    7380 21499 11983 10400 11667 3889 7241 925 2018 -821 -2068 -2236 -4283 ...
    -3406 -5022 -4039 -4842 -4104 -4089 -3582 -2978 -2734 -1805 -1608 -645 ...
    -495 279 471 947 1186 1438 1669 1640 1901 1687 1803 1543 1566 1342 1163 963 ...
    733 665 323 221 -14 -107 -279 -379 -468 -513 -473 -588 -612 -652 -616 -566 ...
    -515 -485 -404 -344 -290 -202 -180 -123];
figure
plot(m5);
title('m5');

L=length(m5);
M5=zeros(1,length(f));
size(M5);

for loop1=0:4095    
    for k=1:L
        M5(loop1+1)=M5(loop5)+m5(k)*exp(-2*1j*pi*loop1*k/8192);
    end
end

figure
plot(M5);
title('M5');

% K5=1/(max(abs(M5)));
K5= 0.933*10^-5;

delay=4000;
gk5=zeros(1,L);
for k=1+delay:L+delay
    gk5(k)=(10^(-6/20)*K5)*m5(k-delay);
end

figure
plot(gk5(delay:delay+L));
title('gk5');

% t=0:0.001:10;
% s=sin(2*pi*50*t);
[s1,fs]=audioread('handel.wav');
s=s1';
Es=conv(gk5,s);

d=zeros(1,length(Es));
for loop1=1:length(s)
    d(loop1)=s(loop1); 
end
d=d+Es;

figure
subplot(3,1,1);
plot(s)
title('语音原声');
subplot(3,1,2)
plot(Es)
title('model5产生的回声');
subplot(3,1,3)
plot(d)
title('听到的声音');

sound(s,8000)
pause(10)
sound(d,8000)


           

模型5的冲击响应

回声模型6

clc
clear all
close all

f=0:1:4095;
i=1:8;
m6=[29 109 -83 198 -249 -135 -415 -202 -444 -337 -313 -450 -105 -503 145 -490 267 -231 340 77 ...
    343 783 158 1341 195 1798 344 1845 629 1604 1182 940 5163 19522 8421 -50953 -9043 18046 ...
    -13553 13336 -3471 -107 1788 -7409 2469 -7994 490 -3860 -837 490 -636 3682 1141 5019 2635 ...
    5025 3946 4414 4026 3005 3380 1616 2007 158 388 -1198 -1117 -2134 -2547 -2589 -3310 -2778 ...
    -3427 -2779 -3116 -2502 -2399 -1956 -1539 -1239 -570 -377 251 331 964 1177 1449 1564 1724 ...
    1871 1767 1802 1630 1632 1379 1271 1063 856 711 482 289 54 -137 -321 -490 -638 -764 -836 ...
    -800 -859 -838 -837 -834 -740 -673 -581 -493 -436 -327 -201];
figure
plot(m6);
title('m6');

L=length(m6);
M6=zeros(1,length(f));
size(M6);

for loop1=0:4095    
    for k=1:L
        M6(loop1+1)=M6(loop6)+m6(k)*exp(-2*1j*pi*loop1*k/8192);
    end
end

figure
plot(M6);
title('M6');

% K6=1/(max(abs(M6)));
K6= 1.51*10^-5;

delay=4000;
gk6=zeros(1,L);
for k=1+delay:L+delay
    gk6(k)=(10^(-6/20)*K6)*m6(k-delay);
end

figure
plot(gk6(delay:delay+L));
title('gk6');

% t=0:0.001:10;
% s=sin(2*pi*50*t);
[s1,fs]=audioread('handel.wav');
s=s1';
Es=conv(gk6,s);

d=zeros(1,length(Es));
for loop1=1:length(s)
    d(loop1)=s(loop1); 
end
d=d+Es;

figure
subplot(3,1,1);
plot(s)
title('语音原声');
subplot(3,1,2)
plot(Es)
title('model6产生的回声');
subplot(3,1,3)
plot(d)
title('听到的声音');

sound(s,8000)
pause(10)
sound(d,8000)
           

模型6的冲击响应

回声模型7

clc
clear all
close all

f=0:1:4095;
i=1:8;
m7=[258 -111 337 347 -434 192 -450 -108 -343 -596 -177 -1187 -52 -1781 -147 -1959 -326 -1601 ...
    -1389 -13620 -720 33818 -10683 -6742 12489 -9862 8950 -1574 758 3526 -3118 2421 -8966 ...
    -4901 11385 18072 -14410 -7473 19836 -16854 -3115 9483 -17799 7399 -4342 -7145 7929 -10726 ...
    6239 -2526 -1317 5345 -4565 6868 -2195 3425 1969 -109 3963 -1275 3087 -892 1239 2 -427 ...
    596 -1184 551 -1244 141 -743 -415 -372 -769 -183 -785 -270 -659 -3778 -523 -325 -245 -255 ...
    -60 35 218 149 340 233 365 303 251 230 209 179 ];
figure
plot(m7);
title('m7');

L=length(m7);
M7=zeros(1,length(f));
size(M7);

for loop1=0:4095    
    for k=1:L
        M7(loop1+1)=M7(loop7)+m7(k)*exp(-2*1j*pi*loop1*k/8192);
    end
end

figure
plot(M7);
title('M7');

% K6=1/(max(abs(M7)));
K6= 2.33*10^-5;

delay=4000;
gk7=zeros(1,L);
for k=1+delay:L+delay
    gk7(k)=(10^(-11.06/20)*K6)*m7(k-delay);
end

figure
plot(gk7(delay:delay+L));
title('gk7');

% t=0:0.001:10;
% s=sin(2*pi*50*t);
[s1,fs]=audioread('handel.wav');
s=s1';
Es=conv(gk7,s);

d=zeros(1,length(Es));
for loop1=1:length(s)
    d(loop1)=s(loop1); 
end
d=d+Es;

figure
subplot(3,1,1);
plot(s)
title('语音原声');
subplot(3,1,2)
plot(Es)
title('model7产生的回声');
subplot(3,1,3)
plot(d)
title('听到的声音');

sound(s,8000)
pause(10)
sound(d,8000)


           

模型7的冲击响应

回声模型8

clc
clear all
close all

f=0:1:4095;
i=1:8;
m8=[80 31 4 42 42 -61 -81 -64 -121 -102 -26 1002 -9250 -22562 39321 35681 -35289 ...
    25312 -1457 -229 15659 -6786 16791 3860 2239 -28730 -11885 33871 -176 -16421 ...
    18173 -9669 -10163 9941 -19365 3592 -5907 -10257 5336 -12933 4348 -4802 -1791 ...
    3035 -4433 5553 -2596 3992 1255 1450 4079 324 4340 1059 3083 1917 1756 2478 ...
    1027 1871 845 1284 813 806 869 471 646 438 449 432 473 394 452 538 717 723 ...
    850 756 753 899 555 669 619 500 650 615 516 492 427 291 356 147 107 -50 -88 ...
    -59 -238 -165 -183 ];
figure
plot(m8);
title('m8');

L=length(m8);
M8=zeros(1,length(f));
size(M8);

for loop1=0:4095    
    for k=1:L
        M8(loop1+1)=M8(loop8)+m8(k)*exp(-2*1j*pi*loop1*k/8192);
    end
end

figure
plot(M8);
title('M8');

% K6=1/(max(abs(M8)));
K6= 1.33*10^-5;

delay=4000;
gk8=zeros(1,L);
for k=1+delay:L+delay
    gk8(k)=(10^(-9.27/20)*K6)*m8(k-delay);
end

figure
plot(gk8(delay:delay+L));
title('gk8');

% t=0:0.001:10;
% s=sin(2*pi*50*t);
[s1,fs]=audioread('handel.wav');
s=s1';
Es=conv(gk8,s);

d=zeros(1,length(Es));
for loop1=1:length(s)
    d(loop1)=s(loop1); 
end
d=d+Es;

figure
subplot(3,1,1);
plot(s)
title('语音原声');
subplot(3,1,2)
plot(Es)
title('model8产生的回声');
subplot(3,1,3)
plot(d)
title('听到的声音');

sound(s,8000)
pause(10)
sound(d,8000)


           

模型8的冲击响应