Artifcats when plotting SCQT results.

[jmerkt]

trafficstars

See here: https://github.com/jmerkt/rt-cqt/pull/4 Requires python bindings and plotting for debugging. Testing of different sample rates and block sizes.

[jumerckx]

To be clear, the screenshot I included in that pr was from my own non-functional implementation. Chances are your implementation does not exhibit this problem. I was mainly curious because of the comment in your code which says there is a bug somewhere.

[jmerkt]

It's worth examining anyways, doing a proper test in python, as I never did this for the sliding implementation. The comment in my code was, if I remember correctly, because of a different issue earlier and I forgot to delete it.

[jmerkt]

I've started testing here. So far I could not reproduce the original issue, but found some output attenuation issue when re-synthesizing:

[jumerckx]

Well that looks a lot better than my implementation :). I've messed around a bit with using all the samples of each block as you mention in your comment. The diagonal line on below spectrogram is from an exponential frequency sweep from 32.7032Hz to 15804.3Hz, the outer most frequencies according to getOctaveBinFreqs. It looks great except for the part after 80000 samples. I'm not sure what's happening there but chances are it's because of some mistake from my part.

When I plotted the spectrum of a sine wave with constant frequency of 15804.3Hz the result is as I'd expect, a line at the bottom of the image.

I didn't use Pybind as I dread having to set it up, instead I wrote some extremely gross code to write the CQT results to stdout and do some text wrangling in Python.

I combined the different sampling rates by just repeating samples for the undersampled cases, which isn't exactly by the rules but it gives a nice enough result I believe.

C++

for (unsigned i_block = 0; i_block < nBlocks; i_block++)
    {
        for (unsigned i = 0; i < blockSize; i++)
        {
            double t = (i+i_block*blockSize)*(1/samplerate);

            // exponential growth from 32.7Hz to 15.804KHz: https://www.desmos.com/calculator/eff93byluq
            double f = std::pow((1+9.4398654e-5), i+i_block*blockSize) + 31.7032;

            audioInputBlock.at(i) = sin(t*2*M_PI*f);// + (t>1)*sin(t*2*M_PI*32.7032) + (t>1.5)*sin(t*2*M_PI*15804.3);
        }

        cqt.inputBlock(audioInputBlock.data(), blockSize);

        for (unsigned i_octave = 0u; i_octave < OctaveNumber; i_octave++)
        {
            const size_t nSamplesOctave = cqt.getSamplesToProcess(i_octave);
            CircularBuffer<std::complex<double>> *octaveCqtBuffer = cqt.getOctaveCqtBuffer(i_octave);
            for (unsigned i_tone = 0u; i_tone < BinsPerOctave; i_tone++)
            {
                octaveCqtBuffer[i_tone].pullBlock(cqtDomainBuffers[i_octave][i_tone].data(), nSamplesOctave);
                for (size_t i_sample = 0u; i_sample < nSamplesOctave; i_sample++)
                {
                    std::complex<double> sample = cqtDomainBuffers[i_octave][i_tone][i_sample];
                    std::cout << sample.real() << "+" << sample.imag() << "j" << " ";
                }
                std::cout << std::endl;
                octaveCqtBuffer[i_tone].pushBlock(cqtDomainBuffers[i_octave][i_tone].data(), nSamplesOctave);
            }
            std::cout << std::endl;
        }
    }

Python

nOctaves = 9

octaves = []
with open("../output.txt") as f: # output from cpp was piped to this file
    lines = f.readlines()
    octave = []
    for l in lines:
        samples = []
        if l == "\n": # next octave
            octaves.append(np.concatenate(octave, axis=0))
            octave = []
        else:
            numbers = l.strip().split(" ")
            for n in numbers:
                samples.append(complex(n.strip().replace("+-", "-")))
            octave.append(np.array(samples).reshape((1, -1)))
new_octaves = []
for octave in range(nOctaves):
    new_octaves.append(np.concatenate(octaves[octave::nOctaves], axis=1))
octaves = new_octaves # a list of arrays, each of shape (tones, samples)

composite = octaves[0]
for i in range(1, nOctaves):
    # https://www.appsloveworld.com/numpy/100/13/how-to-interleave-numpy-ndarrays
    composite = np.concatenate((np.dstack([octaves[i] for _ in range(2**i)]).reshape(octaves[0].shape), composite), axis=0)

plt.imshow(np.abs(composite), aspect='auto', interpolation='bilinear')

[jmerkt]

Thanks for your help with further looking into it! Sine sweeps look good for me too, so I guess the forward transform is working pretty good. I will continue debugging the backward transform.