forked from orcioni/Volterra2.0
-
Notifications
You must be signed in to change notification settings - Fork 0
/
LeeSch5.m
executable file
·79 lines (76 loc) · 3.1 KB
/
LeeSch5.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
%function [k5]=LeeSch5(xn,yn,os,R,A,swap,delay)
%
% k5 is the fifth order kernel of Wiener series according to Lee-Schetzen
% method, which will contain not-a-NaN values ony within the set of
% fundamental points, which is the minimum set of non-diagonal points which
% allows the reconstruction of non-diagonal kernel points by symmetry
% (see symmetrize function).
% For diagonal points refer to WVdiag5 function.
%
% xn is the input sequence.
%
% yn the output sequence.
%
% os is the input/output sequences index from where the cross-correlation is
% started, all the sequence values before os are thrown. In can be used when
% xn and yn have been obtained from an A/D conversion and we the initial
% transient conditions cut away.
%
% R is the length corresponding to length(k5)-1. The lags domain interval
% corresponding to k5 is [0,R]x[0,R]x[0,R]x[0,R]x[0,R].
%
% A is the second order moment of xn (i.e. power).
%
% swap it's an optional parameter for speed optimization (see fastxcorr for
% insights), it depends on the machine you're on, and on having compiled this
% code or not.
% If you don't know what to do use default value.
% (We choose 23 in octave environment and 11 in the standalone version)
%
% delay gives the result restricted to the lags domain interval
% [0+delay,R]x[0+delay,R]x[0+delay,R]x[0+delay,R]x[0+delay,R],
% most useful for higher order kernels.
%
% If you want to contact the authors, please write to [email protected],
% or Simone Orcioni, DII, Università Politecnica delle Marche,
% via Brecce Bianche, 12 - 60131 Ancona, Italy.
% Copyright (C) 2006 Massimiliano Pirani
%
% This program is free software; you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation; either version 2 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License along
% with this program; if not, write to the Free Software Foundation, Inc.,
% 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
function [k5]=LeeSch5(xn,yn,os,R,A,swap,delay)
L=length(xn(os:end-delay));
A5=A*A*A*A*A;
k5=NaNmat(R+1,R+1,R+1,R+1,R+1);
% sgm2=3:R-1
% sgm1=sgm2+1:R
% sgm3=2:sgm2-1
% sgm4=1:sgm3-1
% sgm5=0:sgm4-1
for sgm2=3:R-1
p2=[zeros(sgm2,1);xn(os:end-sgm2-delay)];
for sgm1=sgm2+1:R;
p1=[zeros(sgm1,1);xn(os:end-sgm1-delay)];
for sgm3=2:sgm2-1;
p3=[zeros(sgm3,1);xn(os:end-sgm3-delay)];
for sgm4=1:sgm3-1
k5(sgm1+1,sgm2+1,sgm3+1,sgm4+1,1:sgm4)=frxcorr(xn(os:end-delay),...
[zeros(sgm4,1);xn(os:end-sgm4-delay)].*p3.*p2.*p1.*yn(os+delay:end),...
L,...
sgm4-1,0,swap)/120/A5;
end
end
end
end
end