This patch adds the scripts src_std_int32.m and src_tiny_int16.m those were used to create the conversions matrices for the default and tiny profiles. The scripts call the included src_generate function that does the task with help from the filter design and coefficients export utilities.

The quality and resources consumption of SRC can be customized with this tool.

Signed-off-by: Seppo Ingalsuo seppo.ingalsuo@linux.intel.com --- tune/src/README | 36 +++++ tune/src/src_export_coef.m | 212 ++++++++++++++++++++++++++++ tune/src/src_export_defines.m | 66 +++++++++ tune/src/src_export_table_2s.m | 190 +++++++++++++++++++++++++ tune/src/src_factor1_lm.m | 44 ++++++ tune/src/src_factor2_lm.m | 155 ++++++++++++++++++++ tune/src/src_find_l0m0.m | 71 ++++++++++ tune/src/src_generate.m | 311 +++++++++++++++++++++++++++++++++++++++++ tune/src/src_get.m | 254 +++++++++++++++++++++++++++++++++ tune/src/src_param.m | 87 ++++++++++++ tune/src/src_std_int32.m | 30 ++++ tune/src/src_tiny_int16.m | 41 ++++++ 12 files changed, 1497 insertions(+) create mode 100644 tune/src/README create mode 100644 tune/src/src_export_coef.m create mode 100644 tune/src/src_export_defines.m create mode 100644 tune/src/src_export_table_2s.m create mode 100644 tune/src/src_factor1_lm.m create mode 100644 tune/src/src_factor2_lm.m create mode 100644 tune/src/src_find_l0m0.m create mode 100644 tune/src/src_generate.m create mode 100644 tune/src/src_get.m create mode 100644 tune/src/src_param.m create mode 100644 tune/src/src_std_int32.m create mode 100644 tune/src/src_tiny_int16.m

