1function [l1, m1, l2, m2] = src_factor2_lm(fs1, fs2)
2
3% factor2_lm - factorize input and output sample rates ratio to fraction l1/m2*l2/m2
4%
5% [l1, m1, l2, m2] = factor2_lm(fs1, fs2)
6%
7% fs1 - input sample rate
8% fs2 - output sample rate
9%
10
11% Copyright (c) 2016, Intel Corporation
12% All rights reserved.
13%
14% Redistribution and use in source and binary forms, with or without
15% modification, are permitted provided that the following conditions are met:
16%   * Redistributions of source code must retain the above copyright
17%     notice, this list of conditions and the following disclaimer.
18%   * Redistributions in binary form must reproduce the above copyright
19%     notice, this list of conditions and the following disclaimer in the
20%     documentation and/or other materials provided with the distribution.
21%   * Neither the name of the Intel Corporation nor the
22%     names of its contributors may be used to endorse or promote products
23%     derived from this software without specific prior written permission.
24%
25% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
26% AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
27% IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
28% ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
29% LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
30% CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
31% SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
32% INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
33% CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
34% ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
35% POSSIBILITY OF SUCH DAMAGE.
36%
37% Author: Seppo Ingalsuo <seppo.ingalsuo@linux.intel.com>
38%
39
40
41%% The same sample rate?
42if abs(fs1-fs2) < eps
43        l1 = 1; m1 = 1; l2 = 1; m2 = 1;
44        return
45end
46
47%% Find fraction, and factorize
48k = gcd(fs1, fs2);
49l = fs2/k;
50m = fs1/k;
51[l01, l02] = factor2(l);
52[m01, m02] = factor2(m);
53
54%% Hand fixing for some ratios, guide to reuse some common filters
55if (l==147) && (m==640)
56        l01 = 7; m01 = 8; l02 = l/l01; m02 = m/m01; % 192 to 44.1
57end
58if (l==147) && (m==320)
59        l01 = 7; m01 = 8; l02 = l/l01; m02 = m/m01; % 96 to 44.1
60end
61if (l==147) && (m==160)
62        l01 = 7; m01 = 8; l02 = l/l01; m02 = m/m01; % 48 to 44.1
63end
64if (l==160) && (m==147)
65        l01 = 8; m01 = 7; l02 = l/l01; m02 = m/m01; % 44.1 to 48
66end
67if (l==320) && (m==147)
68        l01 = 8; m01 = 7; l02 = l/l01; m02 = m/m01; % 44.1 to 96
69end
70if (l==4) && (m==3)
71        l01 = 4; m01 = 3; l02 = l/l01; m02 = m/m01; % 24 to 32, no 2 stage
72end
73if (l==3) && (m==4)
74        l01 = 3; m01 = 4; l02 = l/l01; m02 = m/m01; % 24 to 32, no 2 stage
75end
76
77
78r11 = l01/m01;
79r22 = l02/m02;
80r12 = l01/m02;
81r21 = l02/m01;
82fs3 = [r11 r12 r21 r22].*fs1;
83
84
85if fs1 > fs2 % Decrease sample rate, dont go below output rate
86        idx = find(fs3 >= fs2);
87        if length(idx) < 1
88                error('Cant factorise interpolations');
89        end
90        fs3_possible = fs3(idx);
91        delta = fs3_possible-fs2; % Fs3 that is nearest to fs2
92        idx = find(delta == min(delta), 1, 'first');
93        fs3_min = fs3_possible(idx);
94        idx = find(fs3 == fs3_min, 1, 'first');
95        switch idx
96                case 1, l1 = l01; m1 = m01; l2 = l02; m2 = m02;
97                case 2, l1 = l01; m1 = m02; l2 = l02; m2 = m01;
98                case 3, l1 = l02; m1 = m01; l2 = l01; m2 = m02;
99                case 4, l1 = l02; m1 = m02; l2 = l01; m2 = m01;
100        end
101else % Increase sample rate, don't go below input rate
102        idx = find(fs3 >= fs1);
103        if length(idx) < 1
104                error('Cant factorise interpolations');
105        end
106        fs3_possible = fs3(idx);
107        delta = fs3_possible-fs1; % Fs2 that is nearest to fs1
108        idx = find(delta == min(delta), 1, 'first');
109        fs3_min = fs3_possible(idx);
110        idx = find(fs3 == fs3_min, 1, 'first');
111        switch idx
112                case 1, l1 = l01; m1 = m01; l2 = l02; m2 = m02;
113                case 2, l1 = l01; m1 = m02; l2 = l02; m2 = m01;
114                case 3, l1 = l02; m1 = m01; l2 = l01; m2 = m02;
115                case 4, l1 = l02; m1 = m02; l2 = l01; m2 = m01;
116        end
117end
118
119%% If 1st stage is 1:1
120if (l1 == 1) && (m1 == 1)
121        l1 = l2;
122        m1 = m2;
123        l2 = 1;
124        m2 = 1;
125end
126
127f1=l/m;
128f2=l1/m1*l2/m2;
129if abs(f1 - f2) > 1e-6
130        error('Bug in factorization code!');
131end
132
133end
134
135
136function [a, b]=factor2(c)
137x = round(sqrt(c));
138t1 = x:2*x;
139d1 = c./t1;
140idx1 = find(d1 == floor(d1), 1, 'first');
141a1 = t1(idx1);
142t2 = x:-1:floor(x/2);
143d2 = c./t2;
144idx2 = find(d2 == floor(d2), 1, 'first');
145a2 = t2(idx2);
146if (a1-x) < (x-a2)
147        a = a1;
148else
149        a = a2;
150end
151b = c/a;
152end
153
154
155
156