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