comparison toolbox/orthog.m @ 2:c124219d7bfa draft

Re-add the 1995 toolbox after noticing the statement in the ~higham/mctoolbox/ webpage.
author Antonio Pino Robles <data.script93@gmail.com>
date Thu, 07 May 2015 18:36:24 +0200
parents 8f23314345f4
children
comparison
equal deleted inserted replaced
1:e471a92d17be 2:c124219d7bfa
1 function Q = orthog(n, k)
2 %ORTHOG Orthogonal and nearly orthogonal matrices.
3 % Q = ORTHOG(N, K) selects the K'th type of matrix of order N.
4 % K > 0 for exactly orthogonal matrices, K < 0 for diagonal scalings of
5 % orthogonal matrices.
6 % Available types: (K = 1 is the default)
7 % K = 1: Q(i,j) = SQRT(2/(n+1)) * SIN( i*j*PI/(n+1) )
8 % Symmetric eigenvector matrix for second difference matrix.
9 % K = 2: Q(i,j) = 2/SQRT(2*n+1)) * SIN( 2*i*j*PI/(2*n+1) )
10 % Symmetric.
11 % K = 3: Q(r,s) = EXP(2*PI*i*(r-1)*(s-1)/n) / SQRT(n) (i=SQRT(-1))
12 % Unitary, the Fourier matrix. Q^4 is the identity.
13 % This is essentially the same matrix as FFT(EYE(N))/SQRT(N)!
14 % K = 4: Helmert matrix: a permutation of a lower Hessenberg matrix,
15 % whose first row is ONES(1:N)/SQRT(N).
16 % K = 5: Q(i,j) = SIN( 2*PI*(i-1)*(j-1)/n ) + COS( 2*PI*(i-1)*(j-1)/n ).
17 % Symmetric matrix arising in the Hartley transform.
18 % K = -1: Q(i,j) = COS( (i-1)*(j-1)*PI/(n-1) )
19 % Chebyshev Vandermonde-like matrix, based on extrema of T(n-1).
20 % K = -2: Q(i,j) = COS( (i-1)*(j-1/2)*PI/n) )
21 % Chebyshev Vandermonde-like matrix, based on zeros of T(n).
22
23 % References:
24 % N.J. Higham and D.J. Higham, Large growth factors in Gaussian
25 % elimination with pivoting, SIAM J. Matrix Analysis and Appl.,
26 % 10 (1989), pp. 155-164.
27 % P. Morton, On the eigenvectors of Schur's matrix, J. Number Theory,
28 % 12 (1980), pp. 122-127. (Re. ORTHOG(N, 3))
29 % H.O. Lancaster, The Helmert Matrices, Amer. Math. Monthly, 72 (1965),
30 % pp. 4-12.
31 % D. Bini and P. Favati, On a matrix algebra related to the discrete
32 % Hartley transform, SIAM J. Matrix Anal. Appl., 14 (1993),
33 % pp. 500-507.
34
35 if nargin == 1, k = 1; end
36
37 if k == 1
38 % E'vectors second difference matrix
39 m = (1:n)'*(1:n) * (pi/(n+1));
40 Q = sin(m) * sqrt(2/(n+1));
41
42 elseif k == 2
43
44 m = (1:n)'*(1:n) * (2*pi/(2*n+1));
45 Q = sin(m) * (2/sqrt(2*n+1));
46
47 elseif k == 3 % Vandermonde based on roots of unity
48
49 m = 0:n-1;
50 Q = exp(m'*m*2*pi*sqrt(-1)/n) / sqrt(n);
51
52 elseif k == 4 % Helmert matrix
53
54 Q = tril(ones(n));
55 Q(1,2:n) = ones(1,n-1);
56 for i=2:n
57 Q(i,i) = -(i-1);
58 end
59 Q = diag( sqrt( [n 1:n-1] .* [1:n] ) ) \ Q;
60
61 elseif k == 5 % Hartley matrix
62
63 m = (0:n-1)'*(0:n-1) * (2*pi/n);
64 Q = (cos(m) + sin(m))/sqrt(n);
65
66 elseif k == -1
67 % extrema of T(n-1)
68 m = (0:n-1)'*(0:n-1) * (pi/(n-1));
69 Q = cos(m);
70
71 elseif k == -2
72 % zeros of T(n)
73 m = (0:n-1)'*(.5:n-.5) * (pi/n);
74 Q = cos(m);
75
76 else
77
78 error('Illegal value for second parameter.')
79
80 end