1313#include < Eigen/Sparse>
1414#include < Eigen/IterativeLinearSolvers>
1515
16- // using FloatSolver = Eigen::SparseLU<Eigen::SparseMatrix<float>>;
17- using FloatSolver = Eigen::BiCGSTAB<Eigen::SparseMatrix<float >>;
16+ using FloatSolver = Eigen::SparseLU<Eigen::SparseMatrix<float >>;
17+ // using FloatSolver = Eigen::BiCGSTAB<Eigen::SparseMatrix<float>>;
1818
1919
2020namespace audio_plugin {
@@ -34,7 +34,7 @@ void compute_f_optimized(juce::AudioBuffer<float>& f, // size(f) == size(
3434 const juce::AudioBuffer<float >& g, // size(f) == size(x) == size(g)
3535 const std::vector<float >& two_pow_n, // size(two_pow_n) == max_terms
3636 const std::vector<float >& weights, // size(weights) == max_terms
37- int max_terms = 20 /* Must be divisible by 4! */ ) {
37+ int max_terms = 20 ) {
3838 const int N = g.getNumSamples ();
3939 const int numChannels = f.getNumChannels ();
4040 const float inv_g_step = static_cast <float >(N);
@@ -48,62 +48,17 @@ void compute_f_optimized(juce::AudioBuffer<float>& f, // size(f) == size(
4848 const float xi = x_data[i];
4949 float sum = 0 .0f ;
5050
51- for (int n = 0 ; n < max_terms; n += 4 ) {
52- const float arg0 = xi * two_pow_n[n];
53- const float arg1 = xi * two_pow_n[n+1 ];
54- const float arg2 = xi * two_pow_n[n+2 ];
55- const float arg3 = xi * two_pow_n[n+3 ];
56-
57- const int idx0 = static_cast <int >((arg0 - std::floor (arg0)) * inv_g_step + 0 .5f );
58- const int idx1 = static_cast <int >((arg1 - std::floor (arg1)) * inv_g_step + 0 .5f );
59- const int idx2 = static_cast <int >((arg2 - std::floor (arg2)) * inv_g_step + 0 .5f );
60- const int idx3 = static_cast <int >((arg3 - std::floor (arg3)) * inv_g_step + 0 .5f );
61-
62- sum += weights[n] * g_data[idx0];
63- sum += weights[n+1 ] * g_data[idx1];
64- sum += weights[n+2 ] * g_data[idx2];
65- sum += weights[n+3 ] * g_data[idx3];
66- }
67-
68- f_data[i] = sum;
69- }
70- }
71- }
72-
73-
74- /* void compute_f(juce::AudioBuffer<float>& f,
75- const std::vector<float>& x,
76- const juce::AudioBuffer<float>& g,
77- const std::vector<float>& two_pow_n,
78- const std::vector<float>& weights,
79- int max_terms = 20) {
80- const int N = g.getNumSamples();
81- const int numChannels = f.getNumChannels();
82-
83- const float inv_g_step = static_cast<float>(N);
84- const float* x_data = x.data();
85-
86- for (int ch = 0; ch < numChannels; ++ch) {
87- const float* g_data = g.getReadPointer(ch);
88- float* f_data = f.getWritePointer(ch);
89-
90- for (int i = 0; i < N; ++i) {
91- float sum = 0.0f;
92- const float xi = x_data[i];
93-
9451 for (int n = 0 ; n < max_terms; ++n) {
9552 const float arg = xi * two_pow_n[n];
96- const float fractional = arg - std::floor(arg);
97- const int idx = static_cast<int>(fractional * inv_g_step + 0.5f);
98- //const int safe_idx = idx < N ? idx : N - 1;
99-
100- sum += weights[n] * g_data[idx];
53+ const int idx = static_cast <int >((arg - std::floor (arg)) * inv_g_step + 0 .5f );
54+ const int safe_idx = idx < N ? idx : N - 1 ;
55+ sum += weights[n] * g_data[safe_idx];
10156 }
10257
10358 f_data[i] = sum;
10459 }
10560 }
106- }*/
61+ }
10762
10863
10964void fractalize (const std::vector<float > &x_grid,
0 commit comments