pkg load signal;
Imports and Setup
available_graphics_toolkits()
ans = {
[1,1] = fltk
[1,2] = gnuplot
[1,3] = notebook
[1,4] = plotly
}
= "plotly" ; ## Main toolkit
mtk = "gnuplot";
mtk = "gnuplot"; ## Alternative toolkit atk
graphics_toolkit('notebook')
Functions in Octave
function x = square(a)
= a.^2
x endfunction
function x = cube(a)
= a.^3
x endfunction
2); square(
x = 4
= [1 2] a
a | 1 | 2 |
1 | 1 | 2 |
= square(a); s
x | 1 | 2 |
1 | 1 | 4 |
; cube(a)
x | 1 | 2 |
1 | 1 | 8 |
s
s | 1 | 2 |
1 | 1 | 4 |
Signal Creation
Random Signal Creation
= 15
p rand(p,1)
p = 15
ans | 1 |
1 | 0.3882 |
2 | 0.669533 |
3 | 0.391751 |
4 | 0.291379 |
5 | 0.194003 |
6 | 0.717703 |
7 | 0.693793 |
8 | 0.360285 |
9 | 0.444018 |
10 | 0.336876 |
11 | 0.864597 |
12 | 0.493347 |
13 | 0.489672 |
14 | 0.850089 |
15 | 0.217253 |
= 100
n linspace(1, p, n);
n = 100
= 1000; % Hz
srate = 15
p = 5;
noiseamp time = 0:1/srate:3;
= length(time);
n = interp1(rand(p,1)*30, linspace(1,p, n));
ampl = noiseamp*randn(size(time));
noise = ampl+noise; signal
p = 15
function signal = random_signal(srate=1000, p=15, noiseamp=5)
time = 0:1/srate:3;
= length(time);
n = interp1(rand(p,1)*30, linspace(1, p, n));
ampl = noiseamp*randn(size(time));
noise = ampl+noise;
signal endfunction
figure(1), clf
plot(random_signal())
Spike Time Series
= 300
n
= round(exp(randn(n,1))*10);
isi
= 0;
spike_ts for i=1:n
length(spike_ts)+isi(i))=1;
spike_ts(end
figure(8), clf
plot(spike_ts)
n = 300
hist(isi, 30);
set(gcf,'Visible','on')
Filter
Mean Filter
\[ y_t = \frac{1}{2k+1}\sum_{i=t-k}^{t+k}x_i \]
function filtsig = filter_mean(signal, k=20, initzeros=true)
if initzeros
= zeros(size(signal));
filtsig else
= signal;
filtsig end
= length(signal);
n for i=k+1:n-k-1
i) = mean(signal(i-k:i+k));
filtsig(end
endfunction
= random_signal();
signal = filter_mean(signal, k=20, initzeros=false);
filtsig ; filtsig
Compute windowsize in ms
= 1000*(k*2+1)/srate windowsize
windowsize = 41
help dsearchn
'dsearchn' is a function from the file /home/rahuketu/anaconda3/envs/octave/share/octave/7.2.0/m/geometry/dsearchn.m
-- IDX = dsearchn (X, TRI, XI)
-- IDX = dsearchn (X, TRI, XI, OUTVAL)
-- IDX = dsearchn (X, XI)
-- [IDX, D] = dsearchn (...)
Return the index IDX of the closest point in X to the elements XI.
If OUTVAL is supplied, then the values of XI that are not contained
within one of the simplices TRI are set to OUTVAL. Generally, TRI
is returned from ‘delaunayn (X)’.
The optional output D contains a column vector of distances between
the query points XI and the nearest simplex points X.
See also: dsearch, tsearch.
Additional help for built-in functions and operators is
available in the online version of the manual. Use the command
'doc <topic>' to search the manual index.
Help and information about Octave is also available on the WWW
at https://www.octave.org and via the help@octave.org
mailing list.
figure(2), clf, hold on
plot(time, signal, time, filtsig, 'linewidth', 2)
= dsearchn(time',1);
tidx ylim = get(gca,'ylim');
patch(time([ tidx-k tidx-k tidx+k tidx+k ]),ylim([ 1 2 2 1 ]),'k','facealpha',.25,'linestyle','none');
plot(time([tidx tidx]),ylim,'k');
xlabel('Time (sec.)'), ylabel('Amplitude');
title([ 'Running-mean filter with a k=' num2str(round(windowsize)) '-ms filter' ]);
legend({'Signal';'Filtered';'Window';'window center'});
zoom on
= dsearchn(time', 1); tidx tidx
tidx = 1001
help patch
'patch' is a function from the file /home/rahuketu/anaconda3/envs/octave/share/octave/7.2.0/m/plot/draw/patch.m
-- patch ()
-- patch (X, Y, C)
-- patch (X, Y, Z, C)
-- patch ("Faces", FACES, "Vertices", VERTS, ...)
-- patch (..., PROP, VAL, ...)
-- patch (..., PROPSTRUCT, ...)
-- patch (HAX, ...)
-- H = patch (...)
Create patch object in the current axes with vertices at locations
(X, Y) and of color C.
If the vertices are matrices of size MxN then each polygon patch
has M vertices and a total of N polygons will be created. If some
polygons do not have M vertices use NaN to represent "no vertex".
If the Z input is present then 3-D patches will be created.
The color argument C can take many forms. To create polygons which
all share a single color use a string value (e.g., "r" for red), a
scalar value which is scaled by ‘caxis’ and indexed into the
current colormap, or a 3-element RGB vector with the precise
TrueColor.
If C is a vector of length N then the ith polygon will have a color
determined by scaling entry C(i) according to ‘caxis’ and then
indexing into the current colormap. More complicated coloring
situations require directly manipulating patch property/value
pairs.
Instead of specifying polygons by matrices X and Y, it is possible
to present a unique list of vertices and then a list of polygon
faces created from those vertices. In this case the "Vertices"
matrix will be an Nx2 (2-D patch) or Nx3 (3-D patch). The MxN
"Faces" matrix describes M polygons having N vertices—each row
describes a single polygon and each column entry is an index into
the "Vertices" matrix to identify a vertex. The patch object can
be created by directly passing the property/value pairs
"Vertices"/VERTS, "Faces"/FACES as inputs.
Instead of using property/value pairs, any property can be set by
passing a structure PROPSTRUCT with the respective field names.
If the first argument HAX is an axes handle, then plot into this
axes, rather than the current axes returned by ‘gca’.
The optional return value H is a graphics handle to the created
patch object.
Programming Note: The full list of properties is documented at
Patch Properties. Useful patch properties include: "cdata",
"edgecolor", "facecolor", "faces", and "facevertexcdata".
See also: fill, get, set.
Additional help for built-in functions and operators is
available in the online version of the manual. Use the command
'doc <topic>' to search the manual index.
Help and information about Octave is also available on the WWW
at https://www.octave.org and via the help@octave.org
mailing list.
= [ 1 2; 2 1] A
A | 1 | 2 |
1 | 1 | 2 |
2 | 2 | 1 |
figure(3), clf, hold on
= [0 2 4 6; 8 10 12 14; 16 18 20 22];
C imagesc(C)
colorbar
figure(4), clf, hold on
= [0 2 4 6; 8 10 12 14; 16 18 20 22];
C imagesc(C)
colorbar
Uniform PDF
= 20
k ones(1, 2*k+1)/(2*k+1)
k = 20
ans | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 | 9 | 10 | 11 | 12 | 13 | 14 | 15 | 16 | 17 | 18 | 19 | 20 | 21 | 22 | 23 | 24 | 25 | 26 | 27 | 28 | 29 | 30 | 31 | 32 | 33 | 34 | 35 | 36 | 37 | 38 | 39 | 40 | 41 |
1 | 0.0243902 | 0.0243902 | 0.0243902 | 0.0243902 | 0.0243902 | 0.0243902 | 0.0243902 | 0.0243902 | 0.0243902 | 0.0243902 | 0.0243902 | 0.0243902 | 0.0243902 | 0.0243902 | 0.0243902 | 0.0243902 | 0.0243902 | 0.0243902 | 0.0243902 | 0.0243902 | 0.0243902 | 0.0243902 | 0.0243902 | 0.0243902 | 0.0243902 | 0.0243902 | 0.0243902 | 0.0243902 | 0.0243902 | 0.0243902 | 0.0243902 | 0.0243902 | 0.0243902 | 0.0243902 | 0.0243902 | 0.0243902 | 0.0243902 | 0.0243902 | 0.0243902 | 0.0243902 | 0.0243902 |
function g = uniform(k, normalize=false)
= ones(1, 2*k+1);
g if normalize
= g/sum(g);
g end
endfunction
=20) uniform(k
ans | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 | 9 | 10 | 11 | 12 | 13 | 14 | 15 | 16 | 17 | 18 | 19 | 20 | 21 | 22 | 23 | 24 | 25 | 26 | 27 | 28 | 29 | 30 | 31 | 32 | 33 | 34 | 35 | 36 | 37 | 38 | 39 | 40 | 41 |
1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 |
Gaussian Filter
\[ y_t = \sum_{i=t-k}^{t+k}x_ig_i \]
where g is
\[ g = e^{\frac{-4 ln(2) t^2}{w^2}} \]
w = Full width half minimum(width at 0.5 gain)
Notes - Gaussian smoothes out edges even more than mean filter
log(2), exp(2)
= linspace(1, 2, 5)
t = 2
w -4*log(2)*t.^2/ w^2
ans = 0.6931
ans = 7.3891
t | 1 | 2 | 3 | 4 | 5 |
1 | 1 | 1.25 | 1.5 | 1.75 | 2 |
w = 2
ans | 1 | 2 | 3 | 4 | 5 |
1 | -0.693147 | -1.08304 | -1.55958 | -2.12276 | -2.77259 |
Gaussian PDF
function g = gaussian(t, w, normalize=false)
= exp( -4*log(2)*t.^2/w^2);
g if normalize
= g/sum(g)
g end
endfunction
= linspace(-4, 4, 50);
t figure(5), clf
plot(gaussian(t,w=2))
= 1000;
srate = 25.5483213524
fwhm = 50;
k = 1000*(-k:k)/srate;
gtime = gaussian(gtime, fwhm);
gaus_win = k + dsearchn(gaus_win(k+1:end)', 0.5)
pre_peak_half = dsearchn(gaus_win(1:k)', 0.5)
post_peak_half
= gtime(pre_peak_half) - gtime(post_peak_half) exp_fwhm
fwhm = 25.548
pre_peak_half = 64
post_peak_half = 38
exp_fwhm = 26
figure(6), clf , hold on
plot(gtime, gaus_win, 'ko-', 'markerfacecolor', 'w', 'linewidth', 2)
plot(gtime([pre_peak_half post_peak_half]),gaus_win([pre_peak_half post_peak_half]),'m','linewidth',3)
= gaus_win / sum(gaus_win);
gaus_win title([ 'Gaussian kernel with requeted FWHM ' num2str(fwhm) ' ms (' num2str(exp_fwhm) ' ms achieved)' ])
xlabel('Time (ms)'), ylabel('Gain')
Apply Filter Function
function filtsig = apply_filter(signal, filter_pdf, k=20, initzeros=true)
if initzeros
= zeros(size(signal));
filtsig else
= signal;
filtsig end
= length(signal);
n for i=k+1:n-k-1
i) = sum(signal(i-k:i+k).*filter_pdf);
filtsig(end
endfunction
Mean vs Gaussian Filter
= 1000;
srate = 25.5483213524
fwhm = 40;
k = 1000*(-k:k)/srate;
gtime = gaussian(gtime, fwhm, normalize=true);
gaus_win = uniform(k, normalize=true);
uniform_win = apply_filter(signal, uniform_win, k=k, initzeros=false);
filtsig_uniform = apply_filter(signal, gaus_win, k=k, initzeros=false); filtsig_gaussian
fwhm = 25.548
g | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 | 9 | 10 | 11 | 12 | 13 | 14 | 15 | 16 | 17 | 18 | 19 | 20 | 21 | 22 | 23 | 24 | 25 | 26 | 27 | 28 | 29 | 30 | 31 | 32 | 33 | 34 | 35 | 36 | 37 | 38 | 39 | 40 | 41 | 42 | 43 | 44 | 45 | 46 | 47 | 48 | 49 | 50 | 51 | 52 | 53 | 54 | 55 | 56 | 57 | 58 | 59 | 60 | 61 | 62 | 63 | 64 | 65 | 66 | 67 | 68 | 69 | 70 | 71 | 72 | 73 | 74 | 75 | 76 | 77 | 78 | 79 | 80 | 81 |
1 | 4.11089e-05 | 5.75008e-05 | 7.97484e-05 | 0.000109668 | 0.000149537 | 0.000202176 | 0.000271031 | 0.000360263 | 0.000474822 | 0.000620514 | 0.000804051 | 0.00103306 | 0.00131607 | 0.00166242 | 0.00208216 | 0.00258582 | 0.00318414 | 0.00388774 | 0.00470666 | 0.00564987 | 0.00672472 | 0.00793635 | 0.00928705 | 0.0107757 | 0.0123972 | 0.014142 | 0.0159959 | 0.0179399 | 0.0199498 | 0.0219973 | 0.0240497 | 0.0260711 | 0.0280234 | 0.0298671 | 0.0315628 | 0.0330725 | 0.0343614 | 0.0353984 | 0.0361583 | 0.036622 | 0.0367779 | 0.036622 | 0.0361583 | 0.0353984 | 0.0343614 | 0.0330725 | 0.0315628 | 0.0298671 | 0.0280234 | 0.0260711 | 0.0240497 | 0.0219973 | 0.0199498 | 0.0179399 | 0.0159959 | 0.014142 | 0.0123972 | 0.0107757 | 0.00928705 | 0.00793635 | 0.00672472 | 0.00564987 | 0.00470666 | 0.00388774 | 0.00318414 | 0.00258582 | 0.00208216 | 0.00166242 | 0.00131607 | 0.00103306 | 0.000804051 | 0.000620514 | 0.000474822 | 0.000360263 | 0.000271031 | 0.000202176 | 0.000149537 | 0.000109668 | 7.97484e-05 | 5.75008e-05 | 4.11089e-05 |
figure(7), clf, hold on
plot(time, signal, time, filtsig_uniform,'m','linewidth', 2, time, filtsig_gaussian,'g', 'linewidth', 2)
legend('signal', 'mean', 'gaussian')
Gaussian on Spike Time series
= 25
fwhm = 100;
k = -k:k;
gtime = gaussian(gtime, fwhm, normalize=true);
gaus_win2 = apply_filter(spike_ts, gaus_win2, k=k, initzeros=true);
filtsig_spike figure(12)
plot(gaus_win2)
fwhm = 25
g | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 | 9 | 10 | 11 | 12 | 13 | 14 | 15 | 16 | 17 | 18 | 19 | 20 | 21 | 22 | 23 | 24 | 25 | 26 | 27 | 28 | 29 | 30 | 31 | 32 | 33 | 34 | 35 | 36 | 37 | 38 | 39 | 40 | 41 | 42 | 43 | 44 | 45 | 46 | 47 | 48 | 49 | 50 | 51 | 52 | 53 | 54 | 55 | 56 | 57 | 58 | 59 | 60 | 61 | 62 | 63 | 64 | 65 | 66 | 67 | 68 | 69 | 70 | 71 | 72 | 73 | 74 | 75 | 76 | 77 | 78 | 79 | 80 | 81 | 82 | 83 | 84 | 85 | 86 | 87 | 88 | 89 | 90 | 91 | 92 | 93 | 94 | 95 | 96 | 97 | 98 | 99 | 100 | 101 | 102 | 103 | 104 | 105 | 106 | 107 | 108 | 109 | 110 | 111 | 112 | 113 | 114 | 115 | 116 | 117 | 118 | 119 | 120 | 121 | 122 | 123 | 124 | 125 | 126 | 127 | 128 | 129 | 130 | 131 | 132 | 133 | 134 | 135 | 136 | 137 | 138 | 139 | 140 | 141 | 142 | 143 | 144 | 145 | 146 | 147 | 148 | 149 | 150 | 151 | 152 | 153 | 154 | 155 | 156 | 157 | 158 | 159 | 160 | 161 | 162 | 163 | 164 | 165 | 166 | 167 | 168 | 169 | 170 | 171 | 172 | 173 | 174 | 175 | 176 | 177 | 178 | 179 | 180 | 181 | 182 | 183 | 184 | 185 | 186 | 187 | 188 | 189 | 190 | 191 | 192 | 193 | 194 | 195 | 196 | 197 | 198 | 199 | 200 | 201 |
1 | 2.03708e-21 | 4.92493e-21 | 1.18015e-20 | 2.80301e-20 | 6.59867e-20 | 1.5397e-19 | 3.56091e-19 | 8.1627e-19 | 1.85461e-18 | 4.17657e-18 | 9.32251e-18 | 2.06249e-17 | 4.52272e-17 | 9.82999e-17 | 2.11765e-16 | 4.52169e-16 | 9.56962e-16 | 2.00741e-15 | 4.17372e-15 | 8.60118e-15 | 1.75687e-14 | 3.55686e-14 | 7.13744e-14 | 1.41959e-13 | 2.79855e-13 | 5.46824e-13 | 1.05903e-12 | 2.03291e-12 | 3.8679e-12 | 7.2942e-12 | 1.36341e-11 | 2.52594e-11 | 4.63838e-11 | 8.44222e-11 | 1.52298e-10 | 2.72318e-10 | 4.82622e-10 | 8.47782e-10 | 1.47607e-09 | 2.54729e-09 | 4.35709e-09 | 7.38688e-09 | 1.24129e-08 | 2.06743e-08 | 3.41299e-08 | 5.58454e-08 | 9.05703e-08 | 1.4559e-07 | 2.31965e-07 | 3.66321e-07 | 5.73387e-07 | 8.89571e-07 | 1.36792e-06 | 2.0849e-06 | 3.14963e-06 | 4.71606e-06 | 6.99916e-06 | 1.02958e-05 | 1.50113e-05 | 2.16934e-05 | 3.10728e-05 | 4.41145e-05 | 6.20768e-05 | 8.65812e-05 | 0.000119692 | 0.000164003 | 0.000222734 | 0.000299826 | 0.000400035 | 0.000529021 | 0.000693418 | 0.000900874 | 0.00116006 | 0.00148062 | 0.00187306 | 0.00234859 | 0.00291884 | 0.00359551 | 0.00438993 | 0.00531252 | 0.00637222 | 0.00757579 | 0.00892713 | 0.0104266 | 0.0120704 | 0.0138498 | 0.0157513 | 0.0177555 | 0.019838 | 0.021969 | 0.024114 | 0.0262346 | 0.0282895 | 0.030236 | 0.032031 | 0.0336328 | 0.0350028 | 0.0361068 | 0.0369166 | 0.0374112 | 0.0375775 | 0.0374112 | 0.0369166 | 0.0361068 | 0.0350028 | 0.0336328 | 0.032031 | 0.030236 | 0.0282895 | 0.0262346 | 0.024114 | 0.021969 | 0.019838 | 0.0177555 | 0.0157513 | 0.0138498 | 0.0120704 | 0.0104266 | 0.00892713 | 0.00757579 | 0.00637222 | 0.00531252 | 0.00438993 | 0.00359551 | 0.00291884 | 0.00234859 | 0.00187306 | 0.00148062 | 0.00116006 | 0.000900874 | 0.000693418 | 0.000529021 | 0.000400035 | 0.000299826 | 0.000222734 | 0.000164003 | 0.000119692 | 8.65812e-05 | 6.20768e-05 | 4.41145e-05 | 3.10728e-05 | 2.16934e-05 | 1.50113e-05 | 1.02958e-05 | 6.99916e-06 | 4.71606e-06 | 3.14963e-06 | 2.0849e-06 | 1.36792e-06 | 8.89571e-07 | 5.73387e-07 | 3.66321e-07 | 2.31965e-07 | 1.4559e-07 | 9.05703e-08 | 5.58454e-08 | 3.41299e-08 | 2.06743e-08 | 1.24129e-08 | 7.38688e-09 | 4.35709e-09 | 2.54729e-09 | 1.47607e-09 | 8.47782e-10 | 4.82622e-10 | 2.72318e-10 | 1.52298e-10 | 8.44222e-11 | 4.63838e-11 | 2.52594e-11 | 1.36341e-11 | 7.2942e-12 | 3.8679e-12 | 2.03291e-12 | 1.05903e-12 | 5.46824e-13 | 2.79855e-13 | 1.41959e-13 | 7.13744e-14 | 3.55686e-14 | 1.75687e-14 | 8.60118e-15 | 4.17372e-15 | 2.00741e-15 | 9.56962e-16 | 4.52169e-16 | 2.11765e-16 | 9.82999e-17 | 4.52272e-17 | 2.06249e-17 | 9.32251e-18 | 4.17657e-18 | 1.85461e-18 | 8.1627e-19 | 3.56091e-19 | 1.5397e-19 | 6.59867e-20 | 2.80301e-20 | 1.18015e-20 | 4.92493e-21 | 2.03708e-21 |
figure(13), clf, hold on
plot(spike_ts)
plot(filtsig_spike, 'linewidth', 2)
Note
- Filtered Signal/ Red line can be interpreted as probability of observing the spike at any given time point.
- It is sensitive to k and fwhm.
Denoising EMG Signal with TKEO
- When a movement is initiated / muscle twitch there is a spike in the signal
- Computer it might be difficult to differentiate from background noise to spike of signal and background noise
Teager Kaiser energy-tracking operator
\[ y_t = x_t^2 - x_{t-1}x_{t+1} \]
- Substracts background noise and make a signal even larger
load(fullfile(cd, "resources","section_02","emg4TKEO.mat" )); whos
Variables visible from the current scope:
variables in scope: top scope
Attr Name Size Bytes Class
==== ==== ==== ===== =====
A 2x2 32 double
C 3x4 96 double
a 1x2 16 double
ampl 1x3001 24008 double
ans 1x1 8 double
atk 1x7 7 char
doc_file 1x79 79 char
emg 1x1281 5124 single
emgtime 1x1281 10248 double
exp_fwhm 1x1 8 double
filtsig 1x3001 24008 double
filtsig_gaussian 1x3001 24008 double
filtsig_spike 1x5415 43320 double
filtsig_uniform 1x3001 24008 double
fs 1x1 8 double
fwhm 1x1 8 double
gaus_win 1x81 648 double
gaus_win2 1x201 1608 double
gtime 1x201 24 double
i 1x1 8 double
initzeros 1x1 1 logical
isi 300x1 2400 double
k 1x1 8 double
mtk 1x7 7 char
n 1x1 8 double
noise 1x3001 24008 double
noiseamp 1x1 8 double
normalize 1x1 1 logical
p 1x1 8 double
pkg_dir 1x64 64 char
post_peak_half 1x1 8 double
pre_peak_half 1x1 8 double
s 1x2 16 double
signal 1x3001 24008 double
spike_ts 1x5415 43320 double
srate 1x1 8 double
t 1x50 400 double
tidx 1x1 8 double
time 1x3001 24 double
uniform_win 1x81 648 double
w 1x1 8 double
windowsize 1x1 8 double
ylim 1x2 16 double
Total is 35509 elements using 252267 bytes
figure(14), clf, hold on
plot(emgtime, emg)
function tkeo = apply_tkeo(signal)
= signal;
tkeo for i=2:length(signal)-1
i) = signal(i)^2 - signal(i-1)*signal(i+1);
tkeo(end
endfunction
function tkeo = apply_tkeo_vec(signal)
= signal;
tkeo 2:end-1) = signal(2:end-1).^2 - signal(1: end-2).*signal(3:end);
tkeo(endfunction
figure(15), clf, hold on
plot(emgtime, emg)
plot(emgtime, apply_tkeo(emg))
function signalz = to_zscore(signal, signaltime, time_cutoff=0)
= dsearchn(signaltime', time_cutoff);
time0 = (signal - mean(signaltime(1:time0)))/ std(signaltime(1:time0));
signalz endfunction
figure(16), clf
= apply_tkeo(emg);
tkeo subplot(211), hold on
plot(emgtime, emg./max(emg), 'b', 'linewidth', 2)
plot(emgtime, tkeo./max(tkeo), 'm', 'linewidth', 2)
xlabel('Time (ms)'), ylabel('Amplitude or energy')
legend("EMG", "EMG energy (TKEO)")
subplot(212), hold on
plot(emgtime, to_zscore(emg, emgtime, time_cutoff=0), 'b', 'linewidth', 2)
plot(emgtime, to_zscore(tkeo, emgtime, time_cutoff=0), 'm', 'linewidth', 2)
xlabel('Time (ms)'), ylabel('Amplitude or energy')
legend("EMG", "EMG energy (TKEO)")
Note
- time_cutoff defines baseline activity cutoff. ( Remaining signal contains’ more interesting artifacts)
- tkeo signal and original signal are in different units.So we need to apply zscore to make them comparable
- Background noise above in original signal becomes a flatline in tkeo signal
- tkeo in zscore plot very easy to detect a change ( 800 times standard deviation)
Median filtering
- Median to remove spikes
- Non linear filter
- Apply to selected data points, not all data points
, xx] = hist(signal, 30) [nn
nn | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 | 9 | 10 | 11 | 12 | 13 | 14 | 15 | 16 | 17 | 18 | 19 | 20 | 21 | 22 | 23 | 24 | 25 | 26 | 27 | 28 | 29 | 30 |
1 | 2 | 1 | 2 | 5 | 7 | 27 | 42 | 58 | 84 | 117 | 147 | 161 | 191 | 220 | 196 | 193 | 170 | 178 | 179 | 168 | 162 | 158 | 130 | 123 | 110 | 59 | 61 | 22 | 16 | 12 |
xx | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 | 9 | 10 | 11 | 12 | 13 | 14 | 15 | 16 | 17 | 18 | 19 | 20 | 21 | 22 | 23 | 24 | 25 | 26 | 27 | 28 | 29 | 30 |
1 | -14.3157 | -12.5474 | -10.779 | -9.01073 | -7.24241 | -5.47409 | -3.70577 | -1.93745 | -0.169126 | 1.59919 | 3.36751 | 5.13584 | 6.90416 | 8.67248 | 10.4408 | 12.2091 | 13.9774 | 15.7458 | 17.5141 | 19.2824 | 21.0507 | 22.819 | 24.5874 | 26.3557 | 28.124 | 29.8923 | 31.6606 | 33.429 | 35.1973 | 36.9656 |
hist(signal, 30)