dsbf.m

delay and sum beamforming algorithm

% clear all
clear all;
close all;
clc;

% extracted sources' signals
load alz.mat;
load fw.mat;

m=min(length(alz),length(fw));
N1=1;%20000;
N2=m;%N1+3000;
s1=alz(N1:N2);  % Adjust the power of alz
s2=fw(N1:N2);
s1=(s1-min(s1))/(max(s1)-min(s1));
s2=(s2-min(s2))/(max(s2)-min(s2));
fs=16000;
N=8;
auwrite(s1,fs,N,'s1');
auwrite(s2,fs,N,'s2');
s1=1.3*s1;
Enegy1=0;
Enegy2=0;
for k=1:length(s1)
    Enegy1=Enegy1+s1(k)^2;
    Enegy2=Enegy2+s2(k)^2;
end
Ratio=Enegy1/Enegy2;

S1=fft(s1);
S2=fft(s2);

sita1=-pi/3;    % the angle of the first source.  range: from -pi/2 to pi/2
sita2=pi/3;     % the angle of the second source. range: from -pi/2 to pi/2
r1=1.0;         % the radius of the first sound source
r2=1.0;         % the radius of the second sound source
Nfft=length(S1);
D=0.433;        % distance between two microphones, suppose the diameter of circle is 0.5.
V=344.21;       %23 degrees    % velocity of sounds
fs=1.6e4;       % the sampling frequency
Nbeg=0.5;
Nend=1.5;
No=11;

%******************************************
% Variants for Reverberation
%******************************************
T60=0.4;            %  time for degrading up to 60db; seen [1]
L=0.5*T60*fs;       %  length of filter for reverberation

% Mixed sound signals in near--field equation
d11=sqrt(D*D/4+r1*r1+D*r1*sin(sita1));
d12=sqrt(3*D*D/4+r1*r1-sqrt(3)*D*r1*cos(sita1));
d13=sqrt(D*D/4+r1*r1-D*r1*sin(sita1));
d21=sqrt(D*D/4+r2*r2+D*r2*sin(sita2));
d22=sqrt(3*D*D/4+r2*r2-sqrt(3)*D*r2*cos(sita2));
d23=sqrt(D*D/4+r2*r2-D*r2*sin(sita2));
for f=1:Nfft
      X1(f)=exp(-i*2*pi*f/Nfft*d11*fs/V)*d11/d11*S1(f)+exp(-i*2*pi*f/Nfft*d21*fs/V)*d21/d21*S2(f);
      X2(f)=exp(-i*2*pi*f/Nfft*d12*fs/V)*d11/d12*S1(f)+exp(-i*2*pi*f/Nfft*d22*fs/V)*d21/d22*S2(f);
      X3(f)=exp(-i*2*pi*f/Nfft*d13*fs/V)*d11/d13*S1(f)+exp(-i*2*pi*f/Nfft*d23*fs/V)*d21/d23*S2(f);
end
mid_mix1=real(ifft(X1));
mid_mix2=real(ifft(X2));
mid_mix3=real(ifft(X3));
%*********Add reverberation***************

%*********Impulse Response[2]*************

% % addtional noise: 40dB
% save('X1','X1');
% save('X2','X2');
% save('X3','X3');

% load X1;
X1=awgn(X1,20);
% load X2;
X2=awgn(X2,20);
% load X3;
X3=awgn(X3,20);

mixed1=real(ifft(X1));
mixed2=real(ifft(X2));
mixed3=real(ifft(X3));

mixed1=(mixed1-min(mixed1))/(max(mixed1)-min(mixed1));
mixed2=(mixed2-min(mixed2))/(max(mixed2)-min(mixed2));
mixed3=(mixed3-min(mixed3))/(max(mixed3)-min(mixed3));

figure;
subplot(311);
plot(mixed1,'-r');
xlabel('signal points')
ylabel('Amplitude')
title('Mixed signal1')
subplot(312);
plot(mixed2,'-r');
xlabel('signal points')
ylabel('Amplitude')
title('Mixed signal2')
subplot(313);
plot(mixed3,'-r');
xlabel('signal points')
ylabel('Amplitude')
title('Mixed signal3')

