mvar.m

matlab数字信号处理工具箱

        %%%%% multi-channel Levinsion algorithm 
        %%%%% using Nutall-Strand Method [2]
        %%%%% Covariance matrix is normalized by N=length(X) 
        
        %C{1} = C{1}/N;
        F = Y;
        B = Y;
        PEF = C(:,1:M);
        PEB = C(:,1:M);
        for K=1:Pmax,
                [D,n]  = covm(F(K+1:N,:),B(1:N-K,:),'M');
                %D=D/N;

                ARF(:,K*M+(1-M:0)) = D / PEB;	
                ARB(:,K*M+(1-M:0)) = D'/ PEF;	
                
                tmp        = F(K+1:N,:) - B(1:N-K,:)*ARF(:,K*M+(1-M:0)).';
                B(1:N-K,:) = B(1:N-K,:) - F(K+1:N,:)*ARB(:,K*M+(1-M:0)).';
                F(K+1:N,:) = tmp;
                
                for L = 1:K-1,
                        tmp      = ARF(:,L*M+(1-M:0))   - ARF(:,K*M+(1-M:0))*ARB(:,(K-L)*M+(1-M:0));
                        ARB(:,(K-L)*M+(1-M:0)) = ARB(:,(K-L)*M+(1-M:0)) - ARB(:,K*M+(1-M:0))*ARF(:,L*M+(1-M:0));
                        ARF(:,L*M+(1-M:0))   = tmp;
                end;
                
                RCF(:,K*M+(1-M:0)) = ARF(:,K*M+(1-M:0));
                RCB(:,K*M+(1-M:0)) = ARB(:,K*M+(1-M:0));
                
                [PEB,n] = covm(B(1:N-K,:),B(1:N-K,:),'M');
                %PEB = D/N;

                [PEF,n] = covm(F(K+1:N,:),F(K+1:N,:),'M');
                %PEF = D/N;

                PE(:,K*M+(1:M)) = PEF;        
        end;
        
elseif Mode==3, %%%%% multi-channel Levinsion algorithm [2] using Vieira-Morf Method
        fprintf('Warning MDURLEV: It''s not recommended to use this mode\n')        
        C(:,1:M) = C(:,1:M)/N;
        F = Y;
        B = Y;
        PEF = C(:,1:M);
        PEB = C(:,1:M);
        for K=1:Pmax,
                [D, n]  = covm(F(K+1:N,:),B(1:N-K,:),'M');
                D = D./n;

		ARF(:,K*M+(1-M:0)) = (PEF.^-.5)*D*(PEB.^-.5)';	
                ARB(:,K*M+(1-M:0)) = ARF(:,K*M+(1-M:0)); 
                
                tmp        = F(K+1:N,:) - B(1:N-K,:)*ARF(:,K*M+(1-M:0)).';
                B(1:N-K,:) = B(1:N-K,:) - F(K+1:N,:)*ARB(:,K*M+(1-M:0)).';
                F(K+1:N,:) = tmp;
                
                for L = 1:K-1,
                        tmp      = ARF(:,L*M+(1-M:0))   - ARF(:,K*M+(1-M:0))*ARB(:,(K-L)*M+(1-M:0));
                        ARB(:,(K-L)*M+(1-M:0)) = ARB(:,(K-L)*M+(1-M:0)) - ARB(:,K*M+(1-M:0))*ARF(:,L*M+(1-M:0));
                        ARF(:,L*M+(1-M:0))   = tmp;
                end;
                
                %RCF{K} = ARF{K};
                RCF = ARF(:,K*M+(1-M:0));
                
                [PEF,n] = covm(F(K+1:N,:),F(K+1:N,:),'M');
                PEF = PEF./n;

		[PEB,n] = covm(B(1:N-K,:),B(1:N-K,:),'M');
                PEB = PEB./n;

                %PE{K+1} = PEF;        
                PE(:,K*M+(1:M)) = PEF;        
        end;
        
