ltfat-wfilt_db-ccdi.patch

diff --git a/wavelets/wfilt_db.m b/wavelets/wfilt_db.m
index 35c3cada..7016739a 100644
--- a/wavelets/wfilt_db.m
+++ b/wavelets/wfilt_db.m
@@ -1,4 +1,4 @@
-function [h, g, a, info] = wfilt_db(N)
+function [h, g, a, info] = wfilt_db(N, coe = 0, compat = 1, dedeg = 0, il = 1)
 %WFILT_DB    Daubechies FIR filterbank
 %   Usage:  [h,g] = wfilt_db(N);
 %
@@ -15,7 +15,7 @@ function [h, g, a, info] = wfilt_db(N)
 %   spectral factorization of the Lagrange interpolator filter.  *N* also
 %   denotes the number of zeros at $z=-1$ of the lowpass filters of length
 %   $2N$.  The prototype lowpass filter has the following form (all roots of
-%   $R(z)$ are outside of the unit circle):
+%   $R(z)$ are inside of the unit circle):
 %                 
 %   .. H_l(z)=(1+z^-1)^N*R(z),
 %
@@ -53,36 +53,153 @@ end
 if(N<1)
     error('Minimum N is 1.');
 end
-if(N>20)
-    warning('Instability may occur when N is too large.');
-end
+%if(N>20)
+%    warning('Instability may occur when N is too large.');
+%end
 
 
 h = cell(2,1);
 flen = 2*N;
 
 % Calculating Lagrange interpolator coefficients
-sup = [-N+1:N];
-a = zeros(1,N);
-for n = 1:N
-    non  = sup(sup ~= n);
-    a(n) = prod(0.5-non)/prod(n-non);
-end
-P = zeros(1,2*N-1);
-P(1:2:end) = a;
-P = [P(end:-1:1),1,P];
+if (coe == 0);
+	sup = [-N+1:N];
+	a   = zeros(1,N);
+	for n = 1:N
+		non  = sup(sup ~= n);
+		a(n) = prod(0.5-non)/prod(n-non);
+	end;
+	P = zeros(1,2*N-1);
+	P(1:2:end) = a;
+	P = [P(end:-1:1),1,P];
+else
+	sup = [-N+1:N];
+	s   = zeros(1,N);
+	b   = zeros(1,N);
+	p   = zeros(1,N);
+	for n = 1:N
+		non  = sup(sup ~= n);
+		n2 = 1 - non * 2; ln2 = length(n2); sn2 = mod(length(n2(n2 < 0)), 2); an2 = abs(n2);
+		fan2 = [];
+		for i = 1:ln2;
+			fan2 = [fan2, factor(an2(i))];
+		end;
+		fan2 = sort(fan2, 'ascend'); lfan2 = length(fan2);
+		d  = n - non;     ld  = length(d ); sd  = mod(length(d (d  < 0)), 2); ad  = abs(d );
+		fad = [];
+		for i = 1:ld;
+			fad = [fad, factor(ad(i))];
+		end;
+		fad = sort(fad, 'ascend'); lfad = length(fad);
+		rfan2 = [1]; i = 1;
+		rfad  = [1]; j = 1;
+		while (1);
+			if (fan2(i) == 1);
+				i++;
+			elseif (fad(j) == 1);
+				j++;
+			elseif (fan2(i) == fad(j));
+				i++; j++;
+			elseif (fan2(i) < fad(j));
+				rfan2 = [rfan2, fan2(i)]; i++;
+			elseif (fan2(i) > fad(j));
+				rfad  = [rfad, fad(j)]; j++;
+			end;
+			if (i > lfan2 || j > lfad);
+				break;
+			end;
+		end;
+		if (i <= lfan2);
+			rfan2 = [rfan2, fan2(i:end)];
+		end;
+		if (j <= lfad);
+			rfad = [rfad, fad(j:end)];
+		end;
+		lrfan2 = length(rfan2          );
+		lrfad  = length(rfad           );
+		lrfad1 = length(rfad(rfad == 1));
+		lrfad2 = length(rfad(rfad == 2));
+		assert(lrfad == lrfad1 + lrfad2);
+		s(n) = mod(sn2 + sd, 2);
+		t    = floor(log2(rfan2));
+		p(n) = sum(log2(rfan2 .* pow2(-t)));
+		b(n) = ln2 - sum(t) + lrfad2;
+	end;
+	P = zeros(1,2*N-1);
+	P(1:2:end) = ((-1).^s) .* pow2(p - b);
+	P = [P(end:-1:1),1,P];
+end;
 
 R = roots(P);
