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
|
//Compute linear model using Steiglitz-McBride iteration
//calling syntax
//[b,a] = stmcb(h,nb,na)
//[b,a] = stmcb(y,x,nb,na)
//[b,a] = stmcb(h,nb,na,niter)
//[b,a] = stmcb(y,x,nb,na,niter)
//[b,a] = stmcb(h,nb,na,niter,ai)
//[b,a] = stmcb(y,x,nb,na,niter,ai)
//Parameters :
//b,a : coefficients of the system function,nb is number of zeros and na is number of poles
//h:impulse response of the system
//x,y: input and output of same length given to the system
//niter: no of iterations
//ai:initial estimate of the denominator coefficients
//Accepts only real i/ps , imaginary i/ps are not accepted due to limitations of the 'filter' function in Scilab
//Example
//h = fscanfMat("macros/stmcb_h_data.txt");
//stmcb(h,4,4)
//Output :
// ans =
//
// 0.0003 0.0010284 0.0147159 - 0.0077914 0.0316548
function [b,a] = stmcb( x, u_in, q, p, niter, a_in )
narginchk(3, 6, argn(2));
//modify stmcb to handle exceptions when i/p is char
if(type(x)==10 | type(u_in)==10) then
error("Input in stmcb must be double/single, instead it was char");
end
if length(u_in) == 1 then
if argn(2) == 3 then
niter = 5;
p = q;
q = u_in;
a_in = prony(x, 0, p);
elseif argn(2) == 4
niter = p;
p = q;
q = u_in;
a_in = prony(x, 0, p);
elseif argn(2) == 5
a_in = niter;
niter = p;
p = q;
q = u_in;
end
[x_row, x_col] = size(x);
u_in = zeros(x_row, x_col);
u_in(1) = 1;
else
if length(u_in) ~= length(x) then
error('Input Signal x and Output Signal Y must have the same length');
end
if argn(2) < 6 then
[b, a_in] = prony(x, 0, p);
end
if argn(2) < 5 then
niter = 5;
end
end
a = a_in;
N = length(x);
for i = 1:niter
u = filter(1, a, x);
v = filter(1, a, u_in);
C1 = convmtx(u(:),p+1);
C2 = convmtx(v(:),q+1);
T = [ -C1(1:N,:) C2(1:N,:)];
c = T(:,2:p+q+2)\(-T(:,1));
a = [1; c(1:p)];
b = c(p+1:p+q+1);
end
a = a.';
b = b.';
endfunction
function narginchk(min_argin, max_argin, num_of_argin)
if num_of_argin < min_argin then
error('Not enough input arguments')
end
if num_of_argin > max_argin then
error('Too many input arguments')
end
endfunction
|
