Signal processing functionality in DWF

[Guest]

Hi,

I am currently working on some examples for pydwf that show how to use the FDwfSpectrum functionality implemented in the DWF library.

To make those examples I am thoroughly testing the three FDwfSpectrum functions, and I'm seeing several cases where I don't think the behavior is quite correct. I will detail my findings in this post and some followups, in the hope that the behavior observed will be considered for fixing, and for feedback.

Starting with the simplest function: FDwfSpectrumFFT.

First, the way the function returns its result (as separate phase and bin magnitude arrays, rather than real/imaginary components of the complex number results) is quite non-standard. No mainstream FFT implementation that I am aware of does this, and I've seen about a dozen.

Second, the scaling of the FFT (as returned in the bin magnitude vector) deviates from the convention that essentially all FFT implementations (except Mathematica) use nowadays; it is smaller by a factor (2 / n), where n is the length of the vector being FFT'ed. This unusual scaling choice is undocumented.

Third, the function gives an incorrect result of +infinity or -infinity when doing an FFT on a 1-element vector. One could argue that this is not an interesting case, but I see no reason why it shouldn't return the correct result. If this is a "won't fix", at least the documentation should indicate that the function only works correctly for n >= 2.

Lastly, the performance is quite bad compared to a properly optimized function like numpy's rfft function. The graph attached shows how much slower FDwfSpectrumFFT is compared to the numpy function, ranging from 2.5 to 20 times for different FFT sizes. The Numpy function also has the advantage of correctly handling input sizes that are not a power of two, unlike the FDwfSpectrumFFT function.

(My personal opinion is that the FDwfSpectrum functionality should have never made it into DWF, since (1) it simply has no place in an instrument-control library and constitutes a clear-cut example of feature bloat in an already huge library, and (2) dedicated, specialized signal processing libraries exist that handily outperform the FDwfSpectrum implementations. But that ship has sailed, I'm afraid.)

[attila]

Hi @reddish

Like most features and functions, these were added at the request of the users.
C-API can be used with other languages where FFT functions/libraries are not available.

[Guest]

> C-API can be used with other languages where FFT functions/libraries are not available.

I can think of many other things not related to your devices for which the same is true. You have to draw a line somewhere, and the line you guys have drawn here is seriously weird to me. Anyway, my opinion on this is not very relevant; I don't expect you to revert this decision, so we can drop that.

What remains is the question whether you will fix the n==1 bug, and whether you'd consider bringing the scaling in line with the de facto standard as used by all other FFT libraries.

[attila]

Hi @reddish

Yes, it is fixed. 
Thanks for noticing.

[Guest]

Hi @attila

I assume you mean the n==1 bug is fixed, but you're not going to change the scaling? Or are you fixing that as well?

With regard to scaling: I tested the FDwfSpectrumTransform function. It appears to be working correctly but like the FDwfSpectrumFFT function, its scaling deviates from established practice. I compared to Matlab's czt function and Python's scipy.signal.czt functions, for example.

To fix the scaling, i.e., bring it in line with the conventions of Matlab and Python, the "rgdBin" values returned from the function currently need to be divided by (cdData / 2). I urge you to fix or at least document this, as it may trip up users.

What remains is a test of FDwfSpectrumWindow, which still has a couple of issues. I'll gather my findings and post those here.

[Guest]

Hi @attila, I'm still interested to hear what, if anything, you will do with the scaling issues identified in the FDwfSpectrumFFT and FDwfSpectrumTransform functions.

[Guest]

Okay, the last function in the FDswfSpectrum family is the FDswfSpectrumWindow function. I tested the 3.20.1 release and I have the following remarks.

(1) You scale all returned windows such that the sum of the elements in the window is equal to the length of the window. Which is fine, except for windows that come out as all-zero. I think it would be a good idea to document this scaling convention and to document what happens in case a window is requested that is zero everywhere — because then the scaling convention breaks down.

(2) The noise-equivalent bandwidth value that you return is equal to n * sum(w**2) / sum(w)**2, where w is the discrete window scaled according to the convention stated in the previous point. I think it would be useful to document this fact. It will also be useful to document what the value is for an all-zero window because in that case the calculation above once again breaks down.

(3) The behavior of all windows for small window sizes (cdWin parameter) is sometimes weird. I could list all instances where it returns NaNs for the different windows, but it may be best to test and fix this by yourself. In most (perhaps all) cases the NaN window values returned are caused by the convention stated under point (1) above.

(4) The behavior of the Hann and Cosine windows are almost identical (up to numerical position). I think this is not intentional. In SciPy, for example, the "hann" window looks like a cosine from [-pi .. +pi], whereas the "cosine" window looks like a cosine evaluated from [-0.5*pi .. +0.5*pi]. I think that's also the intention of your implementation, but that's not what it does currently.

(5) Your Flattop window is not equal to the similarly-named window in Matlab, Octave, and SciPy. The problem is that "flat-top" windows are not canonically defined in literature and software; Matlab and SciPy follow one definition, Octave follows another, and DWF follows yet another. I think it would be useful to document the precise function your flattop-window definition follows, or even better: change the code to comply with the definition of Matlab/SciPy.

(6) Your implementation of the Kaiser window is incorrect, see the image below. Feel free to copy the correct implementation (and a definition of the Bessel I0 function on which it depends) given here: https://github.com/sidneycadot/WindowFunctions/blob/master/c99/window_functions.c; but make sure to scale in accordance with your scaling convention.

Cheers, Sidney

Edited August 10, 2023 by reddish

[attila]

Hi @reddish

1,2,3,6 should be fixed in the latest version:

Have to look at 4,5 later

Thank you for your posts.

[Guest]

Hi @attila

I tested with the latest beta. Indeed most problems seem to be addressed.

NEBW value is returned as "1" in case an all-zero window is returned, which is the most sane choice I suppose. Might be useful to document that.

The "Hann" window now returns [0 2] for n==2, this seems to be an error.

The "Cosine" window now returns [0 0] for n==2, this seems to be an error.

(Cosine window is still identical to Hann window, otherwise -- my previous point (4)).

Thanks for fixing the issues.

Cheers, Sidney


EDIT: one more thing, I see you're adding public filtering functions to the API -- great. But the enum types and functions all use the term "Fiir" rather than "Fir". I have always and only every seen these filter types referred to as FIR filters (FIR for "finite impulse response"), so that seems to be more standard terminology. I hope you will reconsider.

Edited August 11, 2023 by reddish