jonathanlilly / jLab

A Matlab toolbox for big data analysis, signal processing, mapping, and oceanographic applications.
Other
107 stars 35 forks source link

Analytic Morse may be ill-defined for even lengths #13

Open OverLordGoldDragon opened 3 years ago

OverLordGoldDragon commented 3 years ago

Thank you for open-sourcing this.

Looking at definition of analytic signals, I've noted the transform from real to analytic to be s(t) -> s(t) + js^(t), i.e. the real part is preserved. For even lengths, however, as per Hilbert transform, this requires us to halve the dc and Nyquist bins (or double others):

N = length(x);  # try any even-length x
xf = fft(x);
xfa = horzcat(xf(1), 2*xf(2:N/2), xf(N/2+1), zeros(1, N/2-1));
xa = ifft(xfa);  # real(xa) == real(x)

Here is tail behavior comparison of original and adjusted time-domain wavelets (psi):

image

Adjusted properly decays to zero; this enables stable and accurate computation of wavelet time resolution for low scales (otherwise t**2 in variance integral won't bound).

However, another point of concern: there is no Nyquist bin: psif(N/2+1) == 0. Whether the bin belongs to positives or seems arbitrary, but in terms of basis functions the two are equivalent so one might say "both", where again it makes sense to halve it to "zero the negative half". More importantly the code above fails to yield real(xa) == real(x) without it. Above plots were produced without Nyquist bin; I've yet to investigate behavior with vs without, will update.

OverLordGoldDragon commented 3 years ago

Adding Nyquist also swaps symmetries of time-domain wavelet, peaking real component at center as probability distribution; comparison with and without Nyquist bin, both "adjusted" with halving:

image

This breaks symmetry in the standard sense, as even-length sequences can't have a discrete center, but makes it DFT-symmetric (or "periodic"), which has more desired properties in spectral analysis (stronger decay). Also applies to odd N.

And finally comparing all adjustments vs original:

image

So at least for time-domain representation the adjustments seem desired; unsure about halving in CWT, but adding the Nyquist bin does come with increased "dynamic range" of frequencies (greater max center frequency), and the Nyquist bin is uniquely consistent for every sequence length (i.e. (-1)^[0, 1, ...]). For complex rotation to center time-domain wavelet with Nyquist bin I used *(-1).^[0:N-1].

OverLordGoldDragon commented 3 years ago

Time resolution integrals for original vs fully-adjusted; results still differ significantly from theoretically-predicted values (i.e. compute at greater scale with proper decay, then divide):

image

Note that t is defined a little differently (important) to accommodate differing symmetries in each case. Also I used a high center frequency or low scale in all my examples (i.e. where the issue is relevant), with the frequency-domain wavelet trimmed at negatives rather than decaying smoothly.

OverLordGoldDragon commented 3 years ago

Found concrete benefits to CWT in implementing both adjustments:

  1. Without the Nyquist bin, CWT remains completely blind to frequency components at Nyquist frequency, which is undesired.
  2. The DFT forward-transforms the (input's) Nyquist bin at twice the value of any other bin (take pure tones of same amplitude), so not halving the Nyquist bin yields doubled CWT coefficients for associated rows relative to others.
  3. Not halving results in failure to capture amplitude modulation near Nyquist

All illustrated below: x = cos(2*pi*1024*t) + cos(2*pi*1023*t); t=linspace(0, 2047, 2048)

image

I did not, however, find any time resolution benefits, testing on rapid sine-FM signals, the CWT's were exactly same except for exact expected difference of half a Nyquist bin, suggesting behavior to change only in frequency-domain interactions.

I've used Morlet w/ ssqueezepy.cwt, but key behavior will be same with Morse. (p.s. I'm also porting Morse into ssqueezepy).

OverLordGoldDragon commented 3 years ago

Another artifact evident from a heatmap of time-domain wavelet vs scales, clipped to show small values:

image

Effect diminishes only after both, including Nyquist + halving:

image

Note this is with Morlet which is poorly-behaved for scales near 1, hence the absolute value not fully compacting.