1 /*
2 ==============================================================================
3
4 This file is part of the JUCE library.
5 Copyright (c) 2017 - ROLI Ltd.
6
7 JUCE is an open source library subject to commercial or open-source
8 licensing.
9
10 By using JUCE, you agree to the terms of both the JUCE 5 End-User License
11 Agreement and JUCE 5 Privacy Policy (both updated and effective as of the
12 27th April 2017).
13
14 End User License Agreement: www.juce.com/juce-5-licence
15 Privacy Policy: www.juce.com/juce-5-privacy-policy
16
17 Or: You may also use this code under the terms of the GPL v3 (see
18 www.gnu.org/licenses).
19
20 JUCE IS PROVIDED "AS IS" WITHOUT ANY WARRANTY, AND ALL WARRANTIES, WHETHER
21 EXPRESSED OR IMPLIED, INCLUDING MERCHANTABILITY AND FITNESS FOR PURPOSE, ARE
22 DISCLAIMED.
23
24 ==============================================================================
25 */
26
27 namespace juce
28 {
29 namespace dsp
30 {
31
32 template <typename ElementType>
identity(size_t size)33 Matrix<ElementType> Matrix<ElementType>::identity (size_t size)
34 {
35 Matrix result (size, size);
36
37 for (size_t i = 0; i < size; ++i)
38 result(i, i) = 1;
39
40 return result;
41 }
42
43 template <typename ElementType>
toeplitz(const Matrix & vector,size_t size)44 Matrix<ElementType> Matrix<ElementType>::toeplitz (const Matrix& vector, size_t size)
45 {
46 jassert (vector.isOneColumnVector());
47 jassert (size <= vector.rows);
48
49 Matrix result (size, size);
50
51 for (size_t i = 0; i < size; ++i)
52 result (i, i) = vector (0, 0);
53
54 for (size_t i = 1; i < size; ++i)
55 for (size_t j = i; j < size; ++j)
56 result (j, j - i) = result (j - i, j) = vector (i, 0);
57
58 return result;
59 }
60
61 template <typename ElementType>
hankel(const Matrix & vector,size_t size,size_t offset)62 Matrix<ElementType> Matrix<ElementType>::hankel (const Matrix& vector, size_t size, size_t offset)
63 {
64 jassert(vector.isOneColumnVector());
65 jassert(vector.rows >= (2 * (size - 1) + 1));
66
67 Matrix result (size, size);
68
69 for (size_t i = 0; i < size; ++i)
70 result (i, i) = vector ((2 * i) + offset, 0);
71
72 for (size_t i = 1; i < size; ++i)
73 for (size_t j = i; j < size; ++j)
74 result (j, j - i) = result (j - i, j) = vector (i + 2 * (j - i) + offset, 0);
75
76 return result;
77 }
78
79 //==============================================================================
80 template <typename ElementType>
swapColumns(size_t columnOne,size_t columnTwo)81 Matrix<ElementType>& Matrix<ElementType>::swapColumns (size_t columnOne, size_t columnTwo) noexcept
82 {
83 jassert (columnOne < columns && columnTwo < columns);
84
85 auto* p = data.getRawDataPointer();
86
87 for (size_t i = 0; i < rows; ++i)
88 {
89 auto offset = dataAcceleration.getUnchecked (static_cast<int> (i));
90 std::swap (p[offset + columnOne], p[offset + columnTwo]);
91 }
92
93 return *this;
94 }
95
96 template <typename ElementType>
swapRows(size_t rowOne,size_t rowTwo)97 Matrix<ElementType>& Matrix<ElementType>::swapRows (size_t rowOne, size_t rowTwo) noexcept
98 {
99 jassert (rowOne < rows && rowTwo < rows);
100
101 auto offset1 = rowOne * columns;
102 auto offset2 = rowTwo * columns;
103
104 auto* p = data.getRawDataPointer();
105
106 for (size_t i = 0; i < columns; ++i)
107 std::swap (p[offset1 + i], p[offset2 + i]);
108
109 return *this;
110 }
111
112 //==============================================================================
113 template <typename ElementType>
operator *(const Matrix<ElementType> & other) const114 Matrix<ElementType> Matrix<ElementType>::operator* (const Matrix<ElementType>& other) const
115 {
116 auto n = getNumRows(), m = other.getNumColumns(), p = getNumColumns();
117 Matrix result (n, m);
118
119 jassert (p == other.getNumRows());
120
121 size_t offsetMat = 0, offsetlhs = 0;
122
123 auto* dst = result.getRawDataPointer();
124 auto* a = getRawDataPointer();
125 auto* b = other.getRawDataPointer();
126
127 for (size_t i = 0; i < n; ++i)
128 {
129 size_t offsetrhs = 0;
130
131 for (size_t k = 0; k < p; ++k)
132 {
133 auto ak = a[offsetlhs++];
134
135 for (size_t j = 0; j < m; ++j)
136 dst[offsetMat + j] += ak * b[offsetrhs + j];
137
138 offsetrhs += m;
139 }
140
141 offsetMat += m;
142 }
143
144 return result;
145 }
146
147 //==============================================================================
148 template <typename ElementType>
compare(const Matrix & a,const Matrix & b,ElementType tolerance)149 bool Matrix<ElementType>::compare (const Matrix& a, const Matrix& b, ElementType tolerance) noexcept
150 {
151 if (a.rows != b.rows || a.columns != b.columns)
152 return false;
153
154 tolerance = std::abs (tolerance);
155
156 auto* bPtr = b.begin();
157 for (auto aValue : a)
158 if (std::abs (aValue - *bPtr++) > tolerance)
159 return false;
160
161 return true;
162 }
163
164 //==============================================================================
165 template <typename ElementType>
solve(Matrix & b) const166 bool Matrix<ElementType>::solve (Matrix& b) const noexcept
167 {
168 auto n = columns;
169 jassert (n == n && n == b.rows && b.isOneColumnVector());
170
171 auto* x = b.getRawDataPointer();
172 const auto& A = *this;
173
174 switch (n)
175 {
176 case 1:
177 {
178 auto denominator = A (0,0);
179
180 if (denominator == 0)
181 return false;
182
183 b (0, 0) /= denominator;
184 }
185 break;
186
187 case 2:
188 {
189 auto denominator = A (0, 0) * A (1, 1) - A (0, 1) * A (1, 0);
190
191 if (denominator == 0)
192 return false;
193
194 auto factor = (1 / denominator);
195 auto b0 = x[0], b1 = x[1];
196
197 x[0] = factor * (A (1, 1) * b0 - A (0, 1) * b1);
198 x[1] = factor * (A (0, 0) * b1 - A (1, 0) * b0);
199 }
200 break;
201
202 case 3:
203 {
204 auto denominator = A (0, 0) * (A (1, 1) * A (2, 2) - A (1, 2) * A (2, 1))
205 + A (0, 1) * (A (1, 2) * A (2, 0) - A (1, 0) * A (2, 2))
206 + A (0, 2) * (A (1, 0) * A (2, 1) - A (1, 1) * A (2, 0));
207
208 if (denominator == 0)
209 return false;
210
211 auto factor = 1 / denominator;
212 auto b0 = x[0], b1 = x[1], b2 = x[2];
213
214 x[0] = ( ( A (0, 1) * A (1, 2) - A (0, 2) * A (1, 1)) * b2
215 + (-A (0, 1) * A (2, 2) + A (0, 2) * A (2, 1)) * b1
216 + ( A (1, 1) * A (2, 2) - A (1, 2) * A (2, 1)) * b0) * factor;
217
218 x[1] = -( ( A (0, 0) * A (1, 2) - A (0, 2) * A (1, 0)) * b2
219 + (-A (0, 0) * A (2, 2) + A (0, 2) * A (2, 0)) * b1
220 + ( A (1, 0) * A (2, 2) - A (1, 2) * A (2, 0)) * b0) * factor;
221
222 x[2] = ( ( A (0, 0) * A (1, 1) - A (0, 1) * A (1, 0)) * b2
223 + (-A (0, 0) * A (2, 1) + A (0, 1) * A (2, 0)) * b1
224 + ( A (1, 0) * A (2, 1) - A (1, 1) * A (2, 0)) * b0) * factor;
225 }
226 break;
227
228
229 default:
230 {
231 Matrix<ElementType> M (A);
232
233 for (size_t j = 0; j < n; ++j)
234 {
235 if (M (j, j) == 0)
236 {
237 auto i = j;
238 while (i < n && M (i, j) == 0)
239 ++i;
240
241 if (i == n)
242 return false;
243
244 for (size_t k = 0; k < n; ++k)
245 M (j, k) += M (i, k);
246
247 x[j] += x[i];
248 }
249
250 auto t = 1 / M (j, j);
251
252 for (size_t k = 0; k < n; ++k)
253 M (j, k) *= t;
254
255 x[j] *= t;
256
257 for (size_t k = j + 1; k < n; ++k)
258 {
259 auto u = -M (k, j);
260
261 for (size_t l = 0; l < n; ++l)
262 M (k, l) += u * M (j, l);
263
264 x[k] += u * x[j];
265 }
266 }
267
268 for (int k = static_cast<int> (n) - 2; k >= 0; --k)
269 for (size_t i = static_cast<size_t> (k) + 1; i < n; ++i)
270 x[k] -= M (static_cast<size_t> (k), i) * x[i];
271 }
272 }
273
274 return true;
275 }
276
277 //==============================================================================
278 template <typename ElementType>
toString() const279 String Matrix<ElementType>::toString() const
280 {
281 StringArray entries;
282 int sizeMax = 0;
283
284 auto* p = data.begin();
285
286 for (size_t i = 0; i < rows; ++i)
287 {
288 for (size_t j = 0; j < columns; ++j)
289 {
290 String entry (*p++, 4);
291 sizeMax = jmax (sizeMax, entry.length());
292
293 entries.add (entry);
294 }
295 }
296
297 sizeMax = ((sizeMax + 1) / 4 + 1) * 4;
298
299 MemoryOutputStream result;
300
301 auto n = static_cast<size_t> (entries.size());
302
303 for (size_t i = 0; i < n; ++i)
304 {
305 result << entries[(int) i].paddedRight (' ', sizeMax);
306
307 if (i % columns == (columns - 1))
308 result << newLine;
309 }
310
311 return result.toString();
312 }
313
314 template class Matrix<float>;
315 template class Matrix<double>;
316
317 } // namespace dsp
318 } // namespace juce
319