Line data Source code
1 : /*
2 : * SPDX-License-Identifier: GPL-3.0-only
3 : * MuseScore-CLA-applies
4 : *
5 : * MuseScore Studio
6 : * Music Composition & Notation
7 : *
8 : * Copyright (C) 2023 MuseScore Limited and others
9 : *
10 : * This program is free software: you can redistribute it and/or modify
11 : * it under the terms of the GNU General Public License version 3 as
12 : * published by the Free Software Foundation.
13 : *
14 : * This program is distributed in the hope that it will be useful,
15 : * but WITHOUT ANY WARRANTY; without even the implied warranty of
16 : * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
17 : * GNU General Public License for more details.
18 : *
19 : * You should have received a copy of the GNU General Public License
20 : * along with this program. If not, see <https://www.gnu.org/licenses/>.
21 : */
22 :
23 : #ifndef MUSE_AUDIO_REVERBMATRICES_H
24 : #define MUSE_AUDIO_REVERBMATRICES_H
25 :
26 : #include "log.h"
27 :
28 : namespace muse::audio::fx {
29 : namespace reverb_matrices {
30 : template<int n>
31 : inline void Hadamard(const float i[n], float o[n])
32 : {
33 : UNUSED(i);
34 : UNUSED(o);
35 :
36 : // Hadamard matrix not defined for this size
37 : assert(false);
38 : }
39 :
40 : template<>
41 0 : inline void Hadamard<8>(const float i[8], float o[8])
42 : {
43 : // log2(n)*n adds n muls = 24 adds, 8 muls
44 :
45 0 : const float fact = 1.f / std::sqrt(8.f);
46 :
47 0 : float a[] = { i[7] + i[0], i[6] + i[1], i[5] + i[3], i[2] + i[4] };
48 0 : float b[] = { i[7] - i[0], i[6] - i[1], i[5] - i[3], i[2] - i[4] };
49 :
50 0 : float c[] = { a[0] + a[1], a[2] + a[3], b[0] + b[1], b[2] + b[3] };
51 0 : float d[] = { a[0] - a[1], a[2] - a[3], b[0] - b[1], b[2] - b[3] };
52 :
53 0 : o[0] = fact * (c[0] + c[1]);
54 0 : o[1] = fact * (c[2] + c[3]);
55 0 : o[2] = fact * (d[0] + d[1]);
56 0 : o[3] = fact * (d[2] + d[3]);
57 0 : o[4] = fact * (c[0] - c[1]);
58 0 : o[5] = fact * (c[2] - c[3]);
59 0 : o[6] = fact * (d[0] - d[1]);
60 0 : o[7] = fact * (d[2] - d[3]);
61 0 : }
62 :
63 : template<>
64 0 : inline void Hadamard<12>(const float i[12], float o[12])
65 : {
66 : // adds 12 + 24 + 24 = 60 adds, 12 muls
67 :
68 : // this might seem a bit crazy and it took me 4 attempts to get all signs and numbers right.. I don't have enough
69 : // autism. The 12x12 hadamard matrix is quite irregular (there exists only one).
70 : // I grouped operations as simd-friendly as I could.
71 :
72 0 : const float fact = 1.f / std::sqrt(12.f);
73 :
74 0 : float a[] = { i[7] + i[4], i[10] + i[11], i[2] + i[5], i[0] + i[1], i[8] + i[9], i[6] + i[3] };
75 0 : float b[] = { i[7] - i[4], i[10] - i[11], i[2] - i[5], i[0] - i[1], i[8] - i[9], i[6] - i[3] };
76 :
77 0 : float c[] = { a[0] + a[1], a[0] + b[1], b[0] + a[1], b[0] + b[1], a[0] - a[1], a[0] - b[1], b[0] - a[1], b[0] - b[1] };
78 0 : float d[] = { a[2] + a[3], a[2] + b[3], b[2] + a[3], b[2] + b[3], a[2] - a[3], a[2] - b[3], b[2] - a[3], b[2] - b[3] };
79 0 : float e[] = { a[4] + a[5], a[4] + b[5], b[4] + a[5], b[4] + b[5], a[4] - a[5], a[4] - b[5], b[4] - a[5], b[4] - b[5] };
80 :
81 0 : o[0] = fact * (c[0] + d[0] + e[0]);
82 0 : o[1] = fact * (c[3] + d[1] - e[5]);
83 0 : o[2] = fact * (c[7] - d[6] - e[1]);
84 0 : o[3] = fact * (c[4] + d[2] + e[6]);
85 0 : o[4] = fact * (c[3] - d[3] + e[4]);
86 0 : o[5] = fact * (c[7] - d[5] + e[1]);
87 0 : o[6] = fact * (c[6] + d[7] - e[6]);
88 0 : o[7] = fact * (c[4] - d[2] + e[2]);
89 0 : o[8] = fact * (c[1] - d[5] - e[3]);
90 0 : o[9] = fact * (c[0] - d[1] - e[5]);
91 0 : o[10] = fact * (c[2] + d[6] + e[7]);
92 0 : o[11] = fact * (c[5] + d[4] - e[2]);
93 0 : }
94 :
95 : template<>
96 0 : inline void Hadamard<16>(const float i[16], float o[16])
97 : {
98 : // 16 + 16 + 16 + 16 = 64 adds, 16 muls
99 :
100 0 : const float fact = 1.f / std::sqrt(16.f);
101 :
102 0 : float a[] = { i[10] + i[4], i[0] + i[5], i[15] + i[11], i[12] + i[7],
103 0 : i[2] + i[14], i[1] + i[8], i[3] + i[9], i[13] + i[6] };
104 0 : float b[] = { i[10] - i[4], i[0] - i[5], i[15] - i[11], i[12] - i[7],
105 0 : i[2] - i[14], i[1] - i[8], i[3] - i[9], i[13] - i[6] };
106 :
107 0 : float c[] = { a[0] + a[1], a[2] + a[3], a[4] + a[5], a[6] + a[7], b[0] + b[1], b[2] + b[3], b[4] + b[5], b[6] + b[7] };
108 0 : float d[] = { a[0] - a[1], a[2] - a[3], a[4] - a[5], a[6] - a[7], b[0] - b[1], b[2] - b[3], b[4] - b[5], b[6] - b[7] };
109 :
110 0 : float e[] = { c[0] + c[1], c[2] + c[3], d[0] + d[1], d[2] + d[3], c[4] + c[5], c[6] + c[7], d[4] + d[5], d[6] + d[7] };
111 0 : float f[] = { c[0] - c[1], c[2] - c[3], d[0] - d[1], d[2] - d[3], c[4] - c[5], c[6] - c[7], d[4] - d[5], d[6] - d[7] };
112 :
113 0 : o[0] = fact * (e[0] + e[1]);
114 0 : o[1] = fact * (e[2] + e[3]);
115 0 : o[2] = fact * (e[4] + e[5]);
116 0 : o[3] = fact * (e[6] + e[7]);
117 0 : o[4] = fact * (f[0] + f[1]);
118 0 : o[5] = fact * (f[2] + f[3]);
119 0 : o[6] = fact * (f[4] + f[5]);
120 0 : o[7] = fact * (f[6] + f[7]);
121 0 : o[8] = fact * (e[0] - e[1]);
122 0 : o[9] = fact * (e[2] - e[3]);
123 0 : o[10] = fact * (e[4] - e[5]);
124 0 : o[11] = fact * (e[6] - e[7]);
125 0 : o[12] = fact * (f[0] - f[1]);
126 0 : o[13] = fact * (f[2] - f[3]);
127 0 : o[14] = fact * (f[4] - f[5]);
128 0 : o[15] = fact * (f[6] - f[7]);
129 0 : }
130 :
131 : template<>
132 0 : inline void Hadamard<24>(const float i[24], float o[24])
133 : {
134 : // 24 + 48 + 48 = 120 adds, 24 muls
135 0 : const float fact = 1.f / std::sqrt(24.f);
136 :
137 0 : float a[] = { i[11] + i[12], i[23] + i[14], i[4] + i[2], i[1] + i[22], i[3] + i[19], i[9] + i[0],
138 0 : i[17] + i[7], i[5] + i[10], i[21] + i[13], i[6] + i[18], i[8] + i[15], i[20] + i[16] };
139 0 : float b[] = { i[11] - i[12], i[23] - i[14], i[4] - i[2], i[1] - i[22], i[3] - i[19], i[9] - i[0],
140 0 : i[17] - i[7], i[5] - i[10], i[21] - i[13], i[6] - i[18], i[8] - i[15], i[20] - i[16] };
141 :
142 0 : float c[] = { a[0] + a[1], a[0] + b[1], b[0] + a[1], b[0] + b[1], a[0] - a[1], a[0] - b[1], b[0] - a[1], b[0] - b[1] };
143 0 : float d[] = { a[2] + a[3], a[2] + b[3], b[2] + a[3], b[2] + b[3], a[2] - a[3], a[2] - b[3], b[2] - a[3], b[2] - b[3] };
144 0 : float e[] = { a[4] + a[5], a[4] + b[5], b[4] + a[5], b[4] + b[5], a[4] - a[5], a[4] - b[5], b[4] - a[5], b[4] - b[5] };
145 0 : float f[] = { a[6] + a[7], a[6] + b[7], b[6] + a[7], b[6] + b[7], a[6] - a[7], a[6] - b[7], b[6] - a[7], b[6] - b[7] };
146 0 : float g[] = { a[8] + a[9], a[8] + b[9], b[8] + a[9], b[8] + b[9], a[8] - a[9], a[8] - b[9], b[8] - a[9], b[8] - b[9] };
147 0 : float h[] = { a[10] + a[11], a[10] + b[11], b[10] + a[11], b[10] + b[11],
148 0 : a[10] - a[11], a[10] - b[11], b[10] - a[11], b[10] - b[11] };
149 :
150 0 : o[0] = fact * (+c[0] + d[0] + e[0] + f[0] + g[0] + h[0]);
151 0 : o[1] = fact * (-c[2] - d[5] + e[7] + f[7] + g[3] + h[0]);
152 0 : o[2] = fact * (-c[7] - d[1] - e[2] + f[4] + g[5] - h[6]);
153 0 : o[3] = fact * (-c[6] - d[0] + e[3] - f[7] - g[7] + h[2]);
154 0 : o[4] = fact * (-c[6] + d[6] - e[3] - f[4] - g[4] - h[3]);
155 0 : o[5] = fact * (-c[6] + d[4] - e[5] + f[7] + g[7] + h[3]);
156 0 : o[6] = fact * (-c[4] + d[1] - e[1] - f[2] + g[4] + h[5]);
157 0 : o[7] = fact * (-c[3] + d[0] - e[0] + f[3] - g[7] - h[7]);
158 0 : o[8] = fact * (-c[5] + d[0] + e[6] - f[3] - g[4] - h[4]);
159 0 : o[9] = fact * (-c[3] - d[6] + e[4] - f[5] + g[7] + h[7]);
160 0 : o[10] = fact * (-c[7] + d[2] + e[1] - f[1] - g[2] + h[4]);
161 0 : o[11] = fact * (-c[4] - d[3] + e[0] - f[0] + g[3] - h[7]);
162 0 : o[12] = fact * (-c[1] + d[3] + e[0] + f[6] - g[3] - h[4]);
163 0 : o[13] = fact * (-c[2] + d[5] - e[6] + f[4] - g[5] + h[7]);
164 0 : o[14] = fact * (-c[7] - d[7] + e[2] + f[1] - g[1] - h[2]);
165 0 : o[15] = fact * (-c[4] - d[4] - e[3] + f[0] - g[0] + h[3]);
166 0 : o[16] = fact * (-c[1] + d[7] + e[3] + f[0] + g[6] - h[3]);
167 0 : o[17] = fact * (-c[2] + d[4] + e[5] - f[6] + g[4] - h[5]);
168 0 : o[18] = fact * (-c[5] - d[7] - e[7] + f[2] + g[1] - h[1]);
169 0 : o[19] = fact * (-c[3] - d[4] - e[4] - f[3] + g[0] - h[0]);
170 0 : o[20] = fact * (-c[5] + d[7] + e[7] + f[3] + g[0] + h[6]);
171 0 : o[21] = fact * (-c[1] - d[2] + e[4] + f[5] - g[6] + h[4]);
172 0 : o[22] = fact * (-c[0] + d[3] - e[7] - f[7] + g[2] + h[1]);
173 0 : o[23] = fact * (-c[0] - d[3] - e[4] - f[4] - g[3] + h[0]);
174 0 : }
175 : } // namespace reverb_matrices
176 : } // namespace muse::audio::fx
177 :
178 : #endif // MUSE_AUDIO_REVERBMATRICES_H
|
