Convert showcqt.c to ShowCQT.cs

ShowCQT.cs

/*
 * Copyright (c) 2020 Muhammad Faiz <mfcc64@gmail.com>
 *
 * This library is free software; you can redistribute it and/or
 * modify it under the terms of the GNU Lesser General Public
 * License as published by the Free Software Foundation; either
 * version 3 of the License, or (at your option) any later version.
 *
 * This library is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 * Lesser General Public License for more details.
 *
 * You should have received a copy of the GNU Lesser General Public
 * License along with this library; if not, write to the Free Software
 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
 */

/* Audio visualization based on showcqt mpv/ffmpeg audio visualization */
/* See https://github.com/FFmpeg/FFmpeg/blob/master/libavfilter/avf_showcqt.c */

using System;

using UnityEngine;
using Unity.Mathematics;

using static Unity.Mathematics.math;

public static class _ShowCQT
{
    const int MAX_FFT_SIZE = 32768;
    const int MAX_WIDTH = 7680;
    const int MAX_HEIGHT = 4320;
    const int MAX_KERNEL_SIZE = 6 * 256 * 1024;
    const float MIN_VOL = 1.0f;
    const float MAX_VOL = 100.0f;

    public class ShowCQT
    {
        /* args */
        public double[][] input = new[] { new double[MAX_FFT_SIZE],new double[MAX_FFT_SIZE] };
        public Color[] output = new Color[MAX_WIDTH];

        /* tables */
        public double2[] exp_tbl = new double2[MAX_FFT_SIZE + MAX_FFT_SIZE / 4];
        public int[] perm_tbl = new int[MAX_FFT_SIZE / 4];
        public double[] attack_tbl = new double[MAX_FFT_SIZE / 8];
        public byte[] padding = new byte[1024];

        /* buffers */
        public double2[] fft_buf = new double2[MAX_FFT_SIZE+128];
        public Color[] color_buf = new Color[MAX_WIDTH * 2];
        public float[] rcp_h_buf = new float[MAX_WIDTH];

        /* kernel */
        public (int len, int start)[] kernel_index = new (int len, int start)[MAX_WIDTH * 2];
        public double[] kernel = new double[MAX_KERNEL_SIZE];

        /* props */
        public int width;
        public int height;
        public int aligned_width;
        public int fft_size;
        public int t_size;
        public int attack_size;
        public float sono_v;
        public float bar_v;
        public bool prerender;
    }

    static ShowCQT cqt = new ShowCQT();
    public static int FFTSize => cqt.fft_size;

    public static double[] get_input_array(int index)
    {
        return cqt.input[index];
    }

    public static Color[] get_output_array()
    {
        return cqt.output;
    }

    public static Color[] get_color_array()
    {
        return cqt.color_buf;
    }

    static int revbin(int x, int bits)
    {
        int m = 0x55555555;
        x = ((x & m) << 1) | ((x & ~m) >> 1);
        m = 0x33333333;
        x = ((x & m) << 2) | ((x & ~m) >> 2);
        m = 0x0F0F0F0F;
        x = ((x & m) << 4) | ((x & ~m) >> 4);
        m = 0x00FF00FF;
        x = ((x & m) << 8) | ((x & ~m) >> 8);
        m = 0x0000FFFF;
        x = ((x & m) << 16) | ((x & ~m) >> 16);
        return (x >> (32 - bits)) & ((1 << bits) - 1);
    }

    static void gen_perm_tbl(int bits)
    {
        int n = 1 << bits;
        for (int x = 0; x < n; x++)
            cqt.perm_tbl[x] = revbin(x, bits);
    }

    static double2 ComplexMultiply(in double2 a,in double2 b) => new double2((a.x * b.x) - (a.y * b.y),(a.x * b.y) + (a.y * b.x));

