-
Notifications
You must be signed in to change notification settings - Fork 2
/
absmean_norm.m
56 lines (48 loc) · 1.37 KB
/
absmean_norm.m
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
function [ data ] = absmean_norm( data, delta )
%UNTITLED Summary of this function goes here
% Detailed explanation goes here
para_initial('absmean');
global FREQ_LOW_ABSM
global FREQ_HIGH_ABSM
global FILTER_ORDER
global WATER_LEVEL
global WINLEN_ABSM
global HAMPEL_WIN
global MAX_ITER
freqlow = FREQ_LOW_ABSM;
freqhigh = FREQ_HIGH_ABSM;
filter_order = FILTER_ORDER;
water_level = WATER_LEVEL;
winlen = WINLEN_ABSM;
hampel_win = HAMPEL_WIN;
max_iter = MAX_ITER;
para_initial('clear');
hampel_npts = floor(hampel_win / delta);
if freqhigh >= 0.5 / delta
ME = MException(...
'arguSetting:exceedingNyquist', ...
'Corner freq. %f greater than Nyquist %f!', ...
freqhigh, 0.5 / delta);
throw(ME)
end
[b, a] = butter(filter_order, [freqlow, freqhigh] .* delta .* 2);
while max_iter
data = data / max(abs(data));
dtemp = filter(b, a, data);
weight = smooth(abs(dtemp), winlen);
data = data ./ weight;
data(isnan(data)) = 0;
data = hampel(data, hampel_npts);
trun_bool = abs(data) > (water_level * rms(data));
num_above = sum(trun_bool);
trun_bool = trun_bool + ...
(trun_bool == 0) .* 1 + ...
(trun_bool == 1) .* (1 / water_level - 1);
data = data .* trun_bool;
if num_above == 0
break
else
max_iter = max_iter - 1;
end
end
end