diff --git a/varsamp.c b/varsamp.c index aeef36e..db0dec0 100644 --- a/varsamp.c +++ b/varsamp.c @@ -67,11 +67,15 @@ void calc_varsamp (VARSAMP a) line. Transposing makes each phase contiguous; hshift() interpolates between phases hidx and hidx+1, hence R+1 of them. */ int p, m; - double* h = fir_bandpass(a->ncoef, fc_norm_low, fc_norm_high, (double)a->R, 1, 0, (double)a->R * a->gain); - a->hp = (double *)malloc0 ((size_t)(a->R + 1) * a->rsize * sizeof (double)); - for (p = 0; p <= a->R; p++) - for (m = 0; m < a->rsize; m++) - a->hp[(size_t)p * a->rsize + m] = h[p + (size_t)m * a->R]; + const int R = a->R, rsize = a->rsize; + double* h = fir_bandpass(a->ncoef, fc_norm_low, fc_norm_high, (double)R, 1, 0, (double)R * a->gain); + // every element is written below, so skip malloc0()'s memset of ~1 MB + a->hp = (double *)_aligned_malloc ((size_t)(R + 1) * rsize * sizeof (double), 16); + // walk h forward (p is its fast axis) so the prefetcher sees a linear + // stream; h is cold here, straight from fir_bandpass() + for (m = 0; m < rsize; m++) + for (p = 0; p <= R; p++) + a->hp[(size_t)p * rsize + m] = h[p + (size_t)m * R]; _aligned_free (h); } a->ringI = (double *)malloc0(a->rsize * sizeof(double));