    static void gen_exp_tbl(int n)
    {
        double mul;
        for (int k = 2; k < n; k *= 2)
        {
            mul = 2.0 * PI_DBL / k;
            for (int x = 0; x < k / 2; x++)
                cqt.exp_tbl[k+x] = new double2(cos(mul * x), -sin(mul * x));
            mul = 3.0 * PI_DBL / k;
            for (int x = 0; x < k / 2; x++)
                cqt.exp_tbl[k + k / 2 + x] = new double2( cos(mul * x), -sin(mul * x) );
        }
        mul = 2.0 * PI_DBL / n;
        for (int x = 0; x < n / 4; x++)
            cqt.exp_tbl[n+x] = new double2( cos(mul * x), -sin(mul * x) );
    }


    static void fft_butterfly(Span<double2> v, int n, int q)
    {
        double2 v0, v1, v2, v3;
        double2 a02, a13, s02, s13;

        v0 = v[0];
        v2 = v[q]; /* bit reversed */
        v1 = v[2 * q];
        v3 = v[3 * q];
        a02 = v0 + v2;
        s02 = v0 - v2;
        a13 = v1 + v3;
        s13 = v1 - v3;
        v[0] = a02 + a13;
        v[q] = new double2(s02.x + s13.y,s02.y - s13.x);
        v[2 * q] = a02 - a13;
        v[3 * q] = new double2(s02.x - s13.y,s02.y + s13.x);

        for (int x = 1; x < q; x++)
        {
            v0 = v[x];
            v2 = ComplexMultiply(cqt.exp_tbl[2*q + x], v[q+x]); /* bit reversed */
            v1 = ComplexMultiply(cqt.exp_tbl[4*q + x], v[2*q+x]);
            v3 = ComplexMultiply(cqt.exp_tbl[3*q + x], v[3*q+x]);
            a02 = v0 + v2;
            s02 = v0 - v2;
            a13 = v1 + v3;
            s13 = v1 - v3;
            v[x] = a02 + a13;
            v[q+x] = new double2(s02.x + s13.y,s02.y - s13.x);
            v[2*q+x] = a02 - a13;
            v[3*q+x] = new double2(s02.x - s13.y,s02.y + s13.x);
        }
    }

    static void fft_calc_(int n,Span<double2> v)
    {
        int q = n / 4;
        if(q < 2)
            return;

        if(q == 2)
        {
            var v0 = v[0];
            var v1 = v[1];
            v[0] = v0 + v1;
            v[1] = v0 - v1;
            return;
        }

        fft_calc_(q,v);
        fft_calc_(q,v.Slice(q));
        fft_calc_(q,v.Slice(2*q));
        fft_calc_(q,v.Slice(3*q));
        fft_butterfly(v, n, q);                     
    }

    static void fft_calc(double2[] v, int n)
    {
        switch (n)
        {
            case 1024:
            case 2048:
            case 4096:
            case 8192:
            case 16384:
            case 32768:
                fft_calc_(n,v);
                break;
        }
    }


