tune/src/src_factor2_lm.m

function [l1, m1, l2, m2] = src_factor2_lm(fs1, fs2)

% factor2_lm - factorize input and output sample rates ratio to fraction l1/m2*l2/m2
%
% [l1, m1, l2, m2] = factor2_lm(fs1, fs2)
%
% fs1 - input sample rate
% fs2 - output sample rate
%

% Copyright (c) 2016, Intel Corporation
% All rights reserved.
%
% Redistribution and use in source and binary forms, with or without
% modification, are permitted provided that the following conditions are met:
%   * Redistributions of source code must retain the above copyright
%     notice, this list of conditions and the following disclaimer.
%   * Redistributions in binary form must reproduce the above copyright
%     notice, this list of conditions and the following disclaimer in the
%     documentation and/or other materials provided with the distribution.
%   * Neither the name of the Intel Corporation nor the
%     names of its contributors may be used to endorse or promote products
%     derived from this software without specific prior written permission.
%
% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
% AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
% IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
% ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
% LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
% CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
% SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
% INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
% CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
% ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
% POSSIBILITY OF SUCH DAMAGE.
%
% Author: Seppo Ingalsuo <seppo.ingalsuo@linux.intel.com>
%


%% The same sample rate?
if abs(fs1-fs2) < eps
        l1 = 1; m1 = 1; l2 = 1; m2 = 1;
        return
end

%% Find fraction, and factorize
k = gcd(fs1, fs2);
l = fs2/k;
m = fs1/k;
[l01, l02] = factor2(l);
[m01, m02] = factor2(m);

%% Hand fixing for some ratios, guide to reuse some common filters
if (l==147) && (m==640)
        l01 = 7; m01 = 8; l02 = l/l01; m02 = m/m01; % 192 to 44.1
end
if (l==147) && (m==320)
        l01 = 7; m01 = 8; l02 = l/l01; m02 = m/m01; % 96 to 44.1
end
if (l==147) && (m==160)
        l01 = 7; m01 = 8; l02 = l/l01; m02 = m/m01; % 48 to 44.1
end
if (l==160) && (m==147)
        l01 = 8; m01 = 7; l02 = l/l01; m02 = m/m01; % 44.1 to 48
end
if (l==320) && (m==147)
        l01 = 8; m01 = 7; l02 = l/l01; m02 = m/m01; % 44.1 to 96
end
if (l==4) && (m==3)
        l01 = 4; m01 = 3; l02 = l/l01; m02 = m/m01; % 24 to 32, no 2 stage
end
if (l==3) && (m==4)
        l01 = 3; m01 = 4; l02 = l/l01; m02 = m/m01; % 24 to 32, no 2 stage
end


r11 = l01/m01;
r22 = l02/m02;
r12 = l01/m02;
r21 = l02/m01;
fs3 = [r11 r12 r21 r22].*fs1;


if fs1 > fs2 % Decrease sample rate, dont go below output rate
        idx = find(fs3 >= fs2);
        if length(idx) < 1
                error('Cant factorise interpolations');
        end
        fs3_possible = fs3(idx);
        delta = fs3_possible-fs2; % Fs3 that is nearest to fs2
        idx = find(delta == min(delta), 1, 'first');
        fs3_min = fs3_possible(idx);
        idx = find(fs3 == fs3_min, 1, 'first');
        switch idx
                case 1, l1 = l01; m1 = m01; l2 = l02; m2 = m02;
                case 2, l1 = l01; m1 = m02; l2 = l02; m2 = m01;
                case 3, l1 = l02; m1 = m01; l2 = l01; m2 = m02;
                case 4, l1 = l02; m1 = m02; l2 = l01; m2 = m01;
        end
else % Increase sample rate, don't go below input rate
        idx = find(fs3 >= fs1);
        if length(idx) < 1
                error('Cant factorise interpolations');
        end
        fs3_possible = fs3(idx);
        delta = fs3_possible-fs1; % Fs2 that is nearest to fs1
        idx = find(delta == min(delta), 1, 'first');
        fs3_min = fs3_possible(idx);
        idx = find(fs3 == fs3_min, 1, 'first');
        switch idx
                case 1, l1 = l01; m1 = m01; l2 = l02; m2 = m02;
                case 2, l1 = l01; m1 = m02; l2 = l02; m2 = m01;
                case 3, l1 = l02; m1 = m01; l2 = l01; m2 = m02;
                case 4, l1 = l02; m1 = m02; l2 = l01; m2 = m01;
        end
end

%% If 1st stage is 1:1
if (l1 == 1) && (m1 == 1)
        l1 = l2;
        m1 = m2;
        l2 = 1;
        m2 = 1;
end

f1=l/m;
f2=l1/m1*l2/m2;
if abs(f1 - f2) > 1e-6
        error('Bug in factorization code!');
end

end


function [a, b]=factor2(c)
x = round(sqrt(c));
t1 = x:2*x;
d1 = c./t1;
idx1 = find(d1 == floor(d1), 1, 'first');
a1 = t1(idx1);
t2 = x:-1:floor(x/2);
d2 = c./t2;
idx2 = find(d2 == floor(d2), 1, 'first');
a2 = t2(idx2);
if (a1-x) < (x-a2)
        a = a1;
else
        a = a2;
end
b = c/a;
end