elseif Mode==7 %%%%% multi-channel Levinsion algorithm [2] using Vieira-Morf Method
        fprintf('Warning MDURLEV: It''s not recommended to use this mode\n')        
        C(:,1:M) = C(:,1:M)/N;
        F = Y;
        B = Y;
        PEF = C(:,1:M);
        PEB = C(:,1:M);
        for K=1:Pmax,
                [D,n]  = covm(F(K+1:N,:),B(1:N-K,:),'M');
                D = D./n;

		ARF(:,K*M+(1-M:0)) = (PEF.^-.5)*D*(PEB.^-.5);	
                ARB(:,K*M+(1-M:0)) = (PEF.^-.5)*D'*(PEB.^-.5);	
                
                tmp        = F(K+1:N,:) - B(1:N-K,:)*ARF(:,K*M+(1-M:0)).';
                B(1:N-K,:) = B(1:N-K,:) - F(K+1:N,:)*ARB(:,K*M+(1-M:0)).';
                F(K+1:N,:) = tmp;
                
                for L = 1:K-1,
                        tmp      = ARF(:,L*M+(1-M:0))   - ARF(:,K*M+(1-M:0))*ARB(:,(K-L)*M+(1-M:0));
                        ARB(:,(K-L)*M+(1-M:0)) = ARB(:,(K-L)*M+(1-M:0)) - ARB(:,K*M+(1-M:0))*ARF(:,L*M+(1-M:0));
                        ARF(:,L*M+(1-M:0))   = tmp;
                end;
                
                %RCF{K} = ARF{K};
                RCF = ARF(:,K*M+(1-M:0));
                
                [PEF,n] = covm(F(K+1:N,:),F(K+1:N,:),'M');
                PEF = PEF./n;

                [PEB,n] = covm(B(1:N-K,:),B(1:N-K,:),'M');
                PEB = PEB./n;

                %PE{K+1} = PEF;        
                PE(:,K*M+(1:M)) = PEF;        
        end;
        
elseif Mode==4,  %%%%% nach Kay, not fixed yet. 
        fprintf('Warning MDURLEV: It''s not recommended to use this mode\n')        
        
        C(:,1:M) = C(:,1:M)/N;
        K = 1;
        [C(:,M+(1:M)),n] = covm(Y(2:N,:),Y(1:N-1,:));
        C(:,M+(1:M)) = C(:,M+(1:M))/N;  % biased estimate

        D = C(:,M+(1:M));
        ARF(:,1:M) = C(:,1:M)\D;
        ARB(:,1:M) = C(:,1:M)\D';
        RCF(:,1:M) = ARF(:,1:M);
        RCB(:,1:M) = ARB(:,1:M);
        PEF = C(:,1:M)*[eye(M) - ARB(:,1:M)*ARF(:,1:M)];
        PEB = C(:,1:M)*[eye(M) - ARF(:,1:M)*ARB(:,1:M)];
        
        for K=2:Pmax,
                [C(:,K*M+(1:M)),n] = covm(Y(K+1:N,:),Y(1:N-K,:),'M');
                C(:,K*M+(1:M)) = C(:,K*M+(1:M)) / N; % biased estimate

                D = C(:,K*M+(1:M));
                for L = 1:K-1,
                        D = D - ARF(:,L*M+(1-M:0))*C(:,(K-L)*M+(1:M));
                end;
                ARF(:,K*M+(1-M:0)) = PEB \ D;	
                ARB(:,K*M+(1-M:0)) = PEF \ D';	
                for L = 1:K-1,
                        ARFtmp(:,L*M+(1-M:0))  = ARF(:,L*M+(1-M:0))  - ARB(:,(K-L)*M+(1-M:0)) *ARF(:,K*M+(1-M:0)) ;
                        ARB(:,L*M+(1-M:0))  = ARB(:,L*M+(1-M:0))  - ARF(:,(K-L)*M+(1-M:0)) *ARB(:,K*M+(1-M:0)) ;
                end;
                ARF(:,1:(K-1)*M) = ARFtmp;
                RCF(:,K*M+(1-M:0)) = ARF(:,K*M+(1-M:0)) ;
                RCB(:,K*M+(1-M:0)) = ARB(:,K*M+(1-M:0)) ;
                
                PEF = PEF*[eye(M) - ARB(:,K*M+(1-M:0)) *ARF(:,K*M+(1-M:0)) ];
                PEB = PEB*[eye(M) - ARF(:,K*M+(1-M:0)) *ARB(:,K*M+(1-M:0)) ];
                PE(:,K*M+(1:M))  = PEF;        
        end;
        
