ezmsg-org / ezmsg-sigproc

Timeseries signal processing implementations for ezmsg
MIT License
0 stars 0 forks source link

An error occurred when using wavelet ezmsg.sigproc.wavelets.cwt #41

Open RuolingWu opened 2 weeks ago

RuolingWu commented 2 weeks ago

I came across the following error when I used ezmsg.sigproc.wavelets.CWT to do wavelet transform:

2024-11-03 09:59:44.909 - pid: 36349 - TaskThread - WARNING - windowing: windowing is non-deterministic with zero_pad_until='input' as it depends on the size of the first input. We recommend using 'shift' when window_shift is float-valued.
2024-11-03 09:59:45.005 - pid: 36349 - TaskThread - INFO - on_signal: Traceback (most recent call last):
File "/Users/rwu/anaconda3/envs/umcuezmsg/lib/python3.12/site-packages/ezmsg/sigproc/base.py", line 32, in on_signal
ret = self.STATE.gen.send(message)
^^^^^^^^^^^^^^^^^^^^^^^^^^^^
File "/Users/rwu/anaconda3/envs/umcuezmsg/lib/python3.12/site-packages/ezmsg/sigproc/wavelets.py", line 108, in cwt
template = AxisArray(
^^^^^^^^^^
File "<string>", line 7, in init
File "/Users/rwu/anaconda3/envs/umcuezmsg/lib/python3.12/site-packages/ezmsg/util/messages/axisarray.py", line 72, in post_init
raise ValueError("dims contains repeated dim names")
ValueError: dims contains repeated dim names

The error occurred in line 109-119, wavelets.py (https://github.com/ezmsg-org/ezmsg-sigproc/blob/main/src/ezmsg/sigproc/wavelets.py)

template = AxisArray(
                np.zeros(
                    dummy_shape, dtype=dt_cplx if wavelet.complex_cwt else dt_data
                ),
                dims=msg_in.dims[:ax_idx] + msg_in.dims[ax_idx + 1 :] + ["freq", axis],
                axes={
                    **msg_in.axes,
                    "freq": AxisArray.Axis("Hz", offset=freqs[0], gain=fstep),
                },
            )

I checked the value of dims=msg_in.dims[:ax_idx] + msg_in.dims[ax_idx + 1 :] + ["freq", axis]. The result I got is ['freq', ‘ch', ‘freq', ‘time']. It seems that the two ‘freq’ caused the error.

And here’s the node I created for the wavelet transform:

"GET_WAVELETS": CWT(scales=scales, wavelet="cmor3.0-0.5", min_phase=MinPhaseMode.HOMOMORPHIC, axis='time')

Environment information: Python version: 3.12 Conda version: conda 24.4.0 Os: MacOs

cboulay commented 2 weeks ago

Hi @RuolingWu , I hope you don't mind that I edited your post to make the formatting a little easier to read -- I didn't change any content other than add some backticks.

I would need to know what your pipeline is preceding the CWT node. It seems like the signal coming into the node already has a freq axis which maybe suggests a misunderstanding of what the CWT node is intended for.

RuolingWu commented 2 weeks ago

Hi Chad,

Thanks for changing the format. I checked my code and realised I forgot to delete the node for calculating band power. But after I deleted that node, I got a new error in CWT: 2024-11-04 22:00:59.455 - pid: 51795 - TaskThread - WARNING - windowing: windowing is non-deterministic with zero_pad_until='input' as it depends on the size of the first input. We recommend using 'shift' when window_shift is float-valued. 2024-11-04 22:00:59.550 - pid: 51795 - TaskThread - INFO - on_signal: Traceback (most recent call last): File "/Users/rwu/anaconda3/envs/umcuezmsg/lib/python3.12/site-packages/ezmsg/sigproc/base.py", line 32, in on_signal ret = self.STATE.gen.send(message) ^^^^^^^^^^^^^^^^^^^^^^^^^^^^ File "/Users/rwu/anaconda3/envs/umcuezmsg/lib/python3.12/site-packages/ezmsg/sigproc/wavelets.py", line 122, in cwt conv_msg = fbgen.send(msg_in) ^^^^^^^^^^^^^^^^^^ File "/Users/rwu/anaconda3/envs/umcuezmsg/lib/python3.12/site-packages/ezmsg/sigproc/filterbank.py", line 107, in filterbank sps.minimum_phase( File "/Users/rwu/anaconda3/envs/umcuezmsg/lib/python3.12/site-packages/scipy/signal/_fir_filter_design.py", line 1230, in minimum_phase raise ValueError('Complex filters not supported') ValueError: Complex filters not supported

I checked the Linear-phase FIR filter coefficients (h) in scipy.signal.minimum_phase.minimum_phase(h, method='homomorphic'), which is a complex array array([2.76796018e-12-2.71181886e-27j])

And here is the entire pipeline:

   frequencies = range(60, 120)
   scales = pywt.frequency2scale('cmor3.0-0.5', frequencies)
    comps = {
        "GET_SOURCE": NeoIteratorUnit(
            filepath=base_path / playback_file,
            chunk_dur=0.1,
        ),
        "SELECT_STREAM": FilterOnKey("Signals"), # checked
        "SELECT_CHANNELS": Slicer(selection="0, 1, 2", axis='ch'), # checked
        "FILTER_LINENOISE": ButterworthFilter(axis='time', order=5, cuton=105, cutoff=95),# checked
        "APPLY_CAR": CommonRereference(mode="mean", axis='ch'), # checked
        "GET_WAVELETS": CWT(scales=scales, wavelet="cmor3.0-0.5", min_phase=MinPhaseMode.HOMOMORPHIC, axis='time'),
        "APPLY_MEAN": RangedAggregate(axis='freq', operation=AggregationFunction.MEAN),
        "APPLY_ABS": Abs(),
        "DOWNSAMPLE": Downsample(axis='time', factor=20),
        "APPLY_LOG": Log(),
        "RELEASE_SINK": Rerun() if has_rerun else DebugLog(name="TEST", max_length=80), # checked
    }
    conns = [
        (comps["GET_SOURCE"].OUTPUT_SIGNAL, comps["SELECT_STREAM"].INPUT_SIGNAL),
        (comps["SELECT_STREAM"].OUTPUT_SIGNAL, comps["SELECT_CHANNELS"].INPUT_SIGNAL),
        (comps["SELECT_CHANNELS"].OUTPUT_SIGNAL, comps["FILTER_LINENOISE"].INPUT_SIGNAL),
        (comps["FILTER_LINENOISE"].OUTPUT_SIGNAL, comps["APPLY_CAR"].INPUT_SIGNAL),
        (comps["APPLY_CAR"].OUTPUT_SIGNAL, comps["GET_WAVELETS"].INPUT_SIGNAL),
        (comps["GET_WAVELETS"].OUTPUT_SIGNAL, comps["APPLY_MEAN"].INPUT_SIGNAL)
    ]
cboulay commented 2 weeks ago

I checked the Linear-phase FIR filter coefficients (h) in scipy.signal.minimum_phase.minimum_phase(h, method='homomorphic')

Sorry, what do you mean by "I checked"? Does this mean that you tried running the exact same line of code and it worked for you, so it's only ezmsg that's failing?

I've encountered this problem myself and I was forced to use a different wavelet type when doing a minphase CWT. (Note: minphase may or may not be important depending on your task and how much you care about your lower frequency wavelets' delay.)

RuolingWu commented 2 weeks ago

Sorry for the confusion. I just put a break point in that line to see the actual value of h. But thank you for your suggestion! And I will run another test on scipy.signal.minimum_phase.minimum_phase(h, method='homomorphic') to figure out if this is an ezmsg problem or a scipy problem.