www.gusucode.com > wavelet工具箱matlab源码程序 > wavelet/wavelet/@laurpoly/euclidediv.m
function DEC = euclidediv(A,B) %EUCLIDEDIV Euclidean Algorithm for Laurent polynomials. % DEC = EUCLIDEDIV(A,B) returns a cell array of Laurent polynomials % such that each row of DEC contains an euclidian division of A % by B: A = B*Q + R. Q is the quotient and R the remainder. % For each j from 1 to size(DEC,1): A = B*DEC{j,1} + DEC{j,2} % % The cell array DEC contains at most four rows. % % Example: % % Create two Laurent polynomials % A = laurpoly((1:4),0); % B = laurpoly([1 2],0); % % % Euclidian division of A by B % DEC = euclidediv(A,B); % for k=1:4 % disp(DEC{1,1},'Q'); disp(DEC{1,2},'R') % end % % %---------------------------------------------------------------- % % A(z) = 1 + 2*z^(-1) + 3*z^(-2) + 4*z^(-3) and % % B(z) = 1 + 2*z^(-1) % % There are four decomposition A = B*Q + R: % % Q(z) = 1 + 3*z^(-2) and R(z) = - 2*z^(-3) % % Q(z) = 1 + 2*z^(-2) and R(z) = z^(-2) % % Q(z) = 1 + 0.5*z^(-1) + 2*z^(-2) and R(z) = - 0.5*z^(-1) % % Q(z) = 0.75 + 0.5*z^(-1) + 2*z^(-2) and R(z) = 0.25 % %----------------------------------------------------------------- % M. Misiti, Y. Misiti, G. Oppenheim, J.M. Poggi 20-Mar-2001. % Last Revision 20-Dec-2010. % Copyright 1995-2010 The MathWorks, Inc. if isnumeric(A) && length(A)==1 , A = laurpoly(A,0); elseif isnumeric(B) && length(B)==1 , B = laurpoly(B,0); end maxDEG_A = A.maxDEG; maxDEG_B = B.maxDEG; cA = A.coefs; lA = length(cA); cB = B.coefs; lB = length(cB); minDEG_A = maxDEG_A-lA+1; minDEG_B = maxDEG_B-lB+1; maxDEG_LEFT = maxDEG_A-maxDEG_B; minDEG_RIGHT = minDEG_A-minDEG_B; dL = lA-lB; if dL>0 rLEFT = cA; qLEFT = zeros(1,dL+1); for j=1:dL+1 idxBEG = j; idxEND = idxBEG+lB-1; q = rLEFT(idxBEG)/cB(1); if j==(dL+1) qBIS = rLEFT(idxEND)/cB(end); rLEFT_2 = rLEFT; rLEFT_2(idxBEG:idxEND) = rLEFT_2(idxBEG:idxEND)-qBIS*cB; qLEFT_2 = [qLEFT(1:dL),qBIS]; end % qLEFT = [qLEFT q]; qLEFT(j) = q; rLEFT(idxBEG:idxEND) = rLEFT(idxBEG:idxEND)-q*cB; end rRIGHT = cA; qRIGHT = []; for j=1:dL+1 idxEND = lA-j+1; idxBEG = idxEND-lB+1; q = rRIGHT(idxEND)/cB(end); if j==(dL+1) qBIS = rRIGHT(idxBEG)/cB(1); rRIGHT_2 = rRIGHT; rRIGHT_2(idxBEG:idxEND) = rRIGHT_2(idxBEG:idxEND)-qBIS*cB; qRIGHT_2 = [qBIS,qRIGHT]; end qRIGHT = [q qRIGHT]; rRIGHT(idxBEG:idxEND) = rRIGHT(idxBEG:idxEND)-q*cB; end maxDEG_RIGHT = minDEG_RIGHT+length(qRIGHT)-1; DEC = {... laurpoly(qLEFT,maxDEG_LEFT) , laurpoly(rLEFT,maxDEG_A); ... laurpoly(qLEFT_2,maxDEG_LEFT) , laurpoly(rLEFT_2,maxDEG_A); ... laurpoly(qRIGHT_2,maxDEG_LEFT) , laurpoly(rRIGHT_2,maxDEG_A); ... laurpoly(qRIGHT,maxDEG_RIGHT) , laurpoly(rRIGHT,maxDEG_A) ... }; elseif dL==0 qLEFT = cA(1)/cB(1); qRIGHT = cA(end)/cB(end); rLEFT = cA-qLEFT*cB; rRIGHT = cA-qRIGHT*cB; maxDEG_RIGHT = minDEG_RIGHT; DEC = {... laurpoly(qLEFT,maxDEG_LEFT) , laurpoly(rLEFT,maxDEG_A); ... laurpoly(qRIGHT,maxDEG_RIGHT) , laurpoly(rRIGHT,maxDEG_A) ... }; else Q = laurpoly(0,0); %#ok<*NASGU> R = A; DEC = {laurpoly(0,0),A}; end nbDEC = size(DEC,1); idx = true(1,nbDEC); for j=1:nbDEC for k=j+1:nbDEC if idx(k)==1 idx(k) = (DEC{j,1} ~= DEC{k,1}) | (DEC{j,2} ~= DEC{k,2}); end end end DEC = DEC(idx,:);