MideTechnology / endaq-python

A comprehensive, user-centric Python API for working with enDAQ data and devices
MIT License
25 stars 12 forks source link

Bug with loglog_linear_approx #95

Open S-Hanly opened 2 years ago

S-Hanly commented 2 years ago

@CrepeGoat I tried to use the new loglog_linear_approx function but first hit an error when I inputted a PSD with a 0 Hz value (couldn't divide by 0). That should be removed within the function, the PSD I used was one I got from the welch function.

The next error I'm not so sure as to the cause. Here's a colab: https://colab.research.google.com/drive/1VD96qmdCAw0PMUovcCgC5BxFXmFreGD7?usp=sharing

---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
<ipython-input-7-18cdedbfe22e> in <module>()
----> 1 break_points = endaq.calc.psd.loglog_linear_approx(psd[0.5:], [1, 30,35,45,190,210,250,1375,2000], window=1, freqs_out='knots')
      2 break_points

3 frames
/usr/local/lib/python3.7/dist-packages/endaq/calc/psd.py in loglog_linear_approx(df, knots, window, freqs_out)
    254     tcks = [
    255         scipy.interpolate.splrep(index_log, dlog, k=1, t=knots_log)
--> 256         for dlog in data_log.T
    257     ]
    258 

/usr/local/lib/python3.7/dist-packages/endaq/calc/psd.py in <listcomp>(.0)
    254     tcks = [
    255         scipy.interpolate.splrep(index_log, dlog, k=1, t=knots_log)
--> 256         for dlog in data_log.T
    257     ]
    258 

/usr/local/lib/python3.7/dist-packages/scipy/interpolate/fitpack.py in splrep(x, y, w, xb, xe, k, task, s, t, full_output, per, quiet)
    289 
    290     """
--> 291     res = _impl.splrep(x, y, w, xb, xe, k, task, s, t, full_output, per, quiet)
    292     return res
    293 

/usr/local/lib/python3.7/dist-packages/scipy/interpolate/_fitpack_impl.py in splrep(x, y, w, xb, xe, k, task, s, t, full_output, per, quiet)
    510         else:
    511             try:
--> 512                 raise _iermess[ier][1](_iermess[ier][0])
    513             except KeyError as e:
    514                 raise _iermess['unknown'][1](_iermess['unknown'][0]) from e

ValueError: Error on input data
CrepeGoat commented 2 years ago

Okay so those two bugs: 1) re: the 0-frequency error, that makes sense because it has to calculate the log of the frequency values, and getting the log of 0 gives you a NaN. having said that you're right in that the function should handle that very common case a little more gracefully. 2) re: the other error, it looks like it broke because one of you knots (1) was lower than the lowest frequency sample (1.000015); it works fine if you change the knots from [1, ...] to [2, ...]. but even though the inputs were incorrect, the function should provide a clearer error message.

all-in-all it looks like this function just needs some better input sanitation. thanks @shanlyMIDE! let me know if you find anything else.

S-Hanly commented 2 years ago

Ok nice, I see what I did there! I updated the colab to fix that: https://colab.research.google.com/drive/1VD96qmdCAw0PMUovcCgC5BxFXmFreGD7?usp=sharing

Here's the PSD which matches up well with some break points I selected: newplot - 2021-12-15T155149 463

It does a good job getting the cum RMS but seems to miss things around peaky resonances. If I add more breakpoints there it gets better. Here it is with the same breakpoints in the PSD. newplot - 2021-12-15T155058 168

But now here is when I add more breakpoints near that prominent ~35 Hz point. newplot - 2021-12-15T155237 490

I think this may deserve a blog!

CrepeGoat commented 2 years ago

hmm so I see in the cumsum plots that the approximations are missing some energy content, particularly in the X/Y subchannels... that's a little disheartening :\ I'm glad it's (kind of) working at least!

Re: the missing peaks, it might do a better job of characterizing the data if the "smoothness" parameter was reduced? maybe I should expose that as a function parameter... 🤔

CrepeGoat commented 2 years ago

the loglog_linear_approx function was removed in #113 to avoid presenting users with this issue.. The function should be later added once a more accurate solution is provided.