{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n# Cycle detection from IMFs\nHere we will use the 'cycle' submodule of EMD to identify and analyse individual cycles of an oscillatory signal\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Simulating a noisy signal\nFirstly we will import emd and simulate a signal\n\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"import emd\nimport numpy as np\nimport matplotlib.pyplot as plt\n\n# Define and simulate a simple signal\npeak_freq = 15\nsample_rate = 256\nseconds = 10\nnoise_std = .4\nx = emd.simulate.ar_oscillator(peak_freq, sample_rate, seconds,\n noise_std=noise_std, random_seed=42, r=.96)[:, 0]\nt = np.linspace(0, seconds, seconds*sample_rate)\n\n# Plot the first 5 seconds of data\nplt.figure(figsize=(10, 2))\nplt.plot(t[:sample_rate*3], x[:sample_rate*3], 'k')\n\n# sphinx_gallery_thumbnail_number = 5"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Extract IMFs & find cycles\nWe next run a mask sift with the default parameters to isolate the 15Hz\noscillation. There is only one clear oscillatory signal in this simulation.\nThis is extracted in IMF-3 whilst the remaining IMFs contain low-amplitude\nnoise.\n\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"# Run a mask sift\nimf = emd.sift.mask_sift(x)\n\n# Visualise the IMFs\nemd.plotting.plot_imfs(imf[:sample_rate*4, :])"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Next, we want to identify single cycles of any oscillations that are present\nin our IMFs. There are many ways to do this depending on the signal and\nsignal-features of interest. The EMD package provides a method which extract\ncycle locations based on the instantaneous phase of an IMF. In its simplest\nform this will detect successive cycles within the IMF, though we can also\nrun some additional checks to reject 'bad' cycles or specify time-periods to\nexclude from the cycle detection.\n\nWill will run through some examples of these in the next couple of sections.\nFirst, we compute the instantaneous phase of our IMFs using the\n``frequency_transform`` function before detecting cycle indices from the IP\nusing ``get_cycle_vector``.\n\nThe detection is based on finding large jumps in the instantaneous phase of\neach IMF. By default, we will consider any phase jump greater than 1.5*pi as\na boundary between two cycles. This can be customised using the\n``phase_step`` argument in ``get_cycle_vector``.\n\nWe can optionally run a set of tests on the detected cycles to remove 'bad'\ncycles from the analysis at an early stage. We will go into more detail into\nthis later, for now we will run the function to get ``all_cycles``.\n\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"# Compute frequency domain features using the normalised-Hilbert transform\nIP, IF, IA = emd.spectra.frequency_transform(imf, sample_rate, 'nht')\n\n# Extract cycle locations\nall_cycles = emd.cycles.get_cycle_vector(IP, return_good=False)\n\n# ``all_cycles`` is an array of the same size as the input instantaneous phase.\n# Each row contains a vector of itegers indexing the location of successive\n# cycles for that IMF.\n\nprint('Input IMF shape is - {0}'.format(IP.shape))\nprint('Input all_cycles shape is - {0}'.format(all_cycles.shape))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"For each IMF, ``all_cycles`` stores either a zero or an integer greater than\nzero for each time sample. A value zero indicates that no cycle is occurring\nat that time (perhaps as it has been excluded by our cycle quality checks in\nthe next section) whilst a non-zero value indicates that that time-sample\nbelongs to a specific cycle.\n\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"# Firstly, cycles are detected by looking for phase 'jumps' where the phase\n# resets between cycles. The default threshold for a jump is 1.5*pi. Let's\n# take a look at this in IMF-3\n\nplt.figure(figsize=(8, 6))\nplt.subplots_adjust(hspace=0.3)\nplt.subplot(311)\nplt.plot(t[:sample_rate*2], imf[:sample_rate*2, 2])\nplt.gca().set_xticklabels([])\nplt.title('IMF-3')\nplt.subplot(312)\nplt.plot(t[:sample_rate*2], IP[:sample_rate*2, 2])\nplt.title('IMF-3 Instantaneous Phase')\nplt.ylabel('Radians')\nplt.gca().set_xticklabels([])\nplt.subplot(313)\nplt.plot(t[1:sample_rate*2], np.abs(np.diff(IP[:sample_rate*2, 2])))\nplt.plot((0, 2), (1.5*np.pi, 1.5*np.pi), 'k:')\nplt.xlabel('Time (seconds)')\nplt.title('IMF-3 Instantaneous Phase Abs-Differential')\nplt.legend(['IMF-3 IP Differential', 'Jump threshold'], loc='upper right')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We can see that the large phase jumps occur at the ascending zero-crossing of\neach cycle. In a clear signal, these are very simple to detect using a blunt\nthreshold.\n\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## What makes a 'good' cycle?\nThere are many methods for detecting oscillatory cycles within a signal. Here\nwe provide one approach for identifying whether a signal contains clear and\ninterpretable cycles based on analysis of its instantaneous phase. This\nprocess takes the cycles detected by the phase jumps and runs four\nadditional checks.\n\nWe define a 'good' cycle as one with:\n\n1. A strictly positively increasing phase\n2. A phase starting within phase_step of zero ie the lowest value of IP must be less than phase_step\n3. A phase ending within phase_step of 2pi the highest value of IP must be between 2pi and 2pi-phase_step\n4. A set of 4 unique control points (ascending zero, peak, descending zero & trough)\n\nLets take a look at these checks in IMF-3. Firstly, we run test 1:\n\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"# We use the unwrapped phased so we don't have to worry about jumps between cycles\nunwrapped_phase = np.diff(np.unwrap(IP[:sample_rate*2, 2]))\n\n# Plot the differential of the unwrapped phasee\nplt.figure(figsize=(8, 4))\nplt.subplot(211)\nplt.plot(t[:sample_rate*2], IP[:sample_rate*2, 2])\nplt.legend(['IMF-3 Instantaneous Phase'])\nplt.ylabel('Radians')\nplt.title('Test-1: Check phase is strictly increasing')\nplt.subplot(212)\nplt.plot(t[1:sample_rate*2], unwrapped_phase)\nplt.plot((0, 2), (0, 0), 'k:')\nplt.ylim(-.2, .4)\nplt.legend(['IMF-3 Instantaneous Phase Differential'])"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We can see that the instantaneous phase of most cycles is positive throughout\nthe cycle. Only one cycle (around 1.6 seconds into the simulation) has\nnegative values which correspond to a reversal in the normal wrapped IP.\n\nThe second test looks to make sure that each cycles phase covers the whole\n2pi range. If the phase doesn't reach these limits it indicates that a phase\njump occurred early or late in the cycle (for instance we might have a peak\nwhich is below zero in the raw time-course) leaving an incomplete\noscillation.\n\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"plt.figure(figsize=(8, 4))\nplt.title('Tests 2+3: Check phase covers the full 2pi range')\nplt.plot(t[:sample_rate*2], IP[:sample_rate*2, 2], label='IP')\nplt.plot((0, 2), (0, 0), label='0')\nplt.plot((0, 2), (np.pi*2, np.pi*2), label='2pi')\n\nplt.plot((0, 2), (np.pi*2-np.pi/12, np.pi*2-np.pi/12), ':', label='Upper Thresh')\nplt.plot((0, 2), (np.pi/12, np.pi/12), ':', label='Lower Thresh')\nplt.legend()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We can see that most cycles have instantaneous phase values crossing both the\nupper and lower threshold. Only the first and last cycles in these segments\nare missing these thresholds (as they are incomplete cycles cutt-off at the\nedges of this segment)\n\nFinally, we check that we can detect a complete set of control points from\neach cycle. The control points are the peak, trough, ascending zero-crossing\nand descending zero-crossing. These can be computed from the IMF and a cycles\nvector using ``emd.cycles.get_control_points``\n\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"ctrl = emd.cycles.get_control_points(imf[:, 2], all_cycles[:, 2])\n\n# Define some marker styles and legend labels for the control points.\nmarkers = ['og', '^b', 'oc', 'vb', 'or']\nlabel = ['Asc-Start', 'Peak', 'Desc', 'Trough', 'Asc-End']\n\n# Plot the first 10 cycles with control points\nncycles = 20\nstart = 0\n\nplt.figure()\nplt.plot(111)\nplt.title('Test 4: Control points')\nfor ii in range(ncycles):\n print('Cycle {0:2d} - {1}'.format(ii, ctrl[ii, :]))\n cycle = imf[all_cycles[:, 2] == ii, 2]\n plt.plot(np.arange(len(cycle))+start, cycle, 'k', label='Cycle')\n for jj in range(5):\n if np.isfinite(ctrl[ii, jj]):\n plt.plot(ctrl[ii, jj]+start, cycle[int(ctrl[ii, jj])], markers[jj], label=label[jj])\n start += len(cycle)\n\n # Only plot the legend for the first cycle\n if ii == 1:\n plt.legend()\nplt.ylim(-400, 400)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Most of these cycles have the full set of control points present. Only ones\ncycle (cycle-20 - close to the end) is missing an indicator for its\npeak or trough. This is as a distortion in the cycle means that there are two\npeaks and troughs present. In this case, ``get_control_points`` will return a\n``np.nan`` as the value for that peak.\n\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We run these checks together automatically by setting the ``return_good``\noption to ``True``. This is also the default option in the code. Here we run\ncycle detection with the quality checks on and look at the first four seconds\nof signal.\n\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"good_cycles = emd.cycles.get_cycle_vector(IP, return_good=True, phase_step=np.pi)\n\nplt.figure(figsize=(10, 8))\nplt.subplots_adjust(hspace=.3)\nplt.subplot(311)\nplt.plot(t[:sample_rate*4], imf[:sample_rate*4, 2], 'k')\nplt.gca().set_xticklabels([])\nplt.title('IMF-3')\nplt.subplot(312)\nplt.plot(t[:sample_rate*4], IP[:sample_rate*4, 2], 'b')\nplt.gca().set_xticklabels([])\nplt.plot((0, 4), (0, 0), label='0')\nplt.plot((0, 4), (np.pi*2, np.pi*2), label='2pi')\nplt.plot((0, 4), (np.pi*2-np.pi/12, np.pi*2-np.pi/12), ':', label='Upper Thresh')\nplt.plot((0, 4), (np.pi/12, np.pi/12), ':', label='Lower Thresh')\nplt.title('Instantanous Phase')\nplt.ylabel('Radians')\nplt.subplot(313)\nplt.plot(t[:sample_rate*4], good_cycles[:sample_rate*4, 2])\nplt.title('Good cycles')\nplt.xlabel('Time (seconds)')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Most cycles pass the checks but a few do fail. Closer inspection shows that\nthese cycles tend to have large distortions or have very low amplitudes.\nEither way, the sift has not found a clear oscillation so these cycles\nshould be interpreted with caution.\n\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We can use the information in ``all_cycles`` to find explore the cycle\ncontent of each IMF. For instance, this section prints the number of cycles\nidentified in each IMF\n\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"msg = 'IMF-{0} contains {1:3d} cycles of which {2:3d} ({3}%) are good'\nfor ii in range(all_cycles.shape[1]):\n all_count = all_cycles[:, ii].max()\n good_count = good_cycles[:, ii].max()\n percent = np.round(100*(good_count/all_count), 1)\n\n print(msg.format(ii+1, all_count, good_count, percent))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"IMF-3 contains our simulated oscillation with a spectral peak around 15Hz.\nAs we would expect, the cycle detection finds around 150 cycles in this\n10 second segment. Many of these cycles pass our cycle-quality checks\nindicating that they have well behaved instantaneous phase profiles that can\nbe interpreted in detail. Some cycles do not pass, indicating that parts of\nIMF-3 may not contain a strong oscillatory signal.\n\nThe lower frequency cycles (IMF 4+) have fewer and fewer cycles reflecting\ntheir slowing frequency content (each successive IMF extracts slower dynamics\nthan the previous one). Again, most of these cycles pass the quality checks\non their instantaneous phase.\nHowever, we also see that the the\nhigher frequency IMFs 0 and 1 seem to contain fewer cycles than IMF-3. We\nwould expect these IMFs to capture faster dynamics with more cycles in each\nIMF - so why are there fewer here?\n\nThe answer is in the quality of the instantaneous phase estimation in very\nfast oscillations. Lets plot the IMF and IP for IMF-1\n\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"plt.figure(figsize=(10, 8))\nplt.subplots_adjust(hspace=.3)\nplt.subplot(311)\nplt.plot(t[:sample_rate//2], imf[:sample_rate//2, 0], 'k')\nplt.gca().set_xticklabels([])\nplt.title('IMF-0')\n\nplt.subplot(312)\nplt.plot(t[:sample_rate//2], IP[:sample_rate//2, 0], 'b')\n\nplt.plot((0, .5), (0, 0), label='0')\nplt.plot((0, .5), (np.pi*2, np.pi*2), label='2pi')\nplt.plot((0, .5), (np.pi*2-np.pi/12, np.pi*2-np.pi/12), ':', label='Upper Thresh')\nplt.plot((0, .5), (np.pi/12, np.pi/12), ':', label='Lower Thresh')\n\nplt.gca().set_xticklabels([])\nplt.title('Instantanous Phase')\nplt.ylabel('Radians')\n\nplt.subplot(313)\nplt.plot(t[:sample_rate//2], good_cycles[:sample_rate//2, 0])\nplt.title('Good cycles')\nplt.xlabel('Time (seconds)')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The oscillations here are much faster than in IMF2. We only have a handful of\nsamples for each potential cycle in IMF-1 compared to ~40 for IMF-3. As such,\nmore cycles are showing distortions and failing the quality checks. In this\ncase it is ok as there is no signal in IMF-1 in our simulation. Much of IMF-1\nis noisy for this sift. We could potentially improve this by changing the\nsift parameters to compute more iterations for each IMF. This would increase\nthe number of good cycles in IMF-1 but might lead to over-sifting in other\nIMFs. These parameters should be tuned for the priorities of each analysis.\n\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Further Reading & References\n\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Andrew J. Quinn, V\u00edtor Lopes-dos-Santos, Norden Huang, Wei-Kuang Liang, Chi-Hung Juan, Jia-Rong Yeh, Anna C. Nobre, David Dupret, and Mark W. Woolrich (2001)\nWithin-cycle instantaneous frequency profiles report oscillatory waveform dynamics\nJournal of Neurophysiology 126:4, 1190-1208\nhttps://doi.org/10.1152/jn.00201.2021\n\n"
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.7.9"
}
},
"nbformat": 4,
"nbformat_minor": 0
}