Planck-Taper Window

[xkzl]

Hello,

I have implemented a window, maybe some people might be interested in. https://en.wikipedia.org/wiki/Window_function Maybe it can be implemented in the official repo ? ( I am still using KFR 4.0, but maybe this can still work in 5.0)

template <typename T>
struct expression_planck_taper : input_expression
{
    using value_type = T;

    expression_planck_taper(size_t size, T epsilon = 0.1, window_symmetry symmetry = window_symmetry::symmetric)
        : linspace(size, symmetry), epsilon(epsilon), m_size(size)
    {
    }
    template <size_t N>
    KFR_INTRINSIC friend vec<T, N> get_elements(const expression_planck_taper& self, cinput_t cinput, size_t index,
                                                vec_shape<T, N> y)
    {
        if(index == 0 || index == self.m_size) return 0;

        if(index >= self.m_size - self.epsilon*self.m_size)
            index = self.m_size - index - 1;

        if(index < self.epsilon*self.m_size) {

            return T(1.) / (T(1.) + exp(
                self.epsilon * self.m_size / index - 
                self.epsilon * self.m_size / ( T(self.epsilon * self.m_size) - index) 
            ));
        }

        using TI           = utype<T>;
        const vec<TI, N> i = enumerate(vec_shape<TI, N>()) + static_cast<TI>(index);
        return select(i < static_cast<TI>(self.m_size), T(1), T(0));
    }
    size_t size() const { return m_size; }

private:
    window_linspace_0_1<T> linspace;
    T epsilon;
    size_t m_size;
};

[dancazarin]

Hello,

In just released KFR 5.0.3 there is now an implementation of Planck-taper window. https://github.com/kfrlib/kfr/blob/bae07d0b2936e0f84da5bf7f7a02a140d8c1346e/include/kfr/dsp/window.hpp#L395-L397

@xKZL , as of the code you suggest. There are some significant problems with the code.

Example of its output first (on AVX2 with 8 numbers computed in parallel):

univector<float> f = expression_planck_taper<float>(64, 0.25);
println(f);
// output:
        0         0         0         0         0         0         0         0
      0.5       0.5       0.5       0.5       0.5       0.5       0.5       0.5
        1         1         1         1         1         1         1         1
        1         1         1         1         1         1         1         1
        1         1         1         1         1         1         1         1
        1         1         1         1         1         1         1         1
        1         1         1         1         1         1         1         1
 0.375677  0.375677  0.375677  0.375677  0.375677  0.375677  0.375677  0.375677

Obviously it's not a smooth window function. Below is the reasons. 1) KFR is designed to run computations in parallel using SIMD. It's good for performance but makes writing of correct code harder. In the example above get_elements was called 8 times every time producing vec of 8 numbers, not 64 times producing 1 number as if the code was scalar. 2) Under the right circumstances, index parameter should never be used (it's equal to the index of the first element in the vec to be generated). linspace is the right thing to get a number in the desired range linspace isn't used anywhere in the code you suggest. 3) index == self.m_size and i < static_cast<TI>(self.m_size) checks are not needed. If you use the expression correctly, it will never be called with index >= size. If you need zero padding after expression, you may use padded function. 4) The code doesn't work if size is not even. 5) Window function symmetry parameter is ignored. Instead m_size is used directly which is incorrect.

After all the above errors have been fixed (which is mainly related to parallel computations) Your code prints the following values for (size=8, epsilon=0.25)

        0       0.5         1         1         1         1       0.5         0

But this is incorrect Planck-taper window coefficients.

The reason is incorrect use of m_size in context of window function computation. m_size is not equal to $T$ from the original formula from https://arxiv.org/pdf/1003.2939.pdf (nor $N$ from wikipedia). It's the number of points of the window function we want to get. But, if you want to get symmetric output (zeros in both the first and the last points, in [0] and [m_size-1]), then $T$ is equal to the distance between the first and the last points, so it is equal to m_size-1.

For $T=8$ and $\epsilon=0.25$ the original formula gives the following numbers:

{0, 0.641834, 1, 1, 1, 1, 0.641833, 0}

For $T=9$ and $\epsilon=0.25$ it gives:

{0, 0.5, 1, 1, 1, 1, 1, 0.5, 0}

Here is the demonstration of the problem: https://www.desmos.com/calculator/dd6drms9mf The graph shows nine points with values $[0, 0.25, 1, 1, 1, 1, 1, 0.25, 0]$ while your algorithm never produces the correct values (nor with size=8 nor with size=9).

In KFR 5.0.3 Planck-taper window function is implemented according to the article "A tapering window for time-domain templates and simulated signals in the detection of gravitational waves from coalescing compact binaries", arxiv:1003.2939. Both symmetric and periodic variants are implemented.

[dancazarin]

Closed as completed