    public static int init(int rate, int width, int height, float bar_v, float sono_v, int super)
    {
        if (height <= 0 || height > MAX_HEIGHT || width <= 0 || width > MAX_WIDTH)
            return 0;

        cqt.width = width;
        cqt.height = height;
        cqt.aligned_width = width;

        cqt.bar_v = (bar_v > MAX_VOL) ? MAX_VOL : (bar_v > MIN_VOL) ? bar_v : MIN_VOL;
        cqt.sono_v = (sono_v > MAX_VOL) ? MAX_VOL : (sono_v > MIN_VOL) ? sono_v : MIN_VOL;

        if (rate < 8000 || rate > 100000)
            return 0;

        int bits = (int)ceil(log(rate / 3.0) / LN2_DBL);
        if (bits > 20 || bits < 10)
            return 0;
        cqt.fft_size = 1 << bits;
        if (cqt.fft_size > MAX_FFT_SIZE)
            return 0;

        gen_perm_tbl(bits - 2);
        gen_exp_tbl(cqt.fft_size);

        cqt.attack_size = (int)ceil(rate / 30.0);
        for (int x = 0; x < cqt.attack_size; x++)
        {
            double y = PI_DBL * x * 30 / (double)rate;
            cqt.attack_tbl[x] = 0.355768 + 0.487396 * cos(y) + 0.144232 * cos(2 * y) + 0.012604 * cos(3 * y);
        }

        cqt.t_size = cqt.width * (1 + super);
        double log_base = log(20.01523126408007475);
        double log_end = log(20495.59681441799654);
        for (int f = 0, idx = 0; f < cqt.t_size; f++)
        {
            double freq = exp(log_base + (f + 0.5) * (log_end - log_base) * (1.0 / cqt.t_size));

            if (freq >= 0.5 * rate)
            {
                cqt.kernel_index[f].len = 0;
                cqt.kernel_index[f].start = 0;
                continue;
            }

            double p = 3;
            double ovt = 5;
            double secondHarmonic = 384;

//			double tlen = 384 * 0.33 / (384 / 0.17 + 0.33 * freq / (1 - 0.17)) + 384 * 0.33 / (0.33 * freq / 0.17 + 384 / (1 - 0.17));
//			double tlen = 0.33 / (1/0.17 + 0.33 * freq / (384 * (1 - 0.17))) + 0.33 / (0.33 * freq / (384 * 0.17) + 1 / (1 - 0.17));
//			double tlen = 1/(3/0.17 + freq/(384 * (1 - 0.17))) + 1/(freq/(384 * 0.17) + 3/(1 - 0.17));
//			double tlen = 1/(3/0.17 + freq/384/(1 - 0.17)) + 1/(freq/384/0.17 + 3/(1 - 0.17));
//			double tlen = 1/(3/(1/2 - 1/3) + freq/384/(1/2 + 1/3)) + 1/(freq/384/(1/2 - 1/3) + 3/(1/2 + 1/3));
//			double tlen = 5 * (1/(3 + 5*freq/384) + 1/(3*5 + freq/384)) / 6;
			double tlen = ovt * (1/(p + (ovt * freq / secondHarmonic)) + 1/((p * ovt) + (freq / secondHarmonic))) / 6;
            double flen = 8.0 * cqt.fft_size / (tlen * rate);
            double center = freq * cqt.fft_size / rate;
            int start = (int)ceil(center - 0.5 * flen);
            int end = (int)floor(center + 0.5 * flen);
            int len = end - start + 1;

            if (idx + len + 1000 > MAX_KERNEL_SIZE)
                return 0;
            cqt.kernel_index[f].len = len;
            cqt.kernel_index[f].start = start;

            for (int x = start; x < start + len; x++)
            {
                if (x > end)
                {
                    cqt.kernel[idx + x - start] = 0;
                    continue;
                }
                double sign = (x & 1) > 0 ? -1 : 1;
                double y = 2.0 * PI_DBL * (x - center) * (1.0 / flen);
                double w = 0.355768 + 0.487396 * cos(y) + 0.144232 * cos(2 * y) + 0.012604 * cos(3 * y);
                w *= sign / cqt.fft_size;
                cqt.kernel[idx + x - start] = w;
            }

            idx += len;
        }
        return cqt.fft_size;
    }

    static double2 cqt_calc(Span<double> kernel, int start, int len)
    {
        var a = double2(0);
        var b = double2(0);

        for (int m = 0; m < len; m++)
        {
            double u = kernel[m];
            a += u * cqt.fft_buf[start + m];
            b += u * cqt.fft_buf[cqt.fft_size - start - m];
        }

        double2 v0 = new double2(a.x + b.x, a.y - b.y);
        double2 v1 = new double2(b.y + a.y, b.x - a.x);
        return new double2( lengthsq(v0), lengthsq(v1) );
    }

