1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
|
// Function to implement spectral factorization, as discussed in sec. 13.1.
// 13.2
function [r,b,rbbr] = spec1(A,dA,B,dB,rho)
AA = rho * convol(A,flip(A));
BB = convol(B,flip(B));
diff1 = dA - dB;
dBB = 2*dB;
for i = 1:diff1
[BB,dBB] = polmul(BB,dBB,[0 1],1);
end
[rbbr,drbbr] = poladd(AA,2*dA,BB,dBB);
rts = roots(rbbr); // roots in descending order of magnitude
rts = flip(rts);
rtsin = rts(dA+1:2*dA);
b = 1;
for i = 1:dA,
b = convol(b,[1 -rtsin(i)]);
end
br = flip(b);
bbr = convol(b,br);
r = rbbr(1) / bbr(1);
endfunction;
|
