1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
|
/* For information on usage and redistribution, and for a DISCLAIMER OF ALL
* WARRANTIES, see the file, "LICENSE.txt," in this distribution.
iem_matrix written by Thomas Musil (c) IEM KUG Graz Austria 2002 - 2006 */
/* -------------------------- matrix_orthogonal ------------------------------ */
/*
versucht, eine eingehende matrix (quadratisch) zu orthogonalisiern (Transponierte = Inverse)
und zwar wie folgt:
1. zeile orthonarmalisieren (summe aller zeilen-elemente-quadrate = 1)
K = summe aller a_1i (1 <= i <= N)
fuer alle an_1i = a_1i / sqrt(K) (1 <= i <= N)
2. zeile:
K = summe aller a_2i (2 <= i <= N) - an_12
an_21 = an_12
fuer alle an_2i = a_2i / sqrt(K) (2 <= i <= N)
3. zeile:
K = summe aller a_3i (3 <= i <= N) - an_13 - an_23
an_31 = an_13
an_32 = an_23
fuer alle an_3i = a_3i / sqrt(K) (3 <= i <= N)
usw.
*/
typedef struct _matrix_orthogonal
{
t_object x_obj;
int x_dim;
double *x_work;
t_atom *x_at;
} t_matrix_orthogonal;
static t_class *matrix_orthogonal_class;
static void matrix_orthogonal_matrix(t_matrix_orthogonal *x, t_symbol *s, int argc, t_atom *argv)
{
int oldsize = x->x_dim;
int dim;
int i, j;
int r, c;
double sum1, sum2, el, *v, *u;
t_atom *at=x->x_at;
if(argc > 2)
{
r = (int)(atom_getint(argv++));
c = (int)(atom_getint(argv++));
if(r == c)
{
dim = r;
if(dim < 1)
dim = 1;
if(dim > oldsize)
{
if(oldsize)
{
x->x_work = (double *)resizebytes(x->x_work, oldsize * oldsize * sizeof(double), dim * dim * sizeof(double));
x->x_at = (t_atom *)resizebytes(x->x_at, (oldsize * oldsize + 2) * sizeof(t_atom), (dim * dim + 2) * sizeof(t_atom));
x->x_dim = dim;
}
else
{
x->x_work = (double *)getbytes(dim * dim * sizeof(double));
x->x_at = (t_atom *)getbytes((dim * dim + 2) * sizeof(t_atom));
x->x_dim = dim;
}
}
v = x->x_work;
for(i=0; i<dim; i++) /* init */
{
for(j=0; j<dim; j++)
{
*v++ = (double)(atom_getfloat(argv++));
}
}
for(i=0; i<dim; i++) /* jede zeile */
{
sum1 = 0.0;
v = x->x_work + i*dim + i;
for(j=i; j<dim; j++) /* rest der zeile ab hauptdiagonale quadratisch summieren */
{
el = *v++;
sum1 += el * el;
}
v = x->x_work + i;
u = x->x_work + i*dim;
sum2 = 1.0;
for(j=0; j<i; j++)
{
el = *v;
v += dim;
*u++ = el;
sum2 -= el * el;
}
if(sum1 == 0.0)
sum1 = 1.0;
if(sum2 <= 0.0)
el = 0.0;
else
el = sqrt(sum2 / sum1);
for(j=i; j<dim; j++)
{
*u *= el;
u++;
}
}
at = x->x_at;
SETFLOAT(at, (t_float)dim);
at++;
SETFLOAT(at, (t_float)dim);
at++;
v = x->x_work;
for(i=0; i<dim; i++)
{
for(j=0; j<dim; j++)
{
SETFLOAT(at, (t_float)(*v));
at++;
v++;
}
}
outlet_anything(x->x_obj.ob_outlet, gensym("matrix"), dim*dim+2, x->x_at);
}
}
}
static void matrix_orthogonal_free(t_matrix_orthogonal *x)
{
if(x->x_dim)
{
freebytes(x->x_work, x->x_dim * x->x_dim * sizeof(double));
freebytes(x->x_at, (x->x_dim * x->x_dim + 2) * sizeof(t_atom));
}
}
static void *matrix_orthogonal_new(void)
{
t_matrix_orthogonal *x = (t_matrix_orthogonal *)pd_new(matrix_orthogonal_class);
x->x_dim = 0;
x->x_work = (double *)0;
x->x_at = (t_atom *)0;
outlet_new(&x->x_obj, &s_list);
return (x);
}
static void matrix_orthogonal_setup(void)
{
matrix_orthogonal_class = class_new(gensym("matrix_orthogonal"), (t_newmethod)matrix_orthogonal_new, (t_method)matrix_orthogonal_free,
sizeof(t_matrix_orthogonal), 0, 0);
class_addmethod(matrix_orthogonal_class, (t_method)matrix_orthogonal_matrix, gensym("matrix"), A_GIMME, 0);
class_sethelpsymbol(matrix_orthogonal_class, gensym("iemhelp/matrix_orthogonal-help"));
}
|