diff --git a/tune/src/README b/tune/src/README new file mode 100644 index 0000000..e135f9a --- /dev/null +++ b/tune/src/README @@ -0,0 +1,36 @@ +Sample rate converter (SRC) Setup Tools +======================================= + +This is a tool to set up SRC conversions sample rates list, define +quality related parameters, and test the C implementation for a number +of objective audio quality parameters. + +The tools need GNU Octave version 4.0.0 or later with octave-signal +package. + +src_std_int32.m +--------------- + +This script creates the default coefficient set and contains nothing +else but call for src_generate. + +src_tiny_int16.m +---------------- + +This script creates the tiny coefficient set. The script contains an +example how to customize the input/output rates matrix and in a simple +way the scale conversions quality. More controlled quality adjust can +be done by editing file src_param.m directly. Note that int16 +presentation for SRC coefficients will degrade even the default +quality. + +src_generate.m +-------------- + +Creates the header files to include to C into directory "include". A +report of create modes is written to directory "reports". The +coefficients need to be copied to directory +sof.git/src/include/sof/audio/coefficients/src. + +The default quality of SRC is defined in module src_param.m. The +quality impacts the complexity and coefficents tables size of SRC. diff --git a/tune/src/src_export_coef.m b/tune/src/src_export_coef.m new file mode 100644 index 0000000..f6fdbb4 --- /dev/null +++ b/tune/src/src_export_coef.m @@ -0,0 +1,212 @@ +function success=src_export_coef(src, ctype, vtype, hdir, profile) + +% src_export_coef - Export FIR coefficients +% +% success=src_export_coef(src, ctype, hdir, profile) +% +% src - src definition struct +% ctype - 'float','int32', or 'int24' +% vtype - 'float','int32_t' +% hdir - directory for header files +% profile - string to append to filename +% + +% 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 +% + +if nargin < 5 + profile = ''; +end + +if src.L == src.M + success = 0; +else + pbi = round(src.c_pb*1e4); + sbi = round(src.c_sb*1e4); + if isempty(profile) + hfn = sprintf('%s/src_%s_%d_%d_%d_%d.h', ... + hdir, ctype, src.L, src.M, pbi, sbi); + else + hfn = sprintf('%s/src_%s_%s_%d_%d_%d_%d.h', ... + hdir, profile, ctype, src.L, src.M, pbi, sbi); + end + vfn = sprintf('src_%s_%d_%d_%d_%d_fir', ctype, src.L, src.M, pbi, sbi); + sfn = sprintf('src_%s_%d_%d_%d_%d', ctype, src.L, src.M, pbi, sbi); + + fprintf('Exporting %s ...\n', hfn); + fh = fopen(hfn, 'w'); + + switch ctype + case 'float' + fprintf(fh, 'const %s %s[%d] = {\n', ... + vtype, vfn, src.filter_length); + fprintf(fh,'\t%16.9e', src.coefs(1)); + for n=2:src.filter_length + fprintf(fh, ',\n'); + fprintf(fh,'\t%16.9e', src.coefs(n)); + end + fprintf(fh,'\n\n};'); + case 'int32' + print_int_coef(src, fh, vtype, vfn, 32); + case 'int24' + print_int_coef(src, fh, vtype, vfn, 24); + case 'int16' + print_int_coef(src, fh, vtype, vfn, 16); + otherwise + error('Unknown type %s !!!', ctype); + end + + fprintf(fh, '\n'); + switch ctype + case 'float' + fprintf(fh, 'struct src_stage %s = {\n', sfn); + fprintf(fh, '\t%d, %d, %d, %d, %d, %d, %d, %d, %f,\n\t%s};\n', ... + src.idm, src.odm, src.num_of_subfilters, ... + src.subfilter_length, src.filter_length, ... + src.blk_in, src.blk_out, src.halfband, ... + src.gain, vfn); + case { 'int16' 'int24' 'int32' } + fprintf(fh, 'struct src_stage %s = {\n', sfn); + fprintf(fh, '\t%d, %d, %d, %d, %d, %d, %d, %d, %d,\n\t%s};\n', ... + src.idm, src.odm, src.num_of_subfilters, ... + src.subfilter_length, src.filter_length, ... + src.blk_in, src.blk_out, src.halfband, ... + src.shift, vfn); + otherwise + error('Unknown type %s !!!', ctype); + end + %fprintf(fh, '\n'); + fclose(fh); + success = 1; +end + +end + +function print_int_coef(src, fh, vtype, vfn, nbits) + fprintf(fh, 'const %s %s[%d] = {\n', ... + vtype, vfn, src.filter_length); + + cint = coef_quant(src, nbits); + fprintf(fh,'\t%d', cint(1)); + for n=2:src.filter_length + fprintf(fh, ',\n'); + fprintf(fh,'\t%d', cint(n)); + end + fprintf(fh,'\n\n};'); +end + +function cint = coef_quant(src, nbits) + + sref = 2^(nbits-1); + pmax = sref-1; + nmin = -sref; + + if nbits > 16 + %% Round() is OK + cint0 = round(sref*src.coefs); + else + %% Prepare to optimize coefficient quantization + fs = max(src.fs1, src.fs2); + f = linspace(0, fs/2, 1000); + + %% Test sensitivity for stopband and find the most sensitive + % coefficients + sbf = linspace(src.f_sb,fs/2, 500); + n = src.filter_length; + psens = zeros(1,n); + bq0 = round(sref*src.coefs); + h = freqz(bq0/sref/src.L, 1, sbf, fs); + sb1 = 20*log10(sqrt(sum(h.*conj(h)))); + for i=1:n + bq = src.coefs; + bq(i) = round(sref*bq(i))/sref; + %tbq = bq; %tbq(i) = bq(i)+1; + h = freqz(bq, 1, sbf, fs); + psens(i) = sum(h.*conj(h)); + end + [spsens, pidx] = sort(psens, 'descend'); + + %% Test quantization in the found order + % The impact to passband is minimal so it's not tested + bi = round(sref*src.coefs); + bi0 = bi; + dl = -1:1; + nd = length(dl); + msb = zeros(1,nd); + for i=pidx + bit = bi; + for j=1:nd + bit(i) = bi(i) + dl(j); + h = freqz(bit, 1, sbf, fs); + msb(j) = sum(h.*conj(h)); + end + idx = find(msb == min(msb), 1, 'first'); + bi(i) = bi(i) + dl(idx); + end + h = freqz(bi/sref/src.L, 1, sbf, fs); + sb2 = 20*log10(sqrt(sum(h.*conj(h)))); + + %% Plot to compare + if 0 + f = linspace(0, fs/2, 1000); + h1 = freqz(src.coefs/src.L, 1, f, fs); + h2 = freqz(bq0/sref/src.L, 1, f, fs); + h3 = freqz(bi/sref/src.L, 1, f, fs); + figure; + plot(f, 20*log10(abs(h1)), f, 20*log10(abs(h2)), f, 20*log10(abs(h3))); + grid on; + fprintf('Original = %4.1f dB, optimized = %4.1f dB, delta = %4.1f dB\n', ... + sb1, sb2, sb1-sb2); + end + cint0 = bi; + end + + + %% Re-order coefficients for filter implementation + cint = zeros(src.filter_length,1); + for n = 1:src.num_of_subfilters + i11 = (n-1)*src.subfilter_length+1; + i12 = i11+src.subfilter_length-1; + cint(i11:i12) = cint0(n:src.num_of_subfilters:end); + end + + %% Done check for no overflow + max_fix = max(cint); + min_fix = min(cint); + if (max_fix > pmax) + printf('Fixed point coefficient %d exceeded %d\n.', max_fix, pmax); + error('Something went wrong!'); + end + if (min_fix < nmin) + printf('Fixed point coefficient %d exceeded %d\n.', min_fix, nmax); + error('Something went wrong!'); + end + + +end diff --git a/tune/src/src_export_defines.m b/tune/src/src_export_defines.m new file mode 100644 index 0000000..63fb4c2 --- /dev/null +++ b/tune/src/src_export_defines.m @@ -0,0 +1,66 @@ +function src_export_defines(defs, ctype, hdir, profile) + +% src_export_defines - Exports the constants to header files +% +% src_export_defines(defs, ctype, hdir) +% +% defs - defs struct +% ctype - e.g. 'int24' appended to file name +% hdir - directory for header file +% profile - string to append to file name +% + + +% 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 +% + +if nargin < 4 + profile = ''; +end + +if isempty(profile) + hfn = sprintf('src_%s_define.h', ctype); +else + hfn = sprintf('src_%s_%s_define.h', profile, ctype); +end +fh = fopen(fullfile(hdir,hfn), 'w'); +fprintf(fh, '/* SRC constants */\n'); +fprintf(fh, '#define MAX_FIR_DELAY_SIZE %d\n', defs.fir_delay_size); +fprintf(fh, '#define MAX_OUT_DELAY_SIZE %d\n', defs.out_delay_size); +fprintf(fh, '#define MAX_BLK_IN %d\n', defs.blk_in); +fprintf(fh, '#define MAX_BLK_OUT %d\n', defs.blk_out); +fprintf(fh, '#define NUM_IN_FS %d\n', defs.num_in_fs); +fprintf(fh, '#define NUM_OUT_FS %d\n', defs.num_out_fs); +fprintf(fh, '#define STAGE1_TIMES_MAX %d\n', defs.stage1_times_max); +fprintf(fh, '#define STAGE2_TIMES_MAX %d\n', defs.stage2_times_max); +fprintf(fh, '#define STAGE_BUF_SIZE %d\n', defs.stage_buf_size); +fprintf(fh, '#define NUM_ALL_COEFFICIENTS %d\n', defs.sum_filter_lengths); +fclose(fh); + +end diff --git a/tune/src/src_export_table_2s.m b/tune/src/src_export_table_2s.m new file mode 100644 index 0000000..7772b93 --- /dev/null +++ b/tune/src/src_export_table_2s.m @@ -0,0 +1,190 @@ +function sfl = src_export_table_2s(fs_in, fs_out, l_2s, m_2s, ... + pb_2s, sb_2s, taps_2s, ctype, vtype, ppath, hdir, profile) + +% src_export_table_2s - Export src setup table +% +% src_export_table_2s(fs_in, fs_out, l, m, pb, sb, ctype, vtype, hdir, profile) +% +% The parameters are used to differentiate files for possibly many same +% conversion factor filters with possibly different characteristic. +% +% fs_in - input sample rates +% fs_out - output sample rates +% l - interpolation factors +% m - decimation factors +% pb - passband widths +% sb - stopband start frequencies +% ctype - coefficient quantization +% vtype - C variable type +% ppath - print directory prefix to header file name include +% hdir - directory for header files +% profile - string to append to file name +% + +% 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 +% + +if nargin < 12 + profile = ''; +end + +if isempty(profile) + hfn = sprintf('src_%s_table.h', ctype); +else + hfn = sprintf('src_%s_%s_table.h', profile, ctype); +end +fh = fopen(fullfile(hdir,hfn), 'w'); + +fprintf(fh, '/* SRC conversions */\n'); +sfl = 0; +n_in = length(fs_in); +n_out = length(fs_out); +i=1; +all_modes = zeros(2*n_in*n_out, 7); +for n=1:2 + for b=1:n_out + for a=1:n_in + all_modes(i,:) = [ l_2s(n,a,b) m_2s(n,a,b) ... + pb_2s(n,a,b) sb_2s(n,a,b) n a b ]; + i=i+1; + end + end +end + +all_modes_sub = all_modes(:,1:4); +[unique_modes, ia] = unique(all_modes_sub,'rows'); +sm = size(unique_modes); + +if isempty(profile) + prof_ctype = ctype; +else + prof_ctype = sprintf('%s_%s', profile, ctype); +end +for i=1:sm(1) + um_tmp = unique_modes(i,:); + if isequal(um_tmp(1:2),[1 1]) || isequal(um_tmp(1:2),[0 0]) + else + fprintf(fh, '#include <%ssrc_%s_%d_%d_%d_%d.h>\n', ... + ppath, prof_ctype, um_tmp(1:4)); + + n = all_modes(ia(i), 5); + a = all_modes(ia(i), 6); + b = all_modes(ia(i), 7); + sfl = sfl +taps_2s(n, a, b); % Count sum of filter lengths + end +end +fprintf(fh,'\n'); + +fprintf(fh, '/* SRC table */\n'); +switch ctype + case 'float' + fprintf(fh, '%s fir_one = 1.0;\n', vtype); + fprintf(fh, 'struct src_stage src_double_1_1_0_0 = { 0, 0, 1, 1, 1, 1, 1, 0, 1.0, &fir_one };\n'); + fprintf(fh, 'struct src_stage src_double_0_0_0_0 = { 0, 0, 0, 0, 0, 0, 0, 0, 0.0, &fir_one };\n'); + case 'int16' + scale16 = 2^15; + fprintf(fh, '%s fir_one = %d;\n', vtype, round(scale16*0.5)); + fprintf(fh, 'struct src_stage src_int16_1_1_0_0 = { 0, 0, 1, 1, 1, 1, 1, 0, -1, &fir_one };\n'); + fprintf(fh, 'struct src_stage src_int16_0_0_0_0 = { 0, 0, 0, 0, 0, 0, 0, 0, 0, &fir_one };\n'); + case 'int24' + scale24 = 2^23; + fprintf(fh, '%s fir_one = %d;\n', vtype, round(scale24*0.5)); + fprintf(fh, 'struct src_stage src_int24_1_1_0_0 = { 0, 0, 1, 1, 1, 1, 1, 0, -1, &fir_one };\n'); + fprintf(fh, 'struct src_stage src_int24_0_0_0_0 = { 0, 0, 0, 0, 0, 0, 0, 0, 0, &fir_one };\n'); + case 'int32' + scale32 = 2^31; + fprintf(fh, '%s fir_one = %d;\n', vtype, round(scale32*0.5)); + fprintf(fh, 'struct src_stage src_int32_1_1_0_0 = { 0, 0, 1, 1, 1, 1, 1, 0, -1, &fir_one };\n'); + fprintf(fh, 'struct src_stage src_int32_0_0_0_0 = { 0, 0, 0, 0, 0, 0, 0, 0, 0, &fir_one };\n'); + otherwise + error('Unknown coefficient type!'); +end + +fprintf(fh, 'int src_in_fs[%d] = {', n_in); +j = 1; +for i=1:n_in + fprintf(fh, ' %d', fs_in(i)); + if i < n_in + fprintf(fh, ','); + end + j = j + 1; + if (j > 8) + fprintf(fh, '\n\t'); + j = 1; + end +end +fprintf(fh, '};\n'); + +fprintf(fh, 'int src_out_fs[%d] = {', n_out); +j = 1; +for i=1:n_out + fprintf(fh, ' %d', fs_out(i)); + if i < n_out + fprintf(fh, ','); + end + j = j + 1; + if (j > 8) + fprintf(fh, '\n\t'); + j = 1; + end +end +fprintf(fh, '};\n'); + +for n = 1:2 + fprintf(fh, 'struct src_stage *src_table%d[%d][%d] = {\n', ... + n, n_out, n_in); + i = 1; + for b = 1:n_out + fprintf(fh, '\t{'); + for a = 1:n_in + fprintf(fh, ' &src_%s_%d_%d_%d_%d', ... + ctype, l_2s(n,a,b), m_2s(n,a,b), ... + pb_2s(n,a,b), sb_2s(n,a,b)); + if a < n_in + fprintf(fh, ','); + end + i = i + 1; + if i > 2 + fprintf(fh, '\n\t'); + i = 1; + end + end + fprintf(fh, '}'); + if b < n_out + fprintf(fh, ',\n'); + else + fprintf(fh, '\n'); + end + end + fprintf(fh, '};\n'); +end + +fclose(fh); + +end diff --git a/tune/src/src_factor1_lm.m b/tune/src/src_factor1_lm.m new file mode 100644 index 0000000..54e9173 --- /dev/null +++ b/tune/src/src_factor1_lm.m @@ -0,0 +1,44 @@ +function [l, m] = src_factor1_lm(fs1, fs2) + +% factor1_lm - factorize input and output sample rates ratio to fraction l/m +% +% [l, m] = factor1_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 +% + +k = gcd(fs1, fs2); +l = fs2/k; +m = fs1/k; + +end diff --git a/tune/src/src_factor2_lm.m b/tune/src/src_factor2_lm.m new file mode 100644 index 0000000..6ac8602 --- /dev/null +++ b/tune/src/src_factor2_lm.m @@ -0,0 +1,155 @@ +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 + + + diff --git a/tune/src/src_find_l0m0.m b/tune/src/src_find_l0m0.m new file mode 100644 index 0000000..2a360a3 --- /dev/null +++ b/tune/src/src_find_l0m0.m @@ -0,0 +1,71 @@ +function [l0, m0] = src_find_l0m0(L, M) + +% find_l0m0 - find l0, m0 to meet -l0*L + m0*M == 1 +% +% [l0, m0] = find_l0m0(L, M) +% +% L - interpolatation factor +% M - decimation factor +% + +% 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 +% + +if M == 1 + l0 = 0; + m0 = 1; + return +end + +if L == 1 + l0 = 1; + m0 = 0; + return +end + +l0 = []; +m0 = []; +for lt=1:4*L + mt = (1+lt*L)/M; % Check if -lt*L + mt*M == 1 + if floor(mt) == mt + l0 = [l0 lt]; + m0 = [m0 mt]; + end +end + +s = l0+m0; +idx = find(s == min(s), 1, 'first'); +l0 = l0(idx); +m0 = m0(idx); + +if -l0*L + m0*M ~= 1 + error('Something went wrong!'); +end + +end diff --git a/tune/src/src_generate.m b/tune/src/src_generate.m new file mode 100644 index 0000000..b01b29d --- /dev/null +++ b/tune/src/src_generate.m @@ -0,0 +1,311 @@ +function src_generate(fs_in, fs_out, ctype, fs_inout, profile, qc) + +% src_generate - export src conversions for given fs_in and fs_out +% +% src_generate(fs_in, fs_out <, ctype, fs_inout>>) +% +% fs_in - vector of input sample rates (M) +% fs_out - vector of output sample rates (N) +% ctype - coefficient type, use 'int16','int24', 'int32', or 'float' +% fs_inout - matrix of supported conversions (MxN), +% 0 = no support, 1 = supported +% % +% if ctype is omitted then ctype defaults to 'int16'. +% +% If fs_inout matrix is omitted this script will compute coefficients +% for all fs_in <-> fs_out combinations. +% +% 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 +% + +if (nargin == 1) || (nargin > 6) + error('Incorrect arguments for function!'); +end +if nargin == 0 + %% Default input and output rates + fs_in = [8e3 11025 12e3 16e3 18900 22050 24e3 32e3 44100 48e3 ... + 64e3 88.2e3 96e3 176400 192e3]; + + fs_out = [8e3 11025 12e3 16e3 18900 22050 24e3 32e3 44100 48e3]; + + fs_inout = [ 0 0 0 1 0 0 1 1 0 1 ; ... + 0 0 0 0 0 0 0 0 0 1 ; ... + 0 0 0 0 0 0 0 0 0 1 ; ... + 1 0 0 0 0 0 1 1 0 1 ; ... + 0 0 0 0 0 0 0 0 0 1 ; ... + 0 0 0 0 0 0 0 0 0 1 ; ... + 1 0 0 1 0 0 0 1 0 1 ; ... + 1 0 0 1 0 0 1 0 0 1 ; ... + 0 0 0 0 0 0 0 0 0 1 ; ... + 1 1 1 1 0 1 1 1 1 0 ; ... + 0 0 0 0 0 0 0 0 0 1 ; ... + 0 0 0 0 0 0 0 0 0 1 ; ... + 0 0 0 0 0 0 0 0 0 1 ; ... + 0 0 0 0 0 0 0 0 0 1 ; ... + 0 0 0 0 0 0 0 0 0 1 ]; + + ctype = 'int32'; + profile = 'std'; + qc = 1.0; +else + if nargin < 3 + ctype = 'int16'; + end + if nargin < 4 + fs_inout = ones(length(fs_in), length(fs_out)); + end + if nargin < 5 + profile = ''; + end + if nargin < 6 + qc = 1.0; + end +end + +sio = size(fs_inout); +if (length(fs_in) ~= sio(1)) || (length(fs_out) ~= sio(2)) + error('Sample rates in/out matrix size mismatch!'); +end + +%% Exported coefficients type int16, int24, int32, float + +switch ctype + case 'int16' + coef_label = 'int16'; + coef_ctype = 'int16_t'; + coef_bits = 16; + coef_bytes = 2; + case 'int24' + coef_label = 'int24'; + coef_ctype = 'int32_t'; + coef_bits = 24; + coef_bytes = 4; + case 'int32' + coef_label = 'int32'; + coef_ctype = 'int32_t'; + coef_bits = 32; + coef_bytes = 4; + case 'float' + coef_label = 'float'; coef_ctype = 'float'; + coef_bits = 24; + coef_bytes = 4; + otherwise + error('Request for incorrect coefficient type'); +end +data_bytes = 4; + +hdir = mkdir_check('include'); +rdir = mkdir_check('reports'); + +%% Find fractional conversion factors +nfsi = length(fs_in); +nfso = length(fs_out); +l_2s = zeros(2, nfsi, nfso); +m_2s = zeros(2, nfsi, nfso); +mops_2s = zeros(2, nfsi, nfso); +pb_2s = zeros(2,nfsi, nfso); +sb_2s = zeros(2,nfsi, nfso); +taps_2s = zeros(2,nfsi, nfso); +defs.fir_delay_size = 0; +defs.out_delay_size = 0; +defs.blk_in = 0; +defs.blk_out = 0; +defs.num_in_fs = nfsi; +defs.num_out_fs = nfso; +defs.stage1_times_max = 0; +defs.stage2_times_max = 0; +defs.stage_buf_size = 0; +h = 1; +for b = 1:nfso + for a = 1:nfsi + fs1 = fs_in(a); + fs2 = fs_out(b); + [l1, m1, l2, m2] = src_factor2_lm(fs1, fs2); + fs3 = fs1*l1/m1; + cnv1 = src_param(fs1, fs3, coef_bits, qc); + cnv2 = src_param(fs3, fs2, coef_bits, qc); + if (fs2 < fs1) + % When decimating 1st stage passband can be limited + % for wider transition band + f_pb = fs2*cnv2.c_pb; + cnv1.c_pb = f_pb/min(fs1,fs3); + end + if (fs2 > fs1) + % When interpolating 2nd stage passband can be limited + % for wider transition band + f_pb = fs1*cnv1.c_pb; + cnv2.c_pb = f_pb/min(fs2,fs3); + end + if fs_inout(a,b) > 0 || (a == b) + if cnv2.fs1-cnv2.fs2 > eps + % Allow half ripple for dual stage SRC parts + cnv1.rp = cnv1.rp/2; + cnv2.rp = cnv2.rp/2; + end + src1 = src_get(cnv1); + src2 = src_get(cnv2); + k = gcd(src1.blk_out, src2.blk_in); + stage1_times = src2.blk_in/k; + stage2_times = stage1_times*src1.blk_out/src2.blk_in; + defs.stage1_times_max = max(defs.stage1_times_max, stage1_times); + defs.stage2_times_max = max(defs.stage2_times_max, stage2_times); + l_2s(:,a,b) = [src1.L src2.L]; + m_2s(:,a,b) = [src1.M src2.M]; + mops_2s(:,a,b) = [src1.MOPS src2.MOPS]; + pb_2s(:,a,b) = [round(1e4*src1.c_pb) round(1e4*src2.c_pb)]; + sb_2s(:,a,b) = [round(1e4*src1.c_sb) round(1e4*src2.c_sb)]; + taps_2s(:,a,b) = [src1.filter_length src2.filter_length]; + defs.fir_delay_size = max(defs.fir_delay_size, src1.fir_delay_size); + defs.out_delay_size = max(defs.out_delay_size, src1.out_delay_size); + defs.blk_in = max(defs.blk_in, src1.blk_in); + defs.blk_out = max(defs.blk_out, src1.blk_out); + defs.fir_delay_size = max(defs.fir_delay_size, src2.fir_delay_size); + defs.out_delay_size = max(defs.out_delay_size, src2.out_delay_size); + defs.blk_in = max(defs.blk_in, src2.blk_in); + defs.blk_out = max(defs.blk_out, src2.blk_out); + defs.stage_buf_size = max(defs.stage_buf_size, src1.blk_out*stage1_times); + src_export_coef(src1, coef_label, coef_ctype, hdir, profile); + src_export_coef(src2, coef_label, coef_ctype, hdir, profile); + end + end +end + +%% Export modes table +defs.sum_filter_lengths = src_export_table_2s(fs_in, fs_out, l_2s, m_2s, ... + pb_2s, sb_2s, taps_2s, coef_label, coef_ctype, ... + 'sof/audio/coefficients/src/', hdir, profile); +src_export_defines(defs, coef_label, hdir, profile); + +%% Print 2 stage conversion factors +fn = sprintf('%s/src_2stage.txt', rdir); +fh = fopen(fn,'w'); +fprintf(fh,'\n'); +fprintf(fh,'Dual stage fractional SRC: Ratios\n'); +fprintf(fh,'%8s, ', 'in \ out'); +for b = 1:nfso + fprintf(fh,'%12.1f, ', fs_out(b)/1e3); +end +fprintf(fh,'\n'); +for a = 1:nfsi + fprintf(fh,'%8.1f, ', fs_in(a)/1e3); + for b = 1:nfso + cstr = print_ratios(l_2s, m_2s, a, b); + fprintf(fh,'%12s, ', cstr); + end + fprintf(fh,'\n'); +end +fprintf(fh,'\n'); + +%% Print 2 stage MOPS +fprintf(fh,'Dual stage fractional SRC: MOPS\n'); +fprintf(fh,'%8s, ', 'in \ out'); +for b = 1:nfso + fprintf(fh,'%8.1f, ', fs_out(b)/1e3); +end +fprintf(fh,'\n'); +for a = 1:nfsi + fprintf(fh,'%8.1f, ', fs_in(a)/1e3); + for b = 1:nfso + mops = sum(mops_2s(:,a,b)); + if sum(l_2s(:,a,b)) < eps + mops_str = 'x'; + else + mops_str = sprintf('%.2f', mops); + end + fprintf(fh,'%8s, ', mops_str); + end + fprintf(fh,'\n'); +end +fprintf(fh,'\n'); + +%% Print 2 stage MOPS per stage +fprintf(fh,'Dual stage fractional SRC: MOPS per stage\n'); +fprintf(fh,'%10s, ', 'in \ out'); +for b = 1:nfso + fprintf(fh,'%10.1f, ', fs_out(b)/1e3); +end +fprintf(fh,'\n'); +for a = 1:nfsi + fprintf(fh,'%10.1f, ', fs_in(a)/1e3); + for b = 1:nfso + mops = mops_2s(:,a,b); + if sum(l_2s(:,a,b)) < eps + mops_str = 'x'; + else + mops_str = sprintf('%.2f+%.2f', mops(1), mops(2)); + end + fprintf(fh,'%10s, ', mops_str); + end + fprintf(fh,'\n'); +end +fprintf(fh,'\n'); + +fprintf(fh,'Coefficient RAM %.1f kB\n', ... + defs.sum_filter_lengths*coef_bytes/1024); +fprintf(fh,'Max. data RAM %.1f kB\n', ... + (defs.fir_delay_size + defs.out_delay_size+defs.stage_buf_size) ... + * data_bytes/1024); + +fprintf(fh,'\n'); +fclose(fh); +type(fn); + +end + +function d = mkdir_check(d) +if exist(d) ~= 7 + mkdir(d); +end +end + +function cstr = print_ratios(l_2s, m_2s, a, b) +l1 = l_2s(1,a,b); +m1 = m_2s(1,a,b); +l2 = l_2s(2,a,b); +m2 = m_2s(2,a,b); +if l1+m2+l2+m2 == 0 + cstr = 'x'; +else + if m2 == 1 + if l2 == 1 + cstr2 = ''; + else + cstr2 = sprintf('*%d', l2); + end + else + cstr2 = sprintf('*%d/%d', l2, m2); + end + if m1 == 1 + cstr1 = sprintf('%d', l1); + else + cstr1 = sprintf('%d/%d', l1, m1); + end + cstr = sprintf('%s%s', cstr1, cstr2); +end +end diff --git a/tune/src/src_get.m b/tune/src/src_get.m new file mode 100644 index 0000000..4520af9 --- /dev/null +++ b/tune/src/src_get.m @@ -0,0 +1,254 @@ +function src = src_get(cnv) + +% src_get - calculate coefficients for a src +% +% src = src_get(cnv); +% +% cnv - src parameters +% src - src result +% + +% 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 +% + +use_firpm = 0; +use_remez = 0; +use_kaiser = 0; +switch lower(cnv.design) + case 'firpm' + if exist('OCTAVE_VERSION', 'builtin') + use_remez = 1; + else + use_firpm = 1; + end + case 'kaiser' + use_kaiser = 1; + otherwise + error('Unknown FIR design request!'); +end + +if abs(cnv.fs1-cnv.fs2) < 1 + %% Return minimum needed for scripts to work + src.L=1; + src.M=1; + src.odm=1; + src.idm=1; + src.MOPS=0; + src.c_pb = 0; + src.c_sb = 0; + src.fir_delay_size = 0; + src.out_delay_size = 0; + src.blk_in = 1; + src.blk_out = 1; + src.gain = 1; + src.filter_length = 0; + return +end + +%% fractional SRC parameters +[L, M] = src_factor1_lm(cnv.fs1, cnv.fs2); +src.L = L; +src.M = M; +src.num_of_subfilters = L; +[src.idm src.odm] = src_find_l0m0(src.L, src.M); +min_fs = min(cnv.fs1,cnv.fs2); +src.f_pb = min_fs*cnv.c_pb; +src.f_sb = min_fs*cnv.c_sb; +src.c_pb = cnv.c_pb; +src.c_sb = cnv.c_sb; +src.fs1 = cnv.fs1; +src.fs2 = cnv.fs2; +src.rp = cnv.rp; +src.rs = cnv.rs; + +%% FIR PM design +fs3 = src.L*cnv.fs1; +f = [src.f_pb src.f_sb]; +a = [1 0]; + +% Kaiser fir is not equiripple, can allow more ripple in the passband +% lower freq. +dev = [(10^(cnv.rp/20)-1)/(10^(cnv.rp/20)+1) 10^(-cnv.rs/20)]; +[kn0, kw, kbeta, kftype] = kaiserord(f, a, dev, fs3); +if use_firpm + %dev = [(10^(cnv.rp/20)-1)/(10^(cnv.rp/20)+1) 10^(-cnv.rs/20)]; + [n0,fo,ao,w] = firpmord(f,a,dev,fs3); + n0 = round(n0*0.95); % Decrease order to 95% +end +if use_remez + n0 = round(kn0*0.60); % Decrease order to 70% + fo = [0 2*f/fs3 1]; + ao = [1 1 0 0]; + w = [1 dev(1)/dev(2)]; +end +if use_kaiser + n0 = round(kn0*0.70); % Decrease order to 70% +end + +% Constrain filter length to be a suitable multiple. Multiple of +% interpolation factor ensures that all subfilters are equal length. +% Make each sub-filter order multiple of N helps in making efficient +% implementation. +nfir_increment = src.num_of_subfilters * cnv.filter_length_mult; + +nsf0 = (n0+1)/nfir_increment; +if nsf0 > floor(nsf0) + n = (floor(nsf0)+1)*nfir_increment-1; +else + n = n0; +end + +src.subfilter_length = (n+1)/src.num_of_subfilters; +src.filter_length = n+1; +nfir = n; +f_sb = linspace(src.f_sb, fs3/2, 500); +stopband_ok = 0; +need_to_stop = 0; +delta = 100; +n_worse = 0; +n_worse_max = 20; +n = 1; +n_max = 100; +dn = ones(1, n_max)*1000; +fn = zeros(1, n_max); +while (stopband_ok) == 0 && (n < n_max) + if use_firpm || use_remez + if nfir > 1800 + b = fir1(nfir, kw, kftype, kaiser(nfir+1, kbeta)); + else + if use_firpm + b = firpm(nfir,fo,ao,w); + else + b = remez(nfir,fo,ao,w); + end + end + else + b = fir1(nfir, kw, kftype, kaiser(nfir+1, kbeta)); + end + m_b = max(abs(b)); + %sref = 2^(cnv.coef_bits-1); + %bq = round(b*sref/m_b)*m_b/sref; + bq = b; + h_sb = freqz(bq, 1, f_sb, fs3); + m_sb = 20*log10(abs(h_sb)); + delta_prev = delta; + delta = cnv.rs+max(m_sb); + fprintf('Step=%3d, Delta=%f dB, N=%d\n', n, delta, nfir); + dn(n) = delta; + fn(n) = nfir; + if delta < 0 + stopband_ok = 1; + else + if delta_prev < delta + n_worse = n_worse+1; + if n_worse > n_worse_max + need_to_stop = 1; % No improvement, reverse + idx = find(dn == min(dn), 1, 'first'); + nfir = fn(idx); + else + nfir = nfir + nfir_increment; + end + else + if need_to_stop == 0 + nfir = nfir + nfir_increment; + else + stopband_ok = 1; % Exit after reverse + fprintf('Warning: Filter stop band could not be '); + fprintf('reached.\n', cnv.coef_bits); + end + end + end + n = n + 1; +end + +f_p = linspace(0, fs3/2, 2000); +m_db = 20*log10(abs(freqz(bq, 1, f_p, fs3))); +i_pb = find(f_p < src.f_pb); +g_pb = max(m_db(i_pb)); +g_att_lin = 10^(-g_pb/20); + +if 1 + p_ymin = floor((-cnv.rs-50)/10)*10; + p_ymax = 10; + figure; + clf; + plot(f_p, m_db); + grid on; + xlabel('Frequency (Hz)'); + ylabel('Magnitude (dB)'); + axis([0 fs3/2 p_ymin p_ymax]); + hold on + plot([src.f_sb src.f_sb],[p_ymin p_ymax], 'r--'); + plot([0 fs3],[-cnv.rs -cnv.rs], 'r--'); + hold off + axes('Position', [ 0.58 0.7 0.3 0.2]); + box on; + plot(f_p(i_pb), m_db(i_pb)); + axis([0 src.f_pb -cnv.rp-0.01 cnv.rp+0.01]); + hold on + plot([0 src.f_pb], +0.5*cnv.rp*ones(1,2), 'r--'); + plot([0 src.f_pb], -0.5*cnv.rp*ones(1,2), 'r--'); + hold off + grid on; + box off; +end + +src.subfilter_length = ceil((nfir+1)/src.num_of_subfilters); +src.filter_length = src.subfilter_length*src.num_of_subfilters; +src.b = zeros(src.filter_length,1); +src.gain = 1; +src.b(1:nfir+1) = b*src.L*g_att_lin; +m = max(abs(src.b)); +gmax = (32767/32768)/m; +maxshift = floor(log(gmax)/log(2)); +src.b = src.b * 2^maxshift; +src.gain = 1/2^maxshift; +src.shift = maxshift; + +%% Reorder coefficients +if 1 + src.coefs = src.b; +else + src.coefs = zeros(src.filter_length,1); + for n = 1:src.num_of_subfilters + i11 = (n-1)*src.subfilter_length+1; + i12 = i11+src.subfilter_length-1; + src.coefs(i11:i12) = src.b(n:src.num_of_subfilters:end); + end +end + +src.halfband = 0; +src.blk_in = M; +src.blk_out = L; +src.MOPS = cnv.fs1/M*src.filter_length/1e6; +src.fir_delay_size = src.subfilter_length + ... + (src.num_of_subfilters-1)*src.idm + src.blk_in; +src.out_delay_size = (src.num_of_subfilters-1)*src.odm + 1; + +end diff --git a/tune/src/src_param.m b/tune/src/src_param.m new file mode 100644 index 0000000..49e935c --- /dev/null +++ b/tune/src/src_param.m @@ -0,0 +1,87 @@ +function cnv = src_param(fs1, fs2, coef_bits, q) + +% src_param - get converter parameters +% +% cnv = src_param(fs1, fs2, coef_bits, q) +% +% fs1 - input rate +% fs2 - output rate +% coef_bits - word length identifier +% q - quality scale filter bandwidth and stopband attenuation, +% 1 is default +% + +% 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 +% + +% Note: firpm design with rs 93 dB, or kaiser with 69 dB rs is fair quality. +% Both give about -80 dB THD+N with 24 bit coefficients. With 16 bit +% coefficients THD+N is limited to about -76 dB. + +if nargin < 4 + q = 1.0; +end + +%% Copy input parameters +cnv.fs1 = fs1; +cnv.fs2 = fs2; +cnv.coef_bits = coef_bits; + +%% FIR design +cnv.design = 'kaiser'; % Use firpm or kaiser + +%% Default SRC quality +cnv.c_pb = q * 22/48; % Gives 22 kHz BW @ 48 kHz +cnv.c_sb = 0.5; % Start stopband at Fs/2 +cnv.rs = 70; % Stopband attenuation in dB +cnv.rp = 0.1; % Passband ripple in dB +cnv.rp_tot = 0.1; % Max +/- passband ripple allowed, used in test script only + +%% Constrain sub-filter lengths. Make subfilters lengths multiple of four +% is a good assumption for processors. +cnv.filter_length_mult = 4; + +%% Exceptions for quality +if min(fs1, fs2) > 80e3 + cnv.c_pb = 30e3/min(fs1, fs2); % 30 kHz BW for > 80 kHz +end + +%% Sanity checks +if cnv.c_pb > 0.49 + error('Too wide passband'); +end +if cnv.c_pb < 0.15 + error('Too narrow passband'); +end +if cnv.rs > 160 + error('Too large stopband attenuation'); +end +if cnv.rs < 40 + error('Too low stopband attenuation'); +end diff --git a/tune/src/src_std_int32.m b/tune/src/src_std_int32.m new file mode 100644 index 0000000..125ae96 --- /dev/null +++ b/tune/src/src_std_int32.m @@ -0,0 +1,30 @@ +% 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 +% + +src_generate(); diff --git a/tune/src/src_tiny_int16.m b/tune/src/src_tiny_int16.m new file mode 100644 index 0000000..be310fb --- /dev/null +++ b/tune/src/src_tiny_int16.m @@ -0,0 +1,41 @@ +% 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 +% + +fs1 = [8e3 16e3 24e3 32e3 44.1e3 48e3]; +fs2 = [8e3 16e3 24e3 32e3 44.1e3 48e3]; +fsm = [0 0 0 0 0 1; ... + 0 0 0 0 0 1; ... + 0 0 0 0 0 1; ... + 0 0 0 0 0 1; ... + 0 0 0 0 0 1; ... + 1 1 1 1 1 0; ... + ]; +fmt = 'int16'; profile = 'tiny'; q = 0.75; + +src_generate(fs1, fs2, fmt, fsm, profile, q);