% %DSBF
% B1=0;
% L=zeros(19,No);
% for belta=-pi/2:pi/18:pi/2
%     B1=B1+1;
%     R1=0;
%     z(B1)=0;
%     for focus=Nbeg:0.1:Nend
%         R1=R1+1;
%         Z(B1,R1)=0;
%         d1=sqrt(D*D/4+focus^2+D*focus*sin(belta));
%         d2=sqrt(3*D*D/4+focus^2-sqrt(3)*D*focus*cos(belta));
%         d3=sqrt(D*D/4+focus^2-D*focus*sin(belta));
%         for f=(floor(Nfft*1/16)):(floor(Nfft*3/16))     % frequency: 500 --- 3000
%             B=[exp(i*2*pi*f/Nfft*d1*fs/V),exp(i*2*pi*f/Nfft*d2*fs/V),exp(i*2*pi*f/Nfft*d3*fs/V)];
%             Y(B1,f)=B(1)*X1(f)+B(2)*X2(f)+B(3)*X3(f);
%             Z(B1,R1)=Z(B1,R1)+abs(Y(B1,f))^2;
%         end
%         if Z(B1,R1)>z(B1)
%             z(B1)=Z(B1,R1);
%             r(B1)=R1;
%         end               
%     end    
% end
% 
% % Source localization
% max=0;
% for angle=1:19    
%     if z(angle)>max
%        max=z(angle);
%        sita1=angle;
%        r1=r(angle);
%     end
% end
% submax=0;
% for angle=1:19    
%     if (z(angle)>submax) & (z(angle)<max)
%        submax=z(angle);
%        sita2=angle;
%        r2=r(angle);
%     end
% end
%         
% sita1= (sita1-1)*pi/18-pi/2;
% r1= (r1-1)*0.1+Nbeg;
% sita2= (sita2-1)*pi/18-pi/2;
% r2= (r2-1)*0.1+Nbeg;
% Location1=[sita1*180/pi,r1]
% Location2=[sita2*180/pi,r2]
% 
% % 3D lacalization beam representation based on angle and rangle
% figure;
% X=-90:10:90;
% Y=Nbeg:0.1:Nend;
% % Z=20*log10(Z);
% surf(Y(1:No),X(1:19),Z(1:19,1:No))
% xlabel('range[m]')
% ylabel('Direction[Deg]')
% zlabel('Power')
% 
% % 2D localization beam representation based on angle
% figure;
% x=-90:10:90;
% % z=20*log10(z);
% plot(x,z,'-b');
% xlabel('Direction[Deg]')
% ylabel('Power')

sita1=-pi/3;  % the angle of the first source.  range: from -pi/2 to pi/2
sita2=pi/3;   % the angle of the second source. range: from -pi/2 to pi/2
r1=1.0;         % the radius of the first sound source
r2=1.0;         % the radius of the second sound source
V=344.21;
% sounds Separation based on sound localization
for f=1:Nfft
      d11=sqrt(D*D/4+r1*r1+D*r1*sin(sita1));
      d12=sqrt(3*D*D/4+r1*r1-sqrt(3)*D*r1*cos(sita1));
      d13=sqrt(D*D/4+r1*r1-D*r1*sin(sita1));
      d21=sqrt(D*D/4+r2*r2+D*r2*sin(sita2));
      d22=sqrt(3*D*D/4+r2*r2-sqrt(3)*D*r2*cos(sita2));
      d23=sqrt(D*D/4+r2*r2-D*r2*sin(sita2));
    A=[exp(-i*2*pi*f/Nfft*d11*fs/V)/d11,exp(-i*2*pi*f/Nfft*d21*fs/V)/d21;
       exp(-i*2*pi*f/Nfft*d12*fs/V)/d12,exp(-i*2*pi*f/Nfft*d22*fs/V)/d22;
       exp(-i*2*pi*f/Nfft*d13*fs/V)/d13,exp(-i*2*pi*f/Nfft*d23*fs/V)/d23];
    X=zeros(3,1);
    X(1)=X1(f);
    X(2)=X2(f);
    X(3)=X3(f);
    Sr(1:2,f)=inv(conj(A')*A)*conj(A')*X;
end

St1=real(ifft(Sr(1,:)));
St2=real(ifft(Sr(2,:)));

% plot the separation signals in time domain
figure;
subplot(211);
plot(s1,'-r');
xlabel('Time domain points')
ylabel('Amplitude')
title('Sound source1')
subplot(212);
plot(s2,'-b');
xlabel('Time domain points')
ylabel('Amplitude')
title('Sound source2')

figure;
subplot(211);
plot(St1,'-r');
xlabel('Time domain points')
ylabel('Amplitude')
title('Separated source1')
subplot(212);
plot(St2,'-b');
xlabel('Time domain points')
ylabel('Amplitude')
title('Separated source2')

St1=(St1-min(St1))/(max(St1)-min(St1));
St2=(St2-min(St2))/(max(St2)-min(St2));

fs=16000;
N=8;
auwrite(mixed1,fs,N,'mixture1');
auwrite(mixed2,fs,N,'mixture2');
auwrite(mixed3,fs,N,'mixture3');

auwrite(St1,fs,N,'St1');
auwrite(St2,fs,N,'St2');

% % plot the separation signals in frequency domain
% figure;
% subplot(221);
% plot(20*log10(abs(S1)),'-r');
% xlabel('Frequency domain points')
% ylabel('Amplitude')
% title('Sound source1')
% axis([1 Nfft -100 0])
% grid on;
% subplot(222);
% plot(20*log10(abs(S2)),'-b');
% xlabel('Frequency domain points')
% ylabel('Amplitude')
% title('Sound source2')
% axis([1 Nfft -100 0])
% grid on;
% subplot(223);
% plot(20*log10(abs(Sr(1,:))),'-r');
% xlabel('Frequency domain points')
% ylabel('Amplitude')
% title('Estimate of source1')
% axis([1 Nfft -100 0])
% grid on;
% subplot(224);
% plot(20*log10(abs(Sr(2,:))),'-b');
% xlabel('Frequency domain points')
% ylabel('Amplitude')
% title('Estimate of source2')
% axis([1 Nfft -100 0])
% grid on;

% [1] Sound source localization and signal separatiion for office robot
% "Jijo-2" --- Japan.
% [2] A LEARNING ALGORITHM WITH ADAPTIVE EXPONENTIAL STEPSIZE FOR BLIND
% SOURCE SEPARATION OF CONVOLUTIVE MIXTURES WITH REVERBERATIONS