Re: Convolution, channel-selectiveness of effects and spectrograms

Sergei's analysis seems not to have made it to the sox-users archive,
so herewith:

On Sat, 29 Nov 2025 at 07:23, Sergei Steshenko <[email protected]> wrote:
> Let's start from very basic things.
>
> FT (Fourier Transform - not necessarily DFT - Discrete Fourier Transform) is a subset of transforms in which the function to be transformed is decomposed, i.e. presented as a sum of orthogonal functions - see https://en.wikipedia.org/wiki/Orthogonal_functions .
>
> And on DFT specifically, see, for example, "Nailing Fourier Series and Orthogonal Decompositions" - https://medium.com/@ivethium/nailing-fourier-series-and-orthogonal-decompositions-ff94b9d7476a . And/or https://fftw.org/fftw3_doc/What-FFTW-Really-Computes.html#What-FFTW-Really-Computes .
>
> So, in DFT the 'cos' and 'sin' components are the set of orthogonal functions to do the decomposition on. I emphasize the point - each cos(2 * pi * Bn / N) and sin(2 * pi * Bn / N) (Bn is bin number and N is number of points DFT) is orthogonal each other when Bn is different - this is by definition of orthogonality.
>
> Mathematically the above cos(2 * pi * Bn / N) and sin(2 * pi * Bn / N) are orthogonal, but computationally they are not necessarily orthogonal because of insufficient number of bits.
>
> If we take DFT bin number 0 (Bn is 0), cos(2 * pi * Bn / N) becomes cos(0) which is 1, and if bin number 1 (Bn is 1) cos(2 * pi * Bn / N) becomes cos(2 * pi  / N) . With large enough N because of limited number of bits the 2 * pi  / N will become effectively 0 and thus cos(0) will become equal to cos(2 * pi  / N), thus instead of two orthogonal (i.e. different) functions we'll have just ONE (cos(0)) function.
>
> So, we'll break the orthogonality assumption and thus we will NOT be performing DFT.
>
> A quick check in Julia:
>
> "
>
> julia> let; x::Float32 = 1.0 / 4096; cos(x) end
> 1.0f0
>
> julia> let; x::Float32 = 1.0 / 2048; cos(x) end
> 0.9999999f0
> ",
>
> i.e. N = 4096 is already  too much for 32 bit FP numbers.
>
> ...
>
> I think I'll be able to find other FFT libraries which work with non-power of 2 numbers 64 bit floats.

> To emphasize the point(s).
>
> If/when we perform spectral analysis, we are interested in two (three) things:
>
> 1) presence/absence of spectral components;
>
> 2) magnitude of spectral components;
>
> 3) maybe phase of spectral components.
>
> The issue I described causes lack of resolution when looking for spectral components. I.e. not sufficient number of bits in DFT prevents properly detecting spectral components. With 32 bit FP numbers signal itself has ~23 bit resolution, but frequency domain resolution is less than 12 bits.

> And the same kind of test for 64 bit FP numbers:
>
> "
>
> julia> let; x::Float64 = 1.0 / (2^27); cos(x) end
> 1.0
>
> julia> let; x::Float64 = 1.0 / (2^26); cos(x) end
> 0.9999999999999999
>
> julia> 2.0^26 / 44100
> 1521.742947845805
>
> ",
>
> i.e. 27 bits is too much. Which means that for audio applications 64 bit floats are quite OK as long as buffer size is less than 1522 seconds @ 44100Hz sample rate.

> And back to the issue whether 'pffft' supports 64 bit FP numbers - in https://github.com/cpuimage/pffft/blob/master/fftpack.h one can see (line #52) :
>
> "
>
> // just define FFTPACK_DOUBLE_PRECISION if you want to build it as a double precision fft
>
> #ifndef FFTPACK_DOUBLE_PRECISION
>   typedef float fftpack_real;
>   typedef int   fftpack_int;
> #else
>   typedef double fftpack_real;
>   typedef int    fftpack_int;
> #endif
>
> ".
>
> --Sergei.

On Sat, 29 Nov 2025 at 14:37, Martin Guy <[email protected]> wrote:
>
> Thanks, that's all clear now.
>
> > OK as long as buffer size is less than 1522 seconds @ 44100Hz sample rate
>
> I think we can guarantee that!
>
> > // just define FFTPACK_DOUBLE_PRECISION
>
> Good. SoX's internal "fft4g" routine has the same thing:
>
> #ifdef FFT4G_FLOAT
>   #define double float
>   #define sin   sinf
>   #define cos   cosf
>   #define atan  atanf
>
>   #define cdft  lsx_cdft_f
>   #define rdft  lsx_rdft_f
> ...
>
> which means some jiggery pokery to compile it twice into different object files
> to get both versions (something I was already contemplating for fft4g) and
> that means including the pffft source in SoX' source tree to be able to do this,
> something that would be necessary anyway as it's not included in many distros,
> only ArchLinux, MacPorts and FreeBSD Ports.
> https://repology.org/project/pffft/versions
>
> Thanks for your exhaustive analysis of the necessary precision in DFTs.
> I will be careful with them.
>
>     M