drmikeando (u/drmikeando)

drmikeando commented on Gaussian integration is cool rohangautam.github.io/blo... · Posted by u/beansbeansbeans

drmikeando · 7 months ago

My mind was exploded by this somewhat similar technique https://en.wikipedia.org/wiki/Tanh-sinh_quadrature - it uses a similar transformation of domain, but uses some properties of optimal quadrature for infinite sequences in the complex plane to produce quickly convergent integrals for many simple and pathological cases.

drmikeando commented on Show HN: Resonate – real-time high temporal resolution spectral analysis alexandrefrancois.org/Res... · Posted by u/arjf

drmikeando · 8 months ago

You can view this result as the convolution of the signal with an exponentially decaying sine and cosine.

That is, `y(t') = integral e^kt x(t' - t) dt`, with k complex and negative real part.

If you discretize that using simple integration and t' = i dt, t = j dt you get

    y_i = dt sum_j e^(k j dt) x_{i - j}
    y_{i+1} = dt sum_j e^(k j dt) x_{i+1 - j}
            = (dt e^(k dt) sum_j' e^(k j' dt) x_{i - j'}) + x_i 
            = dt e^(k dt) y_i + x_i

If we then scale this by some value, such that A y_i = z_i we can write this as

    z_{i+1} = dt e^(k dt) z_i + A x_i

Here the `dt e^(k dt)` plays a similar role to (1-alpha) and A is similar to P alpha - the difference being that P changes over time, while A is constant.

We can write `z_i = e^{w dt i} r_i` where w is the imaginary part of k

   e^{w dt (i+1)} r_{i+1} = dt e^(k dt) e^{w dt i} r_i + A x_i
             r_{i+1} = dt e^((k - w) dt) r_i + e^{-w dt (i+1) } A x_i
                     = (1-alpha) r_i + p_i x_i

Where p_i = e^{-w dt (i+1) } A = e^{-w dt ) p_{i-1} Which is exactly the result from the resonate web-page.

The neat thing about recognising this as a convolution integral, is that we can use shaping other than exponential decay - we can implement a box filter using only two states, or a triangular filter (this is a bit trickier and takes more states). While they're tricky to derive, they tend to run really quickly.

drmikeando commented on Going faster by duplicating code voidstar.tech/code_duplic... · Posted by u/voidstarcpp

drmikeando · 2 years ago

IMO the reason the compiler doesn't add special cases for the simplest version is that it doesn't know which of its _many_ special cases to use. If you actually use the unoptimised version of the code like

    void withSwitch(vector<int>& Values, bool v) {
      if (v) {
        multiply1(Values, 2.0);
      } else {
        multiply1(Values, 3.0);
      }
    }

Then it actually inlines the code and optimises each one correctly, as it has context about which special cases are available. (Doesn't even need the `inline` keyword for this at `-O2`)

You can see the code here: https://godbolt.org/z/5beeYe77a