elseif Mode==8,  %%%%% MAXIMUM LIKELIHOOD ESTIMATE [Larbi and Lardies, 2000]
        % step 1;
        [A,RCF,PEF] = mvar(Y,Pmax,2);         % initial values
        
        % step 2
        for i = 1:Pmax,
        %        tmp  = mvfilter(eye(M),[eye(M),-A(:,1:i*M)], Y');
        %        u{i} = tmp;
        end;
	
        ut  = mvfilter(eye(M),[eye(M),-A], Y');
        
	S = covm(ut','M');		% (33)
	SI = inv(S);
	
        % step 3                
        du = 0; 
        R = 1; % ???
	y1 = zeros(M,Pmax);
	I = 0;
        dlL = 0;
        for t = 1:size(Y,1),
		y1(:,1:min(t-1,Pmax)) = Y(t-1:-1:max(t-Pmax,1),:)';
		a = zeros(1,M);
		for k1 = 1:min(Pmax,t), 
		if t>k1,
			a = a+Y(t-k1,:)*A(:,(1-M:0)+k1*M)';
		end;	
                end;
		dut = kron(a,eye(M))*[eye(M*M),zeros(M*M,M*M*Pmax)]*R
		dut1 = kron(y1(:).',eye(M))*[zeros(M*M*Pmax,M*M),eye(M*M*Pmax)]*R
		dut = dut - dut1;
        
		size(dut),size(SI),
		I = I + dut'*SI*dut;		% (36)
		
		%[t,size(dlL),size(ut),size(dut)],
		dlL = dlL + ut(:,t)'*SI*dut;	% (37)
	end;        
        
        % step 4                
        
        % step 5       
	[size(I),size(dlL)]         
	G = inv(I)*dlL'
	
        
	ARF = A;
        
elseif Mode==9,  %%%%% MAXIMUM LIKELIHOOD ESTIMATE [PHAM and TONG, 1994]
        
        if 1, %while 1, 
                [A,RCF,PEF] = mvar(Y,Pmax,6);         % initial values
                U = PEF(:,K*M+(1:M));        
                G = [eye(M), -A]\[PE, zeros(M,M*Pmax)]
                
                [B,RCF,PEF] = mvar(Y(end:-1:1,:),Pmax,6);         % initial values
                V = PEF(:,K*M+(1:M));     
                G = [eye(M), -B]\[zeros(M,M*Pmax), PE]
                
                for K=0:Pmax,
                        [R,n] = covm(Y(K+1:N,:),Y(1:N-K,:),'M');
                end;        
                ARF = A;
        end;
end;

if any(ARF(:)==inf),
% Test for matrix division bug. 
% This bug was observed in LNX86-ML5.3, 6.1 and 6.5, but not in SGI-ML6.5, LNX86-ML6.5, Octave 2.1.35-40; Other platforms unknown.
p = 3;
FLAG_MATRIX_DIVISION_ERROR = ~all(all(isnan(repmat(0,p)/repmat(0,p)))) | ~all(all(isnan(repmat(nan,p)/repmat(nan,p))));

if FLAG_MATRIX_DIVISION_ERROR, 
	%fprintf(2,'### Warning MVAR: Bug in Matrix-Division 0/0 and NaN/NaN yields INF instead of NAN.  Workaround is applied.\n');
	warning('MVAR: bug in Matrix-Division 0/0 and NaN/NaN yields INF instead of NAN.  Workaround is applied.');

	%%%%% Workaround 
	ARF(ARF==inf)=NaN;
	RCF(RCF==inf)=NaN;
end;
end;	

%MAR   = zeros(M,M*Pmax);
DC     = zeros(M);
for K  = 1:Pmax,
%       VAR{K+1} = -ARF(:,K*M+(1-M:0))';
        DC  = DC + ARF(:,K*M+(1-M:0))'.^2; %DC meausure [3]
end;