--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/toolbox/pdtoep.m	Thu May 07 18:36:24 2015 +0200
@@ -0,0 +1,29 @@
+function T = pdtoep(n, m, w, theta)
+%PDTOEP   Symmetric positive definite Toeplitz matrix.
+%         PDTOEP(N, M, W, THETA) is an N-by-N symmetric positive (semi-)
+%         definite (SPD) Toeplitz matrix, comprised of the sum of M rank 2
+%         (or, for certain THETA, rank 1) SPD Toeplitz matrices.
+%         Specifically,
+%                 T = W(1)*T(THETA(1)) + ... + W(M)*T(THETA(M)),
+%         where T(THETA(k)) has (i,j) element COS(2*PI*THETA(k)*(i-j)).
+%         Defaults: M = N, W = RAND(M,1), THETA = RAND(M,1).
+
+%         Reference:
+%         G. Cybenko and C.F. Van Loan, Computing the minimum eigenvalue of
+%         a symmetric positive definite Toeplitz matrix, SIAM J. Sci. Stat.
+%         Comput., 7 (1986), pp. 123-131.
+
+if nargin < 2, m = n; end
+if nargin < 3, w = rand(m,1); end
+if nargin < 4, theta = rand(m,1); end
+
+if max(size(w)) ~= m | max(size(theta)) ~= m
+   error('Arguments W and THETA must be vectors of length M.')
+end
+
+T = zeros(n);
+E = 2*pi*( (1:n)'*ones(1,n) - ones(n,1)*(1:n) );
+
+for i=1:m
+    T = T + w(i) * cos( theta(i)*E );
+end