{
"cells": [
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"%matplotlib inline"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n# Why use EMD?\nThe motivations behind the Empirical Mode Demposition are straightforward but\nsometimes get lost in a relatively complicated and technical literature. This\ntutorial presents a quick summary of what EMD could add to your analysis and\nwhat issues it can help to solve.\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Frequency spectra of dynamic and non-sinusoidal signals\nA core challenge in signal processing is finding an intuitive representation\nof the frequency content of complicated and dynamic signals. The most common\napproach is to use methods based on the Fourier-transform - in which we\nrepresent a signal with a set of sinusoidal basis functions. This is a very\npowerful and flexible approach but has some short-comings when working with\ndynamic or non-sinusoidal signals.\n\nLet's take a look at an example. Here we make two oscillatory signals - one\nwith a sinusoidal waveform and one with a non-linearity.\n\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"import emd\nimport numpy as np\nfrom scipy import signal, ndimage\nimport matplotlib.pyplot as plt\n\nsample_rate = 1000\nseconds = 10\nnum_samples = sample_rate*seconds\n\ntime_vect = np.linspace(0, seconds, num_samples)\n\nfreq = 5\n\n# Sinusoidal signal\nx = np.cos(2*np.pi*freq*time_vect)\n\n# Non-linear signal\ny = np.cos(2*np.pi*freq*time_vect) + 0.25*np.cos(2*np.pi*freq*2*time_vect-np.pi)\n\n# Quick summary figure\nplt.figure(figsize=(8, 4))\nplt.plot(time_vect[:sample_rate], x[:sample_rate])\nplt.plot(time_vect[:sample_rate], y[:sample_rate])\nplt.legend(['Sinsuoidal', 'Non-linear'])\n\n# sphinx_gallery_thumbnail_number = 2"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We can see that our non-linear signal has a wider trough and a sharper peak\nthan the sinusoid - many natural signals might have oscillatory features like\nthis. Unfortunately, the Fourier transform has to use multiple components to\nrepresent this oscillation as its basis set fixed and unchanging over time.\n\nAs an illustration, lets compute a fast-Fourier transform on our two signals\n\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"pxx = np.fft.fft(x) / len(x)\npyy = np.fft.fft(y) / len(x)\nfft_f = np.fft.fftfreq(x.shape[0], d=time_vect[1]-time_vect[0])\n\nplt.figure()\nplt.plot(fft_f, np.abs(pxx))\nplt.plot(fft_f, np.abs(pyy), '--')\nplt.xlim(0, 20)\nplt.legend(['Sinusoidal', 'Non-linear'])\nplt.xlabel('Frequency (Hz)')\nplt.ylabel('Power')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Welch's periodogram shows us that both signals have high power at 5Hz - as we\nwould expect from our signal definition. However, the non-linear case has an\nadditional peak at 10Hz. This arises as the Fourier transform and related\nmethods needs to represent the non-linearity with a combination of strictly\nlinear basis functions. We only need a signal component to represent our\nlinear signal but the non-linear signal needs to use two. This is well\nillustrated in the time domain.\n\nFirst, we can use the parameters of the Fourier transform to recompute our basis signals.\n\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"comp1_amp = 2 * np.abs(pyy[50])\ncomp1_phase = np.arctan2(pyy[50].imag, pyy[50].real)\ncomp1 = comp1_amp * np.cos(2*np.pi*5*time_vect + comp1_phase)\n\ncomp2_amp = 2 * np.abs(pyy[100])\ncomp2_phase = np.arctan2(pyy[100].imag, pyy[100].real)\ncomp2 = comp2_amp * np.cos(2*np.pi*10*time_vect + comp2_phase)\n\nprint('Component 1 amp={0} and phase={1}'.format(comp1_amp, comp1_phase))\nprint('Component 2 amp={0} and phase={1}'.format(comp2_amp, comp2_phase))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The values recovered from the fourier transform are very close to the\nsimulated values from the start. We see that the first component has an\namplitude very close to 1 and a phase close to zero, whilst the second\ncomponent has an amplitude close to 0.25 and a phase close to -pi.\n\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"plt.figure(figsize=(8, 4), frameon=False)\nplt.yticks([])\nplt.plot(time_vect[:sample_rate], y[:sample_rate] + 3)\nplt.plot(time_vect[:sample_rate], comp1[:sample_rate])\nplt.plot(time_vect[:sample_rate], comp2[:sample_rate])\nplt.plot(time_vect[:sample_rate], comp1[:sample_rate] + comp2[:sample_rate] - 3)\nplt.legend(['Original signal', 'Basis component 1', 'Basis component 2', 'Basis component 1 + 2'])"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The top signal in blue is our simulated signal, the orange and green signals\nare the two components identified by the Fourier transform coefficients which\nsum together to create the bottom signal in red.\n\nIt turns out that any signal can be represented by the summation of sets of\nsimple sinusoids - even extremely non-linear or disjoint signals like\ntriangular waves. This is a very powerful tool, but the need to include\nmultiple high frequency components, known as harmonics, can create\ninterpretation issues. In a complex or noisy signal, it might not be clear\nwhether a peak in a power spectrum represents a harmonic of a low frequency\nsignal or a stand-alone high frequency signal. The problem gets worse if we\nhave several oscillations which are all creating separate harmonics. These\nsignals may be untangled using more complicated analysis of the interactions\nbetween different Fourier components but these can get very complicated.\n\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Instantaneous Frequency\n\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The Empirical Mode Decomposition offers a different perspective on this\nproblem. Instead of representing dynamic signals with combinations of static\nbasis functions, the EMD looks to isolate a small number of temporally\nadaptive basis functions and derive dynamics in frequency and amplitude\ndirectly from them.\n\nThese adaptive basis functions are called Intrinsic Mode Functions (IMFs) and\nare isolated from the data using the 'sift' algorithm. This is a time-domain\nalgorithm which doesn't require a Fourier transform to separate out different\nsignals.\n\nIMFs are defined by several features. They must be locally-symmetric around\nzero and contain the same number of extrema as zero-crossings (or differ by\nno more than 1). The 'sift' is a numerical algorithm for isolating signal\ncomponents with exactly these features. Any residual after the sift is a\nnon-oscillatory trend component.\n\nLets run a sift on our non-linear signal and plot the results.\n\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"imf = emd.sift.sift(y)\n\nemd.plotting.plot_imfs(imf[:sample_rate, :])"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The sift has identified two components, the first contains the oscillatory\nsignal with its non-sinusoidal shape intact whilst the second contains a\nconstant mean-term.\n\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We can compute the instantaneous amplitude, frequency and phase of our IMF\nusing the Hilbert transform. This is a function which creates an 'analytic'\nform of a signal that is extremely useful for extracting dynamic information\nfrom a signal. Critically, the result of a Hilbert transform is only valid\nfor signals with very specific features. The definition of the sift ensures\nthat each resulting IMF does meet these criteria (or at least be very close).\n\nThis means that we can use the Hilbert transform to estimate instantaneous\nfrequencies directly from our IMFs. We don't need to use the Fourier\ntransform for this and so can avoid adding harmonic components into our\nresults.\n\nHere, we compute the Hilbert transform from our IMFs\n\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"IP, IF, IA = emd.spectra.frequency_transform(imf, sample_rate, 'hilbert')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"...and plot up the resulting instantaneous frequency.\n\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"plt.figure(figsize=(10, 6))\nplt.subplots_adjust(hspace=0.4)\nplt.subplot(211)\nplt.plot(time_vect[:sample_rate], imf[:sample_rate, 0])\nplt.title('IMF-1')\nplt.subplot(212)\nplt.plot(time_vect[:sample_rate], IF[:sample_rate, 0])\nplt.title('IMF-1 Instantaneous Frequency')\nplt.ylabel('Frequency (Hz)')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"This dynamic representation of frequency is able to adapt to the changing\nwaveform shape and describe the non-sinusoidal shape as a single smoothly\nvarying component rather than as two separate, static components.\n\nWe can compute a power spectrum using these principles called the Hilbert-Huang transform\n\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"freq_range = (0.1, 20, 32)\nhht_f, spec = emd.spectra.hilberthuang(IF, IA, freq_range, scaling='density')\n\nplt.figure()\nplt.plot(fft_f, np.abs(pyy))\nplt.plot(hht_f, spec)\nplt.legend(['Fourier Transform', 'Hilbert-Huang Transform'], frameon=False)\nplt.xlim(0, 20)\nplt.xlabel('Frequency (Hz)')\nplt.ylabel('Power')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The Hilbert-Huang transform gives a summary of the dynamic spectral content\nof the signal. It represents our non-linear signal as a distribution of\nfrequencies around 5Hz without adding an additional signal component.\n\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Time-resolution\n\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Finally, the instantaneous frequency and amplitude from the Hilbert transform\nallow us to see very rapid dynamics in oscillatory signals. Fourier methods\nare able to see dynamic changes by applying the transform in a sliding window\nbut the choice of window size places an important constraint on the\ntime-frequency resolution of the transform. Longer windows have greater\nfrequency resolution as they are able to use a larger number of Fourier basis\ncomponents - however these longer windows necessarily lead to worse time\nresolution.\n\nIn contrast, the instantaneous frequency stats naturally vary over time, so\nno window choice is required. Let's take a closer look with a dynamic\nsimulation.\n\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"# Create a signal with a dynamic oscillation\nz = emd.simulate.ar_oscillator(25, sample_rate, seconds, r=0.975)[:, 0]"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"First we compute the sift, frequency transform and the Hilbert-Huang spectrum\n\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"imf = emd.sift.sift(z)\nIP, IF, IA = emd.spectra.frequency_transform(imf, sample_rate, 'hilbert')\n\nfreq_range = (0.1, 50, 48)\nhht_f, hht = emd.spectra.hilberthuang(IF, IA, freq_range, mode='amplitude', sum_time=False)\nhht = ndimage.gaussian_filter(hht, 1)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Next, we compute a short-time Fourier-transform. The ``nperseg`` variable\ncontrols the window-length and time-frequency resolution. Higher values will\nlead to greater precision in frequency and worse precision in time.\n\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"nperseg = 2048\nftf, ftt, ftZ = signal.stft(z, nperseg=nperseg, fs=sample_rate, noverlap=nperseg-1)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"plt.figure(figsize=(8, 10))\nplt.subplot(311)\nplt.plot(time_vect, z)\nplt.xlim(1, 9)\nplt.title('Signal')\nplt.subplot(312)\nplt.pcolormesh(ftt, ftf, np.abs(ftZ), cmap='hot_r')\nplt.ylim(0, 50)\nplt.xlim(1, 9)\nplt.title('Short-Time Fourier-Transform')\nplt.ylabel('Frequency (Hz)')\nplt.subplot(313)\nplt.pcolormesh(time_vect, hht_f, hht, cmap='hot_r')\nplt.ylim(0, 50)\nplt.title('Hilbert-Huang Transform')\nplt.xlim(1, 9)\nplt.xlabel('Time (seconds)')\nplt.ylabel('Frequency (Hz)')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Changes in oscillatory power are visible in both transforms but are noticeably\nsmoother in the STFT compared to the HHT. The blur in the STFT is due to the\nwindowing that we have applied, we can try to reduce this window size to\nincrease the temporal resolution of our transform\n\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"nperseg = 512\nftf, ftt, ftZ = signal.stft(z, nperseg=nperseg, fs=sample_rate, noverlap=nperseg-1)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Now we plot the two transforms together for another comparison.\n\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"plt.figure(figsize=(8, 10))\nplt.subplot(311)\nplt.plot(time_vect, z)\nplt.xlim(1, 9)\nplt.title('Signal')\nplt.subplot(312)\nplt.pcolormesh(ftt, ftf, np.abs(ftZ), cmap='hot_r')\nplt.ylim(0, 50)\nplt.xlim(1, 9)\nplt.title('Short-Time Fourier-Transform')\nplt.ylabel('Frequency (Hz)')\nplt.subplot(313)\nplt.pcolormesh(time_vect, hht_f, hht, cmap='hot_r')\nplt.ylim(0, 50)\nplt.title('Hilbert-Huang Transform')\nplt.xlim(1, 9)\nplt.xlabel('Time (seconds)')\nplt.ylabel('Frequency (Hz)')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Shortening the window-size in the STFT has reduced the temporal blur but\nincreased the blur in the frequency dimension. The HHT produces a sharper\nimage in both dimensions.\n\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The remaining tutorials go into all of the above analysis in much more detail\nand try to introduce the motivations and concepts behind EMD alongside a\nrange of practical examples. We also introduce some problems with the EMD and\nhow to recognise & reduce them - there is no free-lunch in signal processing!\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
}