1
$\begingroup$

I'm reading up on fourier theory, especially the transforms. I implement the math as spectrograms in C++ to get a better understanding of what is going on.

I've made an implementation of the short time discrete fourier transform which works well, however as for transforms derived from the fourier transform, they seem to be bounded by the so called uncertainty principle which basically is a tradeoff between time and frequency resolution - which in effect means, that if you want better frequency resolution, you need a larger sample size - which again means you lose time resolution.

Now, I always wanted a spectrum analyzer with increased resolution at lower frequencies, so that a logarithmic display of the spectrum would have equal resolution over the entire frequency range - this is obviously not the case with FFTs or short time DFTs.

Knowing this, I went ahead and researched constant Q transforms and wavelets (from which the constant Q should be a special case of) which promises geometric spacing and constant q of the bins compared to your standard fast fourier transform.

I tried to implement the algorithm on this page, but i find the definitions poorly explained and similarly confusing, so it doesn't currently work. Here's my current implementation - i would appreciate if anyone could find the error(s):

typedef std::uint_fast32_t fint_t;

template<typename Scalar, class Vector>
    std::complex<Scalar> constantQTransform(Vector input, 
        size_t size,
        fint_t kbin, 
        Scalar lowestFreq, 
        Scalar numFiltersPerOctave, 
        Scalar sampleRate)
    {
        // "spectral width" per filter
        Scalar filterWidth = std::pow(2, Scalar(1.0) / numFiltersPerOctave); // r = 2 ^ (1 / n)
        filterWidth = std::pow(filterWidth, kbin) * lowestFreq; // r ^ k * fmin == fk

        // centre frequency
        Scalar centreFrequency = std::pow(2, kbin * 1 / numFiltersPerOctave) * lowestFreq;

        // window length for bin
        Scalar windowLength = sampleRate / filterWidth; // fs / fk == N[k]
        fint_t end = static_cast<fint_t>(std::floor(windowLength));
        const Scalar Q = centreFrequency / filterWidth;
        // bounds check
        if (end > size)
            return 0;
        auto hammingWindow = [&](fint_t n) // W[k, n]
        {
            auto const a = 25.0 / 46.0;
            return a - (1 - a) * std::cos((TAU * n) / windowLength);
        };

        std::complex<Scalar> acc;
        Scalar real, imag, sample;

        for (fint_t n = 0; n < end; ++n)
        {
            sample = input[n] * hammingWindow(n);
            real = std::cos((TAU * Q * n) / windowLength) * sample;
            imag = std::sin((TAU * Q * n) / windowLength) * sample;
            acc += std::complex<Scalar>(real, -imag);

        }
        return acc / windowLength; // normalize
    }

Here's the spectrum of this function: Spectrum of constantQ for a square wave

And this is how it should look like, using my short time dft with logarithmic spacing: Spectrum of short time dft for a square wave

So obviously, something is wrong - I'm not sure what, but it 'reacts' to different frequencies, ie. the first two bubbles move accordingly.

For reference, here's my short time (unoptimized) dft:

template<typename Scalar, class Vector, bool scale = false>
    std::complex<Scalar> fourierTransform(const Vector & data, 
        Scalar frequency, 
        Scalar sampleRate, 
        size_t a, 
        size_t b)
    {
        std::complex<Scalar> acc = 0;
        Scalar real = 0, imag = 0;
        //std::complex<Scalar> z = 0; 

        Scalar w = frequency * 2 * PI / sampleRate;
        double i = 1.0;
        for (int t = a; t < b; ++t)
        {
            // e^(-i*w*t) == cos(w*t) + -i * sin(w * t) ?
            //real = std::cos(w * t);
            //imag = -i * std::sin(w * t);
            //z = data[t] * std::complex<Scalar>(std::cos(w * t), -std::sin(w * t));
            // integrate(a, b, { f(t) * e^(-i*w*t) }) <- fourier transform
            acc += Scalar(data[t]) * std::complex<Scalar>(std::cos(w * t), -std::sin(w * t));
        }
        if (scale)
        {
            return acc / std::complex<Scalar>((b - a) / 2.0);
        }
        else
            return acc;

    }

