[ENH] Complete last cycle of bursts

[ryanhammonds]

For bursty signals, the last cycle of each burst was incomplete by one sample. This updates fixes this by checking if the next cycle is not oscillating, and if so, adding the last sample to complete the cycle.

[ryanhammonds]

This issue extends to non-bursty / continuous oscillations as well. In addition to the last sample not completing a full cycle, there is a second issue: If a cycle is expected to complete in n_samples, it will actually complete in n_samples + 1 (i.e. all the blue points, plus the first orange point below), and may result in a tiny frequency shift. So the two related issues are:

1. Full cycles complete in n + 1 samples, when n is intended
2. Tiling will always leave out the last sample

import numpy as np
import matplotlib.pyplot as plt
from neurodsp.sim import sim_oscillation

n_seconds = 1
fs = 60
freq = 2

ys = sim_oscillation(n_seconds, fs, freq)
xs = np.arange(len(ys))

plt.axhline(0, color='k', ls='--', alpha=.25)
plt.plot(xs[:30], ys[:30], ls='', marker='.')
plt.plot(xs[30:], ys[30:], ls='', marker='.')

I don't think there is a clean solution to this problem using tiling. The alternative fix for sine waves is simulate a full time array to pass into np.sin. This could also be extend to bursty oscillations, fixing the last sample issue as well. However, this wouldn't work for non-sine waves and would require us to write quite a bit of code to update simulation functions for all cycles to work similar to np.sin.

import numpy as np
import matplotlib.pyplot as plt

n_cycles = 2
n_samples = 60

times = np.linspace(0, n_cycles * 2 * np.pi, n_samples)
sig = np.sin(times)

plt.plot(sig, ls='', marker='.')
plt.axhline(0, color='k', alpha=.25, ls='--')

tldr : There are two issues with tiling and would take a lot of work to fully fix both issues. At the end of the day this issue are probably relatively minor. If we decide to rework tiling, it probably should be post 2.2 release.

[TomDonoghue]

These are very good points / notes. I think the main thing here ("missing" last sample) is basically a property of the tiling approach. It's less clear to me to it's a problem per se. In order to tile, we need to skip the last sample value - if it didn't then at the "attachment" point we would have two adjacent points at the same y-value. So the sims as they are as intended, I think, the important question is if there is an issue with this.

For continuous tiling, I think this is basically fine - there is never actually a "missing last sample". This does have the slightly odd "missing last sample" with bursts. I guess I'm not super clear this is a real problem though. We could add the "last cycle" fix - but does it really change anything, in practice? I guess I would be curious to see if one or other have a "better" power spectrum, or something. If not, I'm not sure this really matters.

The "off by one" frequency aspect is interesting. I'm not sure there's a problem in practice, in part because we define the cycles in terms of the estimated time length of the cycle, and I'm not totally we do end up with an off by one-error. And I think that last sample shouldn't be counted in the "cycle", as it's effectively part of the next cycle? Otherwise the next cycle would be "missing" the first sample? Even if there is an off by one sample issue here - I think for any reasonable combination of frequency + sampling rate, this effect should be quite negligible.

More generally, I don't think there's an "uh oh, this is a problem" issue here - unless I'm missing something? I think some of the benefits of tiling (different cycles + burst elements) make it worth keeping. For some contexts, it probably does make sense to simulate a continuous oscillation with np.sin, but I think I would suggest that's a possible extension / addition, rather than a replacement. Thoughts?

[ryanhammonds]

I think the "missing last sample" does apply to continuous simulations. The first plot above is a continuous simulation, and the last sample doesn't end on zero. However, for bursty simulation, this is occurs at the end of every burst, whereas in continuous simulations this only happens once at the end of the signal.

The off by one issue is present because cycles are simulated to not fully complete to allow for tiling. Instead, cycles fully complete on the first sample of the next cycle. If a cycle is simulated at 10hz and 1000 fs, a cycle is simulated with 100 samples. But cycles will start/end at zero in 101 samples. This is negligible at a high enough sampling rate.

A better example of this issue with incomplete cycles iswhen an odd number of samples per cycle are required (#223):

import numpy as np
import matplotlib.pyplot as plt
from neurodsp.sim import sim_oscillation
from neurodsp.spectral import compute_spectrum
from neurodsp.utils.norm import normalize_sig

n_seconds = 10
fs = 1000
freq = 9

# Tiling
sig = sim_oscillation(n_seconds, fs, freq)

# np.sin
times = np.arange(0, n_seconds, 1/fs)
sig_np_sin = np.sin(freq * times * 2 * np.pi)
sig_np_sin = normalize_sig(sig_np_sin, mean=0, variance=1)

# PSD
freqs, powers = compute_spectrum(sig, fs)
freqs_np_sin, powers_np_sin = compute_spectrum(sig_np_sin, fs)

# Plot
plt.figure(figsize=(8, 5))
plt.loglog(freqs, powers, label='Tiling')
plt.loglog(freqs_np_sin, powers_np_sin, label='np.sin')
plt.legend();

For most cases, this small frequency mismatch isn't a huge deal. And you're right that there are benefits to tiling in cases of mixing cycles. I think we could potentially develop an extension for continuous / non-tiled simulations, with an optional kwarg in simulations to choose from either 'tiling' or 'continuous' - but this probably isn't a huge priority and can be explored after the next release.