diff -r b6c2035313e6 INDEX
--- a/INDEX	Sun May 19 17:15:31 2013 +0200
+++ b/INDEX	Sat Jun 15 20:23:36 2013 +0100
@@ -130,6 +130,7 @@
  rectwin
  welchwin
  parzenwin
+ ultrwin
 System identification
  arburg
  aryule
diff -r b6c2035313e6 inst/chebwin.m
--- a/inst/chebwin.m	Sun May 19 17:15:31 2013 +0200
+++ b/inst/chebwin.m	Sat Jun 15 20:23:36 2013 +0100
@@ -1,93 +1,79 @@
-## Copyright (C) 2002 André Carezia <acarezia@uol.com.br>
-##
-## This program is free software; you can redistribute it and/or modify it under
-## the terms of the GNU General Public License as published by the Free Software
-## Foundation; either version 3 of the License, or (at your option) any later
-## version.
-##
-## This program is distributed in the hope that it will be useful, but WITHOUT
-## ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-## FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
-## details.
-##
-## You should have received a copy of the GNU General Public License along with
-## this program; if not, see <http://www.gnu.org/licenses/>.
+%% Copyright (c) 2013 Rob Sykes <robs@users.sourceforge.net>
+%%
+%% This program is free software; you can redistribute it and/or modify it
+%% under the terms of the GNU General Public License as published by the Free
+%% Software Foundation: either version 3 of the License, or (at your option)
+%% any later version.  See http://gnu.org/licenses for details.
+%%
+%% This program is distributed in the hope that it will be useful, but WITHOUT
+%% ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+%% FITNESS FOR A PARTICULAR PURPOSE.
 
