1# test_SVD.praat
2# djmw 20171208
3
4appendInfoLine: "test_SVD.praat"
5
6n = 4
7eps = 1.5e-5
8appendInfoLine: tab$, "Hilbert matrix order 4"
9t = Create TableOfReal: "hilbert", n, n
10for i to n
11	for j to n
12		Set value: i, j, 1 / (i + j - 1)
13	endfor
14endfor
15svd = To SVD
16
17utab = Extract left singular vectors
18selectObject: svd
19vtab = Extract right singular vectors
20selectObject: svd
21dvec = Extract singular values
22u$= " -0.792608   0.582076  -0.179186  -0.029193 "
23... + " -0.451923  -0.370502   0.741918   0.328712 "
24... + " -0.322416  -0.509579  -0.100228  -0.791411 "
25... + " -0.252161  -0.514048  -0.638283   0.514553 "
26v$ = "-0.792608   0.582076  -0.179186  -0.029193 "
27... + " -0.451923  -0.370502   0.741918   0.328712 "
28... + " -0.322416  -0.509579  -0.100228  -0.791411 "
29... + " -0.252161  -0.514048  -0.638283   0.514553"
30d# = {1.5002, 1.6914e-01, 6.7383e-03, 9.6702e-05}
31
32appendInfoLine: tab$, tab$, "UDV' test equality of U"
33@string_to_table: u$, n, n
34utab_octave = selected ("TableOfReal")
35appendInfoLine: tab$, tab$, "UDV' test equality of V"
36@string_to_table: v$, n, n
37vtab_octave = selected ("TableOfReal")
38
39@check_tors: utab, utab_octave, eps
40@check_tors: vtab, vtab_octave, eps
41
42removeObject: utab, utab_octave, vtab, vtab_octave, svd, t, dvec
43
44appendInfoLine: tab$, "reconstruct 6x2 matrix"
45@test_reconstruction: 6, 2, eps
46
47appendInfoLine: tab$, "reconstruct 2x6 matrix"
48@test_reconstruction: 2, 6, eps
49
50appendInfoLine: tab$, "reconstruct 30x500 matrix"
51@test_reconstruction: 30, 500, eps
52
53appendInfoLine: "test_SVD.praat OK"
54
55procedure test_reconstruction: .nrows, .ncols, .eps
56	.t = Create TableOfReal: "t", .nrows, .ncols
57	Formula: "randomUniform (-1,1)"
58	.svd = To SVD
59	.tr = To TableOfReal: 1, 0
60	@check_tors: .t, .tr, .eps
61	removeObject: .t, .svd, .tr
62endproc
63
64procedure check_tors: .tor1, .tor2, .eps
65	selectObject: .tor1
66	.nrows = Get number of rows
67	.ncols = Get number of columns
68	for .i to .nrows
69		for .j to .ncols
70			.d =abs ((object[.tor1, .i, .j]-object[.tor2, .i, .j])/object[.tor2, .i, .j])
71			assert .d < .eps; ['.i','.j']: '.tor1', '.tor2', '.d'
72		endfor
73	endfor
74endproc
75
76procedure string_to_table: .string$, .nrows, .ncols
77	.s = Create Strings as tokens: .string$
78	.numberOfTokens = Get number of strings
79	assert .numberOfTokens == .nrows * .ncols
80	.t = Create TableOfReal: "t", .nrows, .ncols
81	for .i to .nrows
82		for .j to .ncols
83			selectObject: .s
84			.val$ = Get string: (.i-1)*.ncols + .j
85			selectObject: .t
86			Set value: .i, .j, number (.val$)
87		endfor
88	endfor
89	removeObject: .s
90	selectObject: .t
91endproc
92
93
94