However, if we take a look at the constant Q transform:

Constant Q transform

Where this is the standard short time fourier transform:

Short time fourier transform

What basically seems to happen is that it runs a simple short time dft on a windowed subset of the signal, where the subsets size is a function of the frequency being analyzed, thereby artificially generating a small bandwidth for lower frequencies and vice-versa.

This leads me to a couple of conclusions and questions, i hope you can help me with...

  • Given that the constant q effectively just is a dft, it is limited by the very same uncertainty principle that the fourier transform suffers from?
  • is this the case with wavelet transforms like morlet as well?
  • and in general, is it just completely impossible to achieve higher resolution without directly altering the sample size of the transform?
  • Since the constant Q transform only takes a subset of the samples for higher frequencies, am i correct in assuming that for the following sample, it may very well report no high frequency content?:

sample

  • If the previous 4 questions are true, why would anyone ever want to compute the constant q transform given that it can be calculated directly, and more precise, using the short time fourier transform? The only answer i can come up with, is that you can compute it at fft speed, where your bins are geometrically spaced - instead of linearly - but still with no additional resolution.
  • If i am correct in the previous assumptions, then don't bother helping with the constant-Q code, because i might as well use my DFT then!

Thanks for your interest in an newbie's confused mind - i hope you can help with some of my questions.

$\endgroup$
  • $\begingroup$ The merhod of spectrograms is low-resolution spectral estimation. If you want to achieve much better resolution without dramatically increasing the number of observed samples, you should try more sophisticated methods. Read about Capon or MUSIC spectral estimations techniques. user.it.uu.se/~ps/SAS-new.pdf will be a good start. $\endgroup$ – Serj Aug 18 '14 at 3:22
  • $\begingroup$ What is the Big O complexity of Short Time Fourier Transform and Constant Q Transform ? $\endgroup$ – user16198 Jun 11 '15 at 11:45
  • $\begingroup$ This should be a separate question, or maybe should be asked on stackoverflow. $\endgroup$ – JRE Jun 11 '15 at 12:13
  • $\begingroup$ This does not provide an answer to the question. To critique or request clarification from an author, leave a comment below their post - you can always comment on your own posts, and once you have sufficient reputation you will be able to comment on any post. $\endgroup$ – jojek Jun 11 '15 at 13:10
0
$\begingroup$

Given that the constant q effectively just is a dft, it is limited by the very same uncertainty principle that the fourier transform suffers from?

Yes

is this the case with wavelet transforms like morlet as well?

Yes

and in general, is it just completely impossible to achieve higher resolution without directly altering the sample size of the transform?

Yes

Since the constant Q transform only takes a subset of the samples for higher frequencies, am i correct in assuming that for the following sample, it may very well report no high frequency content?:

With sufficient overlap between consecutive transforms, you can avoid that.

If the previous 4 questions are true, why would anyone ever want to compute the constant q transform given that it can be calculated directly, and more precise, using the short time fourier transform? The only answer i can come up with, is that you can compute it at fft speed, where your bins are geometrically spaced - instead of linearly - but still with no additional resolution.

On low frequency, Constant Q transform has higher frequency-domain resolution but lower time-domain resolution.

On high frequency, Constant Q transform has lower frequency-domain resolution but higher time-domain resolution.

It has more similarity with human hearing compared to plain STFT.

Maybe you want to take a look at my (not so) constant Q transform implementation. Take a look at cqt_calc() and init_cqt(). Here is a screenshot:

showcqt

taken from my mpv-scripts.

$\endgroup$
  • $\begingroup$ Thank you for the clarifications. "On low frequency, Constant Q transform has higher frequency-domain resolution" should still be no greater than an equivalent interpolated dft, right? For my particular problem at that point I instead went with a complex resonator bank, where I could tune spacing and resonator bandwidth to obtain logarithmic resolution. $\endgroup$ – Shaggi Jul 20 at 13:32
  • $\begingroup$ @Shaggi But constant Q transform has higher time-domain resolution at high frequency compared to the equivalent interpolated dft. $\endgroup$ – mfcc64 Jul 20 at 16:24

Not the answer you're looking for? Browse other questions tagged or ask your own question.