-## Usage:  chebwin (L, at)
-##
-## Returns the filter coefficients of the L-point Dolph-Chebyshev window
-## with at dB of attenuation in the stop-band of the corresponding
-## Fourier transform.  The default attenuation value is 100 dB.
-##
-## For the definition of the Chebyshev window, see
-##
-## * Peter Lynch, "The Dolph-Chebyshev Window: A Simple Optimal Filter",
-##   Monthly Weather Review, Vol. 125, pp. 655-660, April 1997.
-##   (http://www.maths.tcd.ie/~plynch/Publications/Dolph.pdf)
-##
-## * C. Dolph, "A current distribution for broadside arrays which
-##   optimizes the relationship between beam width and side-lobe level",
-##   Proc. IEEE, 34, pp. 335-348.
-##
-## The window is described in frequency domain by the expression:
-##
-##          Cheb(L-1, beta * cos(pi * k/L))
-##   W(k) = -------------------------------
-##                 Cheb(L-1, beta)
-##
-## with
-##
-##   beta = cosh(1/(L-1) * acosh(10^(at/20))
-##
-## and Cheb(m,x) denoting the m-th order Chebyshev polynomial calculated
-## at the point x.
-##
-## Note that the denominator in W(k) above is not computed, and after
-## the inverse Fourier transform the window is scaled by making its
-## maximum value unitary.
-##
-## See also: kaiser
+%% -*- texinfo -*-
+%% @deftypefn {Function File} {@var{w} =} chebwin (@var{l}, @var{att})
+%% Return as column-vector @var{w}, the coefficients of the @var{l}-point
+%% Dolph-Chebyshev window with side-lobe magnitude @var{att} dB below the
+%% main-lobe magnitude in its Fourier transform; normalise @var{w} such that the
+%% central coefficient(s) value is unitary.
+%%
+%% The un-normalised window can be described in the frequency domain as:
+%%
+%% @verbatim
+%%             W(k) = T   [ beta cos(pi k/L) ],  k = 0, 1, ..., L-1,
+%%                     L-1
+%% where
+%%                            ATT/20
+%%             T   (beta) = 10
+%%              L-1
+%% @end verbatim
+%%
+%% and T is the Chebyshev polynomial of the first kind, which can be defined as:
+%%
+%% @verbatim
+%%                      /  cos[ m arccos(x) ],   |x| <= 1 
+%%             T (x) = (   
+%%              m       \ cosh[ m arccosh(x) ],  |x| > 1 
+%% @end verbatim
+%%
+%% The Dolph-Chebyshev window is a special case of the Ultraspherical window.
+%% For further information about the window, see
+%%
+%% @itemize @bullet
+%% @item
+%% Lynch, P., 1997: The Dolph-Chebyshev Window: A Simple Optimal Filter.
+%% Monthly Weather Review, 125, pp. 655-660.
+%% @item
+%% Dolph, C. L., 1946: A current distribution for broadside arrays which
+%% optimizes the relationship between beam width and side-lobe level.
+%% Proc. IRE, 34, pp. 335-348.
+%% @end itemize
+%% @seealso{ultrwin, kaiser}
+%% @end deftypefn
 
-function w = chebwin (L, at)
-
-  if (nargin < 1 || nargin > 2)
-    print_usage;
-  elseif (nargin == 1)
-    at = 100;
+function w = chebwin (L, att = 100)
+  if !(nargin >= 1 && nargin <= 2 && nargout <= 1)
+    print_usage
+  elseif !(isscalar(L) && L == round(L) && L > 0)
+    error ('chebwin: L must be a positive integer');
+  elseif !(isscalar(att) && att == real(att))
+    error ('chebwin: ATT must be a real scalar');
   endif
 
-  if !(isscalar (L) && (L == round(L)) && (L > 0))
-    error ("chebwin: L has to be a positive integer");
-  elseif !(isscalar (at) && (at == real (at)))
-    error ("chebwin: at has to be a real scalar");
-  endif
-  
-  if (L == 1)
-    w = 1;
-  else
-				# beta calculation
-    gamma = 10^(-at/20);
-    beta = cosh(1/(L-1) * acosh(1/gamma));
-				# freq. scale
-    k = (0:L-1);
-    x = beta*cos(pi*k/L);
-				# Chebyshev window (freq. domain)
-    p = cheb(L-1, x);
-				# inverse Fourier transform
-    if (rem(L,2))
-      w = real(fft(p));
-      M = (L+1)/2;
-      w = w(1:M)/w(1);
-      w = [w(M:-1:2) w]';
-    else
-				# half-sample delay (even order)
-      p = p.*exp(j*pi/L * (0:L-1));
-      w = real(fft(p));
-      M = L/2+1;
-      w = w/w(2);
-      w = [w(M:-1:2) w(2:M)]';
-    endif
-  endif 
-  
-  w = w ./ max (w (:)); 
-endfunction
+  w = ultrwin(L, 0, att, 'a');
+end
+
+%!test
+%! L=201; m=L-1;
+%! for att=40:10:120
+%!   x=cosh(acosh(10^(att/20))/m)*cos([0:m]*pi/L);
+%!   b=real(ifft(cheb(L-1, x))); b=b(m/2+2:L)/b(1);
+%!   assert(chebwin(L, att)', [b 1 fliplr(b)], 1e-12);
+%! end
+
+%!demo
+%! w=chebwin(81, 60); [W,f]=freqz(w); clf
+%! subplot(2,1,1); plot(f/pi, 20*log10(W/abs(W(1)))); grid; axis([0 1 -90 0])
+%! subplot(2,1,2); plot(0:length(w)-1, w); grid
+%! %-----------------------------------------------------------
+%! % Figure shows a Dolph-Chebyshev window with 60dB attenuation:
+%! % frequency domain above, time domain below.
diff -r b6c2035313e6 inst/ultrwin.m
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/inst/ultrwin.m	Sat Jun 15 20:23:36 2013 +0100
@@ -0,0 +1,198 @@
+%% Copyright (c) 2013 Rob Sykes <robs@users.sourceforge.net>
+%%
+%% This program is free software; you can redistribute it and/or modify it
+%% under the terms of the GNU General Public License as published by the Free
+%% Software Foundation: either version 3 of the License, or (at your option)
+%% any later version.  See http://gnu.org/licenses for details.
+%%
+%% This program is distributed in the hope that it will be useful, but WITHOUT
+%% ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+%% FITNESS FOR A PARTICULAR PURPOSE.
+
+%% -*- texinfo -*-
+%% @deftypefn {Function File} {[@var{w}, @var{xmu}] =} ultrwin (@var{l}, @var{mu}, @var{beta})
+%% @deftypefnx {Function File} {[@var{w}, @var{xmu}] =} ultrwin (@var{l}, @var{mu}, @var{att}, 'att')
+%% @deftypefnx {Function File} {[@var{w}, @var{xmu}] =} ultrwin (@var{l}, @var{mu}, @var{latt}, 'latt')
+%% @deftypefnx {Function File} {@var{w} =} ultrwin (@var{l}, @var{mu}, @var{xmu}, 'xmu')
+%% Return as column-vector @var{w}, the coefficients of the @var{l}-point
+%% Ultraspherical window, where @var{mu} controls the window's Fourier
+%% transform's side-lobe to side-lobe ratio, and the third given parameter
+%% controls the transform's main-lobe width/side-lobe-ratio; normalise @var{w}
+%% such that the central coefficient(s) value is unitary.
+%%
+%% By default, the third parameter is @var{beta}, which sets the main lobe width
+%% to @var{beta} times that of a rectangular window.  Alternatively, giving
+%% @var{att} or @var{latt} sets the ripple ratio at the first or last side-lobe
+%% respectively, or giving @var{xmu} sets the (un-normalised) window's Fourier
+%% transform according to its canonical definition:
+%%
+%% @verbatim
+%%              (MU)
+%%      W(k) = C   [ XMU cos(pi k/L) ],  k = 0, 1, ..., L-1,
+%%              L-1
+%% @end verbatim
+%%
+%% where C is the Ultraspherical (a.k.a. Gegenbauer) polynomial, which can be
+%% defined using the recurrence relationship:
+%%
+%% @verbatim
+%%       (l)    1                  (l)                    (l)
+%%      C (x) = - [ 2x(m + l - 1) C   (x) - (m + 2l - 2) C   (x) ]
+%%       m      m                  m-1                    m-2
+%%
+%%                                 (l)        (l)
+%%      for m an integer > 1, and C (x) = 1, C (x) = 2lx.
+%%                                 0          1
+%% @end verbatim
+%%
+%% For given @var{beta}, @var{att}, or @var{latt}, the corresponding
+%% (determined) value of @var{xmu} is also returned.
+%%
+%% The Dolph-Chebyshev and Saramaki windows are special cases of the
+%% Ultraspherical window, with @var{mu} set to 0 and 1 respectively.  Note that
+%% when not giving @var{xmu}, stabilty issues may occur with @var{mu} <= 1.5.
+%% For further information about the window, see
+%%
+%% @itemize @bullet
+%% @item
+%% Kabal, P., 2009: Time Windows for Linear Prediction of Speech.
+%% Technical Report, Dept. Elec. & Comp. Eng., McGill University.
+%% @item
+%% Bergen, S., Antoniou, A., 2004: Design of Ultraspherical Window
+%% Functions with Prescribed Spectral Characteristics. Proc. JASP, 13/13,
+%% pp. 2053-2065.
+%% @item
+%% Streit, R., 1984: A two-parameter family of weights for nonrecursive
+%% digital filters and antennas. Trans. ASSP, 32, pp. 108-118.
+%% @end itemize
+%% @seealso{chebwin, kaiser}
+%% @end deftypefn
+
+function [w, xmu] = ultrwin (L, mu, par, key = 'BETA', norm = 0)
+  type = [strmatch(upper(key), types = ['XMU';'BETA';'ATT';'LATT'])' 0](1);
+  if !(type > 0 && nargin >= 3 && nargin <= 5 && nargout <= 2)
+    print_usage
+  elseif !(isscalar(L) && L == round(L) && L > 0)
+    error ('ultrwin: L must be a positive integer');
+  elseif !(isscalar(mu) && mu == real(mu))
+    error ('ultrwin: MU must be a real scalar');
+  elseif !(isscalar(par) && par == real(par))
+    error (horzcat('ultrwin: ', types(type,:), ' must be a real scalar'));
+  elseif !(isscalar(norm) && norm == round(norm) && norm >= 0) % Alt. norms; WIP
+    error ('ultrwin: NORM must be a non-negative integer');
+  endif
+
+  [w xmu] = ultrwin_(L, mu, par, type-1, norm);
+end
+
+%!test
+%! assert(ultrwin(100, 1, 1), ones(100, 1), 1e-14);
+
+%!test
+%! L = 201; xmu = 1.01; m = L-1;
+%! for mu = -1.35:.3:1.35
+%!   x = xmu*cos([0:m]*pi/L);
+%!   C(2,:) = 2*mu*x; C(1,:) = 1;
+%!   for k = 2:m; C(k+1,:) = 2*(k+mu-1)/k*x.*C(k,:) - (k+2*mu-2)/k*C(k-1,:); end
+%!   b = real(ifft(C(m+1,:))); b = b(m/2+2:L)/b(1);
+%!   assert(ultrwin(L, mu, xmu, 'x')', [b 1 fliplr(b)], 1e-12);
+%! end
+
+%!test
+%! b = [
+%!   5.7962919401511820e-03
+%!   1.6086991349967078e-02
+%!   3.6019014684117417e-02
+%!   6.8897525451558125e-02
+%!   1.1802364384553447e-01
+%!   1.8566749737411145e-01
+%!   2.7234740630826737e-01
+%!   3.7625460141456091e-01
+%!   4.9297108901880221e-01
+%!   6.1558961695849457e-01
+%!   7.3527571856983598e-01
+%!   8.4222550739092694e-01
+%!   9.2688779484512085e-01
+%!   9.8125497127708561e-01]';
+%! [w xmu] = ultrwin(29, 0, 3);
+%! assert(w', [b 1 fliplr(b)], 1e-14);
+%! assert(xmu, 1.053578297819277, 1e-14);
+
+%!test
+%! b = [
+%!   2.9953636903962466e-02
+%!   7.6096450051659603e-02
+%!   1.5207129867916891e-01
+%!   2.5906995366355179e-01
+%!   3.9341065451220536e-01
+%!   5.4533014012036929e-01
+%!   6.9975915071207051e-01
+%!   8.3851052636906720e-01
+%!   9.4345733548690369e-01]';
+%! assert(ultrwin(20, .5, 50, 'a')', [b 1 1 fliplr(b)], 1e-14);
+
+%!test
+%! b = [
+%!   1.0159906492322712e-01
+%!   1.4456358609406283e-01
+%!   2.4781689516201011e-01
+%!   3.7237015168857646e-01
+%!   5.1296973026690407e-01
+%!   6.5799041448113671e-01
+%!   7.9299087042967320e-01
+%!   9.0299778924260576e-01
+%!   9.7496213649820296e-01]';
+%! assert(ultrwin(19, -.4, 40, 'l')', [b 1 fliplr(b)], 1e-14);
+
+%!demo
+%! w=ultrwin(120, -1, 40, 'l'); [W,f]=freqz(w); clf
+%! subplot(2,1,1); plot(f/pi, 20*log10(W/abs(W(1)))); grid; axis([0 1 -90 0])
+%! subplot(2,1,2); plot(0:length(w)-1, w); grid
+%! %-----------------------------------------------------------
+%! % Figure shows an Ultraspherical window with MU=-1, LATT=40:
+%! % frequency domain above, time domain below.
+
+%!demo
+%! c='krbm'; clf; subplot(2, 1, 1)
+%! for beta=2:5
+%!   w=ultrwin(80, -.5, beta); [W,f]=freqz(w);
+%!   plot(f/pi, 20*log10(W/abs(W(1))), c(1+mod(beta, length(c)))); hold on
+%! end; grid; axis([0 1 -140 0]); hold off
+%! subplot(2, 1, 2);
+%! for n=2:10
+%!   w=ultrwin(n*20, 1, 3); [W,f]=freqz(w,1,2^11);
+%!   plot(f/pi, 20*log10(W/abs(W(1))), c(1+mod(n, length(c)))); hold on
+%! end; grid; axis([0 .2 -100 0]); hold off
+%! %--------------------------------------------------
+%! % Figure shows transfers of Ultraspherical windows:
+%! % above: varying BETA with fixed N & MU,
+%! % below: varying N with fixed MU & BETA.
+
+%!demo
+%! c='krbm'; clf; subplot(2, 1, 1)
+%! for j=0:4
+%!   w=ultrwin(80, j*.6-1.2, 50, 'a'); [W,f]=freqz(w);
+%!   plot(f/pi, 20*log10(W/abs(W(1))), c(1+mod(j, length(c)))); hold on
+%! end; grid; axis([0 1 -100 0]); hold off
+%! subplot(2, 1, 2);
+%! for j=4:-1:0
+%!   w=ultrwin(80, j*.75-1.5, 50, 'l'); [W,f]=freqz(w);
+%!   plot(f/pi, 20*log10(W/abs(W(1))), c(1+mod(j, length(c)))); hold on
+%! end; grid; axis([0 1 -100 0]); hold off
+%! %--------------------------------------------------
+%! % Figure shows transfers of Ultraspherical windows:
+%! % above: varying MU with fixed N & ATT,
+%! % below: varying MU with fixed N & LATT.
+
+%!demo
+%! clf; a=[.8 2 -115 5]; fc=1.1/pi; l='labelxy';
+%! for k=1:3; switch (k); case 1; w=kaiser(L=159, 7.91);
+%!   case 2; w=ultrwin(L=165, 0, 2.73); case 3; w=ultrwin(L=153, .5, 2.6); end
+%!   subplot(3, 1, 4-k); f=[1:(L-1)/2]*pi;f=sin(fc*f)./f; f=[fliplr(f) fc f]';
+%!   [h,f]=freqz(w.*f,1,2^14); plot(f,20*log10(h)); grid; axis(a,l); l='labely';
+%! end
+%! %-----------------------------------------------------------
+%! % Figure shows example lowpass filter design (Fp=1, Fs=1.2
+%! % rad/s, att=80 dB) and comparison with other windows.  From
+%! % top to bottom: Ultraspherical, Dolph-Chebyshev, and Kaiser
+%! % windows, with lengths 153, 165, and 159 respectively.
diff -r b6c2035313e6 src/Makefile
--- a/src/Makefile	Sun May 19 17:15:31 2013 +0200
+++ b/src/Makefile	Sat Jun 15 20:23:36 2013 +0100
@@ -1,6 +1,6 @@
 MKOCTFILE ?= mkoctfile
 
-all: __fwht__.oct cl2bp.oct remez.oct medfilt1.oct sosfilt.oct upfirdn.oct
+all: __fwht__.oct cl2bp.oct remez.oct medfilt1.oct sosfilt.oct upfirdn.oct ultrwin_.oct
 
 %.oct: %.cc
 	$(MKOCTFILE) -Wall -s $<
diff -r b6c2035313e6 src/ultrwin.c
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/src/ultrwin.c	Sat Jun 15 20:23:36 2013 +0100
@@ -0,0 +1,215 @@
+/* Copyright (c) 2013 Rob Sykes <robs@users.sourceforge.net>
+ *
+ * This program is free software; you can redistribute it and/or modify it
+ * under the terms of the GNU General Public License as published by the Free
+ * Software Foundation: either version 3 of the License, or (at your option)
+ * any later version.  See http://gnu.org/licenses for details.
+ *
+ * This program is distributed in the hope that it will be useful, but WITHOUT
+ * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+ * FITNESS FOR A PARTICULAR PURPOSE.  */
+
+#include <math.h>
+#include <stdlib.h>
+#include <stdio.h>
+
+#if !defined M_PI
+  #define M_PI 3.14159265358979323846
+#endif
+#if defined _MSC_VER && _MSC_VER < 1700
+  #define acosh(x) log((x) + sqrt(1. * (x) * (x) - 1))
+#endif
+
+#include "ultrwin.h"
+
+#if defined NDEBUG
+  #define fprintf if (0) fprintf
+#endif
+#define SHOW1(x) fprintf(stderr, "%c%+.3e", " 1"[(x) > .5], (x) - ((x) > .5))
+
+#define EPSILON     (1./0x100000) /* For caller rounding error. */
+#define BETA_MAX    (12*(1+EPSILON))
+#define XMU_MIN     (.99*(1-EPSILON))
+
+
+
+static double * ultraspherical_win(int n, double mu, double xmu)
+{
+  double * w = NULL;
+  int bad = n < 1 || xmu < XMU_MIN || (!mu && xmu == 1) ||
+    (n > BETA_MAX * 2 && xmu * cos(M_PI * BETA_MAX / n) > 1);
+  if (!bad && (w = (double *)malloc(sizeof(*w)* n))) {
+    int i, j, k, l = 0, m = (n + 1) / 2, met;
+    double * divs = w + m - 1, c = 1 - 1 / (xmu * xmu), t, u, v[64], vp, s;
+    for (i = 0; i < (int)(sizeof(v) / sizeof(v[0])); v[i++] = 0);
+    if (n > 1) for (i = 0; i < m; l = j - (j <= i++)) {
+      vp = *v, s = *v = i? (*v + v[1]) * mu * (divs[i] = 1./i) : 1;
+      for (met = 0, j = 1, u = 1; ; ++l, v[l] = vp * (i - l) / (mu + l - 1)) {
+        #define _ t = v[j], v[j] += vp, vp = t, t = s, s += \
+            v[j] * (u *= c * (n - i - j) * divs[j]), met = s && s == t, ++j,
+        for (k = ((l-j+1) & ~7) + j; j < k && !met; _ _ _ _ _ _ _ _ (void)0);
+        for (; j <= l && !met; _ (void)0);
+        #undef _
+        if (met || !(j <= i)) break;
+      }
+      w[i] = s / (n - i - 1);
+    }
+    else w[0] = 1;
+    u = 1 / w[i = m - 1], w[i] = 1;
+    for (--i ; i >= 0; u *= (n - 2 - i + mu) / (n - 2 - i), w[i] *= u, --i);
+    for (i = 0; i < m; w[n - 1 - i] = w[i], ++i);
+  }
+  return w;
+}
+
+
+
+typedef struct {double f, fp;} uspv_t;
+
+static uspv_t ultraspherical_polyval(
+    int n, double mu, double x, double const * divs)
+{
+  double fp = n > 0? 2 * x * mu : 1, fpp = 1, f;
+  uspv_t result;
+  int i, k;
+  #define _ f = (2*x*(i+mu)*fp - (i+2*mu-1)*fpp) * divs[i+1], fpp=fp, fp=f, ++i,
+  for (i = 1, k = i + ((n - i) & ~7); i < k; _ _ _ _ _ _ _ _ (void)0);
+  for (; i < n; _ (void)0);
+  #undef _
+  result.f = fp, result.fp = fpp;
+  return result;
+}
+
+
+
+#define MU_EPSILON  (1./0x4000)
+#define EQ(mu,x)    (fabs((mu)-(x)) < MU_EPSILON)
+
+static uspv_t ultraspherical_polyval2(     /* With non-+ve integer protection */
+    int n, double mu, double x, double const * divs)
+{
+  int sign = (~(int)floor(mu) & ~(~2u/2))? 1:-1; /* -ve if floor(mu) <0 & odd */
+  uspv_t r;
+  if (mu < MU_EPSILON && EQ(mu,(int)mu))
+    mu = floor(mu + .5) + MU_EPSILON * ((int)mu > mu? -1:1);
+  r = ultraspherical_polyval(n, mu, x, divs);
+  r.f *= sign, r.fp *= sign;
+  return r;
+}
+
+
+
+static double find_zero(int n, double mu, int l, double extremum_mag,
+    double ripple_ratio, double lower_bound, double const * divs)
+{
+  double dx, x0, t, x, epsilon = 1e-10, one_over_deriv, target = 0;
+  int i, met = 0;
+  if (!divs)
+    return 0;
+  if (!l) {
+    double r = ripple_ratio;    /* TODO: factor in weighted extremum_mag here */
+    x = r > 1 ? cosh(acosh(r) / n) : cos(acos(r) / n); /* invert chebpoly-1st */
+    x0 = x *= lower_bound / cos(M_PI * .5 / n) + epsilon;
+    target = log(extremum_mag * ripple_ratio);
+  }
+  else {
+    double cheb1 = cos(M_PI * (l - .5) / n), cheb2 = cos(M_PI * l / (n + 1));
+    if (mu < 1 - l && EQ((int)(mu+.5),mu+.5)) x = met = 1;
+    else if (EQ(mu,0)) x = cheb1, met = 1;               /* chebpoly-1st-kind */
+    else if (EQ(mu,1)) x = cheb2, met = 1;               /* chebpoly-2nd-kind */
+    else x = (cheb1 * cheb2) / (mu * cheb1 + (1 - mu) * cheb2);
+    x0 = x;
+  }
+  for (i = 0; i < 24 && !met; ++i, met = fabs(dx) < epsilon) {/*Newton-Raphson*/
+    uspv_t r = ultraspherical_polyval2(n, mu, x, divs);
+    if (!(t = ((2*mu + n-1) * r.fp - n*x * r.f)))         /* Fail if div by 0 */
+      break;
+    one_over_deriv = (1 - x*x) / t;    /* N-R slow for deriv~=1, so take log: */
+    if (!l) {                               /* d/dx(f(g(x))) = f'(g(x)).g'(x) */
+      one_over_deriv *= r.f;                             /* d/dx(log x) = 1/x */
+      if (r.f <= 0)                                 /* Fail if log of non-+ve */
+        break;
+      if (x + (dx = (target - log(r.f)) * one_over_deriv) <= lower_bound)
+        dx = (lower_bound - x) * .875;
+      x += dx;
+    }
+    else x += dx = -r.f * one_over_deriv;
+    fprintf(stderr, "1/deriv=%9.2e dx=%9.2e x=", one_over_deriv, dx);
+    SHOW1(x); fprintf(stderr, "\n");
+  }
+  fprintf(stderr, "find_zero(n=%i mu=%g l=%i target=%g r=%g x0=",
+      n, mu, l, target, ripple_ratio);
+  SHOW1(x0); fprintf(stderr, ") %s ", met? "converged to" : "FAILED at");
+  SHOW1(x); fprintf(stderr, " in %i iterations\n", i);
+  return met? x : 0;
+}
+
+
+static double * make_divs(int n, double * * divs)
+{
+  int i;
+  if (!*divs) {
+    *divs = (double *)malloc(n * sizeof(**divs));
+    if (*divs)
+      for (i = 0; i < n; (*divs)[i] = 1./(i+1), ++i);
+  }
+  return *divs? *divs - 1 : 0;
+}
+
+#define DIVS make_divs(n, &divs)
+
+
+
+double * ultraspherical_window(
+    int n, double mu, double par, uswpt_t type, int even_norm, double * xmu_)
+{
+  double * w = 0, xmu = 0, * divs = 0, last_extremum_pos = 0;
+
+  if (n > 0 && fabs(mu) <= (8*(1+EPSILON))) switch (type) {
+    case uswpt_Beta:
+      xmu = mu == 1 && par == 1? 1 : par < .5 || par > BETA_MAX? 0 :
+        find_zero(n-1, mu, 1, 0, 0, 0, DIVS) / cos(M_PI * par / n);
+      break;
+
+    case uswpt_AttFirst: if (par < 0) break; /* Falling... */
+
+    case uswpt_AttLast:
+      if (type == uswpt_AttLast && mu >= 0 && par < 0);
+      else if (!EQ(mu,0)) {
+        int extremum_num =
+          type == uswpt_AttLast? (int)((n-2)/2 +.5) : 1 + EQ(mu,-1.5);
+        double extremum_pos =
+          find_zero(n-2, mu+1, extremum_num, 0, 0, 0, DIVS);
+        double extremum_mag = !extremum_pos? 0 :
+          fabs(ultraspherical_polyval2(n-1, mu, extremum_pos, DIVS).f);
+        double xmu_lower_bound = !extremum_mag? 0 :
+          find_zero(n-1, mu, 1, 0, 0, 0, DIVS); /* 1st null */
+        xmu = !xmu_lower_bound? 0 : find_zero(
+            n-1, mu, 0, extremum_mag, pow(10, par/20), xmu_lower_bound, DIVS);
+        last_extremum_pos =
+          type == uswpt_AttLast? extremum_pos : last_extremum_pos;
+      }
+      else xmu = cosh(acosh(pow(10, par/20))/(n-1)); /* Cheby 1st kind */
+      break;
+
+    default: case uswpt_Xmu: xmu = par; break;
+  }
+  fprintf(stderr, "n=%i mu=%.3f xmu=%.16g\n", n, mu, xmu);
+
+  if (xmu > 0)
+    w = ultraspherical_win(n, mu, xmu);
+
+  if (w && (~n & !!even_norm) && n > 2 && !(mu == 1 && xmu == 1)) {
+    int i = n / 2 - 1, j = 1;
+    double * d = DIVS, t = 0, s = -1, x = even_norm == 1? 0 : last_extremum_pos?
+      last_extremum_pos : find_zero(n-2, mu+1, i, 0, 0, 0, d);
+    x = x? M_PI/2 - acos(x/xmu) : 0;
+    for (; i >= 0; t += w[i] * d[j] * (s=-s) * (x?cos(j*x):1), --i, j += 2);
+    for (t = M_PI/4 / t, i = 0; t < 1 && i < n; w[i] *= t, ++i);
+    fprintf(stderr, "%snorm DFT(w.sinc πx) @ %g %.16g\n", t<1? "":"NO ", 2*x,t);
+  }
+  free(divs);
+  if (xmu_)
+    *xmu_ = xmu;
+  return w;
+}
diff -r b6c2035313e6 src/ultrwin.h
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/src/ultrwin.h	Sat Jun 15 20:23:36 2013 +0100
@@ -0,0 +1,15 @@
+/* Copyright (c) 2013 Rob Sykes <robs@users.sourceforge.net>
+ *
+ * This program is free software; you can redistribute it and/or modify it
+ * under the terms of the GNU General Public License as published by the Free
+ * Software Foundation: either version 3 of the License, or (at your option)
+ * any later version.  See http://gnu.org/licenses for details.
+ *
+ * This program is distributed in the hope that it will be useful, but WITHOUT
+ * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+ * FITNESS FOR A PARTICULAR PURPOSE.  */
+
+typedef enum {uswpt_Xmu, uswpt_Beta, uswpt_AttFirst, uswpt_AttLast} uswpt_t;
+
+double * ultraspherical_window(
+  int n, double mu, double par, uswpt_t par_type, int even_norm, double * xmu_);
diff -r b6c2035313e6 src/ultrwin_.cc
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/src/ultrwin_.cc	Sat Jun 15 20:23:36 2013 +0100
@@ -0,0 +1,38 @@
+/* Copyright (c) 2013 Rob Sykes <robs@users.sourceforge.net>
+ *
+ * This program is free software; you can redistribute it and/or modify it
+ * under the terms of the GNU General Public License as published by the Free
+ * Software Foundation: either version 3 of the License, or (at your option)
+ * any later version.  See http://gnu.org/licenses for details.
+ *
+ * This program is distributed in the hope that it will be useful, but WITHOUT
+ * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+ * FITNESS FOR A PARTICULAR PURPOSE.  */
+
+#include <octave/oct.h>
+#define NDEBUG
+#include "ultrwin.c"
+
+#define USAGE "internal function usage: ultrwin_(L, MU, PAR, PAR_TYPE, NORM)"
+
+DEFUN_DLD (ultrwin_, args, , USAGE) {
+  for (int i = 0; i < 5; ++i)
+    if (args.length() <= i || !args(i).is_real_scalar()) {
+      error(USAGE);
+      return octave_value();
+    }
+  int L;
+  double xmu, * w = ultraspherical_window(L = args(0).scalar_value(),
+      args(1).scalar_value(), args(2).scalar_value(),
+      (uswpt_t)args(3).scalar_value(), args(4).scalar_value(), &xmu);
+  if (!w) {
+    error("ultrwin: parameter(s) out of range");
+    return octave_value();
+  }
+  octave_value_list result;
+  ColumnVector ww(L);
+  for (int i = 0; i < L; ww(i) = w[i], ++i);
+  free(w);
+  result(0) = ww, result(1) = xmu;
+  return result;
+}
