/
slidingS4SigmaPHI.m
executable file
·122 lines (112 loc) · 5.39 KB
/
slidingS4SigmaPHI.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
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
function [S4out, SPHout] = slidingS4SigmaPHI(Powin, PHin, N)
% derived from slidingavg
% OUTPUT_ARRAY = SLIDINGAVG(INPUT_ARRAY, N)
%
% The function 'slidingavg' implements a one-dimensional filtering,
% applying a sliding window to a sequence. Such filtering replaces the
% center value in the window with the average value of all the points
% within the window. When the sliding window is exceeding the lower or
% upper boundaries of the input vector INPUT_ARRAY, the average is
% computed among the available points. Indicating with nx the length of
% the the input sequence, we note that for values of N larger or equal to
% 2*(nx - 1), each value of the output data array are identical and equal
% to mean(in).
%
% * The input argument INPUT_ARRAY is the numerical data array to be
% processed. * The input argument N is the number of neighboring data
% points to average over for each point of IN.
%
% * The output argument OUTPUT_ARRAY is the output data array.
%
% © 2002 - Michele Giugliano, PhD and Maura Arsiero
% (Bern, Friday July 5th, 2002 - 21:10)
% (http://www.giugliano.info) (bug-reports to michele@giugliano.info)
%
% Two simple examples with second- and third-order filters are
% slidingavg([4 3 5 2 8 9 1],2)
% ans =
% 3.5000 4.0000 3.3333 5.0000 6.3333 6.0000 5.0000
%
% slidingavg([4 3 5 2 8 9 1],3)
% ans =
% 3.5000 4.0000 3.3333 5.0000 6.3333 6.0000 5.0000
%
if (isempty(Powin)) || (isempty(PHin)) || (N<=0)
% If the input array is empty or N is non-positive,% an error is reported
% to the standard output and the % execution of the routine is stopped.
disp(sprintf('SlidingAvg: (Error) empty input data or N null.'));
return;
end % if
%find amplitude
amp = sqrt(Powin);
% if (N==1) If the number of neighbouring points over which the sliding
% S4out = Powin; SPHout = PHin;
% out = in;
% average will be performed is '1', then no average actually occur and
% return;
% OUTPUT_ARRAY will be the copy of INPUT_ARRAY and the execution of the
% routine end % if is stopped.
nx = length(amp); %=length(PHin)
% The length of the input data structure is acquired to later evaluate the
% 'mean' over the appropriate boundaries.
% If the number of neighbouring points over which the sliding average
% will be performed is large enough, then the average actually covers
% all the points
% of INPUT_ARRAY, for each index of OUTPUT_ARRAY and some CPU time can be
% gained by such an approach.
% The execution of the routine is stopped.
if (N>=(2*(nx-1)))
S4out= sqrt((mean(Powin) - mean(amp).^2)/ mean(amp).^2)*ones(size(Powin));
SPHout = std(PHin)*ones(size(PHin)); %sigma_phi
return;
end % if
S4out = zeros(size(Powin));
% In all the other situations, the initialization of the output data
% structure is performed.
SPHout = zeros(size(PHin));
% When N is even, then we proceed in taking the half of it:% m = N
% / 2.% Otherwise (N >= 3, N odd), N-1 is even ( N-1 >= 2) and we
% proceed taking the half of it:
% m = (N-1) / 2.
if rem(N,2)~=1
m = N/2;
else
m = (N-1)/2;
end % if
for i=1:nx,
% For each element (i-th) contained in the input numerical array, a
% check must be performed:
if ((i-m) < 1) && ((i+m) <= nx)
% If not enough points are available on the left of the i-th element
% then we proceed to evaluate the mean from the first element to the
% (i + m)-th.
S4out(i)= sqrt((mean(Powin(1:i+m)) - mean(amp(1:i+m)).^2)/ ...
mean(amp(1:i+m)).^2);
SPHout(i) = std(PHin(1:i+m));
% out(i) = mean(in(1:i+m));
elseif ((i-m) >= 1) && ((i+m) <= nx)
% If enough points are available on the left and on the right of the
% i-th element then we proceed to evaluate the mean on 2*m
% elements centered on the i-th position.
S4out(i)= sqrt((mean(Powin(i-m:i+m)) - mean(amp(i-m:i+m)).^2)/ ...
mean(amp(i-m:i+m)).^2);
SPHout(i) = std(PHin(i-m:i+m));
% out(i) = mean(in(i-m:i+m));
elseif ((i-m) >= 1) && ((i+m) > nx)
% If not enough points are available on the rigth of the i-th element
% then we proceed to evaluate the mean from the element (i - m)-th to
% the last one.
S4out(i)= sqrt((mean(Powin(i-m:nx)) - mean(amp(i-m:nx)).^2)/ ...
mean(amp(i-m:nx)).^2);
SPHout(i) = std(PHin(i-m:nx));
% out(i) = mean(in(i-m:nx));
elseif ((i-m) < 1) && ((i+m) > nx)
% If not enough points are available on the left and on the rigth of
% the i-th element then we proceed to evaluate the mean from the first
% element to the last.
S4out(i)= sqrt((mean(Powin(1:nx)) - mean(amp(1:nx)).^2)/ ...
mean(amp(1:nx)).^2);
SPHout(i) = std(PHin(1:nx));
% out(i) = mean(in(1:nx));
end % if
end % for i