    public static void calc()
    {
        int fft_size_h = cqt.fft_size >> 1;
        int fft_size_q = cqt.fft_size >> 2;
        int shift = fft_size_h - cqt.attack_size;

        for (int x = 0; x < cqt.attack_size; x++)
        {
            int i = 4 * cqt.perm_tbl[x];
            cqt.fft_buf[i] = new double2(cqt.input[0][shift+x], cqt.input[1][shift+x]);
            cqt.fft_buf[i+1] = cqt.attack_tbl[x] * new double2(cqt.input[0][fft_size_h+shift+x], cqt.input[1][fft_size_h+shift+x]);
            cqt.fft_buf[i+2] = new double2(cqt.input[0][fft_size_q+shift+x], cqt.input[1][fft_size_q+shift+x]);
            cqt.fft_buf[i+3] = 0;
        }

        for (int x = cqt.attack_size; x < fft_size_q; x++)
        {
            int i = 4 * cqt.perm_tbl[x];
            cqt.fft_buf[i] = new double2(cqt.input[0][shift+x], cqt.input[1][shift+x]);
            cqt.fft_buf[i+1] = 0;
            cqt.fft_buf[i+2] = new double2(cqt.input[0][fft_size_q+shift+x], cqt.input[1][fft_size_q+shift+x]);
            cqt.fft_buf[i+3] = 0;
        }

        fft_calc(cqt.fft_buf, cqt.fft_size);

        for (int x = 0, m = 0; x < cqt.t_size; x++)
        {
            int len = cqt.kernel_index[x].len;
            int start = cqt.kernel_index[x].start;
            if (len <= 0)
            {
                cqt.color_buf[x] = Color.black;
                continue;
            }

            var r = cqt_calc(cqt.kernel.AsSpan(m), start, len);

            double a = 0.5 * (r.x + r.y);

            cqt.color_buf[x] = (Vector4)(float4)double4(sqrt(cqt.sono_v * sqrt(double3(r.x,a,r.y))),cqt.bar_v * sqrt(a));

            m += len;
        }

        if (cqt.t_size != cqt.width)
        {
            for (int x = 0; x < cqt.width; x++)
            {
                cqt.color_buf[x] = 0.5f * (cqt.color_buf[2*x] + cqt.color_buf[2*x+1]);
            }
        }

        cqt.prerender = true;
    }

    static void prerender()
    {
        for (int x = 0; x < cqt.width; x++)
        {
            cqt.color_buf[x] = (Vector4)clamp((Vector4)cqt.color_buf[x],0,1);
        }

        for (int x = 0; x < cqt.aligned_width; x++)
            cqt.rcp_h_buf[x] = 1.0f / (cqt.color_buf[x].a + 0.0001f);

        cqt.prerender = false;
    }

    public static void render_line_alpha(int y, float a)
    {
        if (cqt.prerender)
            prerender();


        if (y >= 0 && y < cqt.height)
        {
            float ht = (cqt.height - y) / (float)cqt.height;
            for (int x = 0; x < cqt.width; x++)
            {
                if (cqt.color_buf[x].a <= ht)
                {
                    cqt.output[x] = new Color(a,a,a,a);
                }
                else
                {
                    var c = (cqt.color_buf[x].a - ht) * cqt.rcp_h_buf[x] * new float4((Vector4)cqt.color_buf[x]).xyz;
                    cqt.output[x] = new Color(c.x,c.y,c.z,a) * 100;
                }
            }
        }
        else
        {
            for (int x = 0; x < cqt.width; x++)
            {
                var c = new float4((Vector4)cqt.color_buf[x]).xyz;
                cqt.output[x] = new Color(c.x,c.y,c.z,a) * 100;
            }
        }
    }

    public static void render_line_opaque(int y)
    {
        render_line_alpha(y,1);
    }

    public static void set_volume(float bar_v, float sono_v)
    {
        cqt.bar_v = (bar_v > MAX_VOL) ? MAX_VOL : (bar_v > MIN_VOL) ? bar_v : MIN_VOL;
        cqt.sono_v = (sono_v > MAX_VOL) ? MAX_VOL : (sono_v > MIN_VOL) ? sono_v : MIN_VOL;
    }

    public static void set_height(int height)
    {
        cqt.height = (height > MAX_HEIGHT) ? MAX_HEIGHT : (height > 1) ? height : 1;
    }

    public static bool detect_silence(float threshold)
    {
        for (int x = 0; x < cqt.fft_size; x++)
            if (dot(cqt.input[0][x],cqt.input[1][x]) > threshold)
                return false;
        return true;
    }
}