summaryrefslogtreecommitdiff
path: root/macros/poly2lsf.sci
blob: 69288f35ef78d6420f881cae30a5c79a4942512b (plain)
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
function lsf=poly2lsf(a)


//poly2lsf Prediction polynomial to line spectral frequencies.

// Calling Sequence
// lsf = poly2lsf(a)
// Parameters
// k: define the prediction polynomial.
// lsf: returns corresponding line spectral frequencies.
// Examples
//X = [0.5 0.3 0.8 0.9 0.4 0.05];
// lsf = poly2lsf(X)
// See also
//
// Author
// Jitendra Singh
//modified to match MATLAb o/p
if(find(a(1,1)==0)==1) then
    error ("Input to ROOTS must not contain NaN or Inf");
end //end  of modification


if or(type(a)==10) then
    error ('The polynomial must have all roots inside of the unit circle.')
end



if (size(a,1) > 1) & (size(a,2) > 1)
    error('The prediction polynomial must be stored in a vector.')
end


if ~isreal(a) then
    error('Line spectral frequencies are not defined for complex polynomials.')
end

if (max(abs(roots(a))) >= 1) then
    error('The polynomial must have all roots inside of the unit circle.');
end

if a(1) ~= 1 then
    a = a./a(1);
end

a=a(:);
p  = length(a)-1;
a1 = [a;0];
a2 = a1($:-1:1);
pp = a1-a2;
qq = a1+a2;

if (p-fix(p./2).*2)~=0 then  // Odd order
    //////////////////////////////
    aa=[1 0 -1];
    [m,n] = size(pp);
    n = max(m,n);
    na = length(aa);
    if na > n
        P = 0;
        
    else
        P= filter(pp, aa, [1 zeros(1,n-na)]);
        if m ~= 1
            P = P(:);
        end
    end
    
    Q = qq;
else          // Even order
    
    aa=[1 -1];
    [m,n] = size(pp);
    n = max(m,n);
    na = length
    (aa);
    if na > n
        P = 0;
        
    else
        P= filter(pp, aa, [1 zeros(1,n-na)]);
        if m ~= 1 
            P = P(:);
        end
    end
    
    aa=[1 1];
    [m,n] = size(qq);
    n = max(m,n);
    na = length(aa);
    
    if na > n
        Q = 0;
        
    else
        Q= filter(qq, aa, [1 zeros(1,n-na)]);
        if m ~= 1
            Q = Q(:);
        end
    end
    
end

r_p  = roots(P); r_q  = roots(Q);

ap  =atan(imag(r_p(1:2:$)),real(r_p(1:2:$)));
aq  = atan(imag(r_q(1:2:$)),real(r_q(1:2:$)));

lsf = mtlb_sort([ap;aq]);
endfunction