-% Roots outside of the unit circle and in the right halfplane
-R = R(abs(R)<1 & real(R)>0);
-
-% roots of the 2*conv(lo_a,lo_r) filter
-hroots = [R(:);-ones(N,1)];
-
-
-% building synthetizing low-pass filter from roots
-h{1}= real(poly(sort(hroots,'descend')));
-% normalize
+% Roots inside of the unit circle and in the right halfplane
+Rur = R(abs(R)<1 & real(R)>0);
+
+% Building synthesizing low-pass filter from roots
+if (compat == 1);
+	c = poly(sort([Rur(:);-ones(N,1)],'descend'));
+else
+	% Build c is not having problematic imaginary parts, especially for large N.
+	% c is polynomial w/ complex coefficients if w/o any consideration of Rur.
+	% Consider Rur is N-1 length complex vector, but for all 1 <= i <= N-1,
+	%  Rur(i) is real, or
+	%  Rur(j) exists where 1 <= j <= N-1 and j != i and Rur(j) == conj(Rur(i))
+	% is satisfied.
+	% This means c is actually polynominal w/ real coefficients.
+	c = [1];
+	for i = 1:N;
+		c = conv(c, [1, 1]);
+	end;
+	Rzp = Rur(imag(Rur) >= 0); Rzp = sort(Rzp, 'ascend'); lRzp = length(Rzp);
+	% FIXME: dedeg by unused R?
+	%if (dedeg == 1);
+	%	Ruu = R(!(abs(R)<1 & real(R)>0)); lRuu = length(Ruu);
+	%	% FIXME: Add Ruu refining?
+	%	if (lRuu > 0);
+	%		P = deconv(P, poly(sort(Ruu, 'descend')));
+	%	end;
+	%end;
+	% Rzp refininig by Newton-Raphson method bounded by macihne epsilon or iteration limit.
+	for i = 1:lRzp;
+		if (i == 1 || dedeg == 1);
+			P  = P / P(1);   P(1) = 1; lP  = length(P );
+			Pd = polyder(P);           lPd = length(Pd);
+		end;
+		if (lP >= 2 && lPd >= 1);
+			d  = eps * max(abs(P (2:lP ) .* (Rzp(i) .^ (lP -2:-1:0)))); d  *= lP ;
+			dd = eps * max(abs(Pd(1:lPd) .* (Rzp(i) .^ (lPd-1:-1:0)))); dd *= lPd;
+			p  = polyval(P,  Rzp(i));
+			pd = polyval(Pd, Rzp(i)); %assert(abs(pd) >  dd); % FIXME: pd actually should be far from 0.
+			e  = d / (abs(pd) - dd);  %assert(    e   >=  0);
+			Rzpi = Rzp(i);
+			j = 0;
+			 while (abs(p) > e && j < il) % FIXME: e is too pessimistic?
+			%while (              j < il)
+				Rzpi = [Rzpi, Rzpi(end) - p / pd];
+				d  = eps * max(abs(P (2:lP ) .* (Rzpi(end) .^ (lP -2:-1:0)))); d  *= lP ;
+				dd = eps * max(abs(Pd(1:lPd) .* (Rzpi(end) .^ (lPd-1:-1:0)))); dd *= lPd;
+				p  = polyval(P,  Rzpi(end));
+				pd = polyval(Pd, Rzpi(end)); %assert(abs(pd) >  dd); % FIXME: pd actually should be far from 0.
+				e  = d / (abs(pd) - dd);     %assert(    e   >=  0);
+				j++;
+			end;
+			[~, kmin] = min(abs(polyval(P, Rzpi))); Rzp(i) = Rzpi(kmin);
+		end;
+		if (imag(Rzp(i)) == 0);
+			Q = [1, -real(Rzp(i))];
+		else
+			Q = [1, -2*real(Rzp(i)), abs(Rzp(i))^2];
+		end;
+		c = conv(c, Q);
+		if (dedeg == 1);
+			P = deconv(P, Q);
+		end;
+	end;
+end;
+% Following assertion may fail if compat = 1 and N >= 5.
+%maic = max(abs(imag(c))); printf('maic = %.15e.\n', maic); assert(maic <= eps);
+h{1} = real(c);
+
+% Normalize
 h{1}= h{1}/norm(h{1});
 % QMF modulation low-pass -> highpass
 h{2}= (-1).^(0:flen-1).*h{1}(end:-1:1);