利用经典联合对角化方法实现盲卷积语音分离，可处理超定及欠定的盲卷积语音分离问题matlab.zip

function [y,W] = alcobs(x, Ns, nFFT, nP, nSample, ovl) % Alternating Least Squares Convolutive Blind Source Separation % Version 1.0 % Author: Kamran Rahbar % Institue: McMaster University (Statistical Signal processing Lab) % This function does blind source separation of % convolutive mixtures useing a frequency domain approach. % For more info. see: % Kamran Rahbar, James P. Reilly, % "A Frequency Domain Method for Blind Source Separation % of Convolutive Audio Mixtures", % IEEE Trans. on Speech and Audio Processing (to be published). % % Please report any bug or problem to reillyj@mcmaster.ca % % usage: % [y,w]=alcobs(x, Ns, nFFT, nP, nSample, ovl); % % x Observed signals from convolutive mixtures % Ns Number of sources % nFFT Number of FFT points (usually greather than the length of channel) % nP Number of cross spectral density matrices % nSample Number of samples in each epoch % ovl Overlap percetange (0<ovl<1) represents the amount of overlap between adjacent FFT windows % % Note that Lx > (nP*nFFT)/ovl where Lx is the length of observed data % % % y Separated output signal with y(i,:) corresponds to ith signal % W Demixing matrix in time domain; W(i,j,t) representsb the ij element of demixing matrix % at time t % Example of Usage for data collected in a real room environmnet: % [y,W]=alcobs(x, 2, 4096, 100, 10000, 0.5); % % Copyright Kamran Rahbar & James P. Reilly % %THIS SOFTWARE IS PROVIDED BY THE AUTHOR "AS IS" AND ANY EXPRESS OR %IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED %WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE %HEREBY DISCLAIMED. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR %ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL %DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS %OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) %HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, %STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN %ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE %POSSIBILITY OF SUCH DAMAGE. % %----------- 2002 McMaster University - Kamran Rahbar ----------- % [mx, nx]=size(x); if mx > nx x=x'; [mx nx]=size(x); end if (mx < Ns) error('Ns should be less than or equal to num. of observed signals')' return; end if (ovl >1 | ovl<=0 ) error('0<ovl<1'); return; end if (fix(log2(nFFT))-log2(nFFT)~=0) error('nFFT should be in power of 2') return end if (nFFT> nSample) error('nSample should be >= nFFT') end shift=ceil((1-ovl)*nFFT); % if (nx < (nSample*nP+ceil(ovl*nSample))) % % error(' Lx > nSample*nP+ovl*nSample') % return; % end % Calculate the Cross Spectral Density Matrices [Px, x]=CalPx(x, shift, nFFT, nP, nSample); % Calculate H(w) at each frequency bin using ALS [H,D,Binv,Nit]=AlsBss(Px,Ns,nFFT); % remove permutations [HH,DD,per]=PermRemove(H,D); % Calculate the final output W=real(ifft(HH,nFFT,3)); y=mat_conv(W, x); %$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ function [Px, x]=CalPx(x, shift, nFFT, M, nSample) [mx, nx]=size(x); Len=nSample*M; if (Len > nx) error('nSample*M should be leq to number of observed data') return; end x=x(:,1:Len); %=========================remove the mean and normalize============================= for i=1:mx x(i,:)=x(i,:)-mean(x(i,:)); end x=diag(1./std(x'))*x; %=========Main program starts from here================================ Px=zeros(mx,mx,nFFT/2+1,M); kk=0; for z=0:nSample:Len-nSample kk=kk+1; ShowProgress(kk/M, 'Calculating CPSD matrices.....'); for v=0:shift:nSample-nFFT-1 i=z+v; for mmm=1:mx X(mmm,:)=fft(x(mmm,i+1:i+nFFT).*hanning(nFFT)',nFFT); end for k=1:1:nFFT/2+1 Px(:,:,k,kk)=(1/nFFT)*X(:,k)*X(:,k)'+Px(:,:,k,kk); end end end %$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ function [B,D,Binv,Nit]=AlsBss(Px,nSource,nFFT); [N1,N2,Kr,M]=size(Px); kk=0; eps=1e-8; delta=0.005; Kp=2*(Kr-1)/nFFT; K=nFFT/2+1; H_init=zeros(N1,nSource); H_init(1:nSource,1:nSource)=eye(nSource); for k=1:Kp:Kr kk=kk+1; ShowProgress(k/Kr, 'Calculating B(w) using ALS.....'); R=[]; P_sum=zeros(N1,N1); for i=1:M R(:,:,i)=Px(:,:,k,i)/(norm(Px(:,:,k,i),'fro')+1e-8); end R(:,:,1)=R(:,:,1)+delta*eye(N1); if kk==1 [Binv(:,:,kk),D(:,:,:,kk),Nit(kk),cost]=als(R,nSource,eps,H_init); else [Binv(:,:,kk),D(:,:,:,kk),Nit(kk),cost]=als(R,nSource,eps,Binv(:,:,kk-1)); end B(:,:,kk)=pinv(Binv(:,:,kk)); B(:,:,kk)=B(:,:,kk)/norm(B(:,:,kk),'fro'); end clc disp('Final Adjusments......'); for k=2:(K-1) Binv(:,:,2*K-k)=conj(Binv(:,:,k)); B(:,:,2*K-k)=conj(B(:,:,k)); end Delay=K-1; for k=1:2*(K-1) B(:,:,k)=B(:,:,k)*exp(2*pi*1j*Delay*(k-1)/(2*(K-1))); end %$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ function [H,D,Nit,cost]=als(P,nSource,thresh,H_init); [N,N,M]=size(P); MaxNumItr=60; Psi=zeros(N*N,nSource); E=zeros(N*N,M); d=zeros(nSource,M); tmp=zeros(N*N,M); for i=1:nSource G(:,:,i)=H_init(:,i)*H_init(:,i)'; G1=G(:,:,i); Psi(:,i)=G1(:); end for m=1:M P1=P(:,:,m); Pv(:,m)=P1(:); d(:,m)=pinv(Psi)*Pv(:,m); end H=H_init; cost(1)=0; for ii=1:MaxNumItr cost(ii+1)=0; for i=1:nSource for m=1:M tmp(:,m)=Pv(:,m)-Psi*d(:,m); E(:,m)=tmp(:,m)+Psi(:,i)*d(i,m); end Psi(:,i)=E*d(i,:)'/norm(d(i,:))^2; q_init=H(:,i); G1=mat(Psi(:,i)); z_init=G1*q_init; V_max=z_init/norm(z_init); H(:,i)=V_max; G1=V_max*V_max'; Psi(:,i)=G1(:); end for m=1:M d(:,m)=pinv(Psi)*Pv(:,m); D(:,:,m)=diag(d(:,m)); cost(ii+1)=norm(tmp(:,m))^2+cost(ii+1); end if abs(cost(ii+1)-cost(ii))/(cost(2)+1e-8)<thresh break; end end Nit=ii; %$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ function [H,D,per]=PermRemove(H,D) [N1,N2,K]=size(H); [N1,N2,M,KKK]=size(D); Pindx=perms(1:N1); I=eye(N1); for i=1:length(Pindx) Permute_Table(:,:,i)=I(:,Pindx(i,:)); end KK=K/2+1; Level=log2(K/2); Htmp=H; per=zeros(1,KK); ii=0; D_ave=zeros(N1,N2,M,2^(Level+1)-1); D_ave(:,:,:,1:KK-1)=D(:,:,:,1:KK-1); indx=0; %**************Level Zero********************** indx=KK-1; for k=1:2:KK-1 ShowProgress(k/(KK-1), 'removing permutations in level--0 ......'); indx=indx+1; [PI,flag]=FindPermute(D_ave(:,:,:,k),D_ave(:,:,:,k+1),Permute_Table); if flag>0 for i=1:M D_ave(:,:,i,k)=PI*D_ave(:,:,i,k)*PI'; end H(:,:,k)=PI*H(:,:,k); per(k)=abs(per(k)-1); end D_ave(:,:,:,indx)=D_ave(:,:,:,k)+D_ave(:,:,:,k+1); end %******************************************************************************** %**************Level >0 **************************** n=0; for level=1:Level-1 n=2^(Level-level+1)+n; for k=1:2:(KK-1)/(2^level) ShowProgress(k/((KK-1)/2^level), strcat('Removing permutations in level--', num2str(level), '......')); kk=n+k; indx=indx+1; if level==5 stop1=1; end [PI flag]=FindPermute(D_ave(:,:,:,kk),D_ave(:,:,:,kk+1),Permute_Table); if flag>0 [D_ave,H,per]=ChangePermute(kk,level,KK-1,D_ave,H,PI,per); end D_ave(:,:,: