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
|
function [val_ar t_ar f_val_ar EC] = AR_estimate( t, val, freq, nu, mu, expand_factor)
% [val_ar t_ar f_val_ar EC] = AR_estimate( t, val, freq, < nu, mu, expand_factor >)
%
% apply autoregressive signal model to improve dft results
%
% t : time vector
% val : time domain values
% freq : frequency vector for dft
%
% optional
% nu : AR order (default 40)
% mu : number of timesteps to train the model (default 3*nu)
% expand_factor : increase signal length by this factor (default 5)
%
% return values:
% val_ar: AR estimated time signal
% t_ar: time vector
% f_val_ar: FD transformed AR estimated signal
% EC: error code
% 0 --> no error
% 1 --> input error: t and val mismatch
% 2 --> input error: mu has to be larger than 2*nu
% 3 --> inout error: expand_factor has to be larger than 1
% 10 --> AR error: signal is to short for AR estimate --> decrease AR order
% 11 --> AR error: estimated signal appears to be unstable --> use a different mu
%
% openEMS matlab interface
% -----------------------
% Author: Thorsten Liebig, 2011
%
% See also ReadUI, DFT_time2freq
EC = 0;
val_ar = [];
t_ar = [];
f_val_ar = [];
if numel(t) ~= numel(val)
if (nargout<4)
error 'numel(t) ~= numel(val)'
else
EC = 1;
return
end
end
if (nargin<4)
nu = 40;
end
if (nargin<5)
mu = 3*nu;
end
if (nargin<6)
expand_factor=5;
end
if (mu<=2*nu)
if (nargout<4)
error 'mu has to be larger than 2*nu'
else
EC = 2;
return
end
end
if (expand_factor<=1)
if (nargout<4)
error 'expand_factor has to be larger than 1'
else
EC = 3;
return
end
end
dt = t(2)-t(1);
M = numel(t);
if (M<0.6*mu)
if (nargout<4)
error 'signal is to short for AR estimate --> decrease AR order'
else
EC = 10;
return
end
end
for n=1:mu-nu
b(n) = val(end-n+1);
for m=1:nu
A(n,m)=val(end-n+1-m);
end
end
a = ((A'*A)\A')*b';
val_ar = val;
t_ar = t;
for k=M:expand_factor*M
val_ar(k) = 0;
t_ar(k) = t_ar(k-1)+dt;
val_ar(k) = sum(a.*val_ar(k-(1:nu))');
end
if (max(val_ar(M:end)) > max(val))
if (nargout<4)
error 'estimated signal appears to be unstable --> use a different mu'
else
EC = 11;
return
end
end
f_val_ar = DFT_time2freq(t_ar, val_ar, freq);
|
