1function success=src_export_coef(src, ctype, vtype, hdir, profile)
2
3% src_export_coef - Export FIR coefficients
4%
5% success=src_export_coef(src, ctype, hdir, profile)
6%
7% src     - src definition struct
8% ctype   - 'float','int32', or 'int24'
9% vtype   - 'float','int32_t'
10% hdir    - directory for header files
11% profile - string to append to filename
12%
13
14% SPDX-License-Identifier: BSD-3-Clause
15%
16% Copyright (c) 2016-2020, Intel Corporation. All rights reserved.
17%
18% Author: Seppo Ingalsuo <seppo.ingalsuo@linux.intel.com>
19
20if nargin < 5
21        profile = '';
22end
23
24if src.L == src.M
25        success = 0;
26else
27        pbi =  round(src.c_pb*1e4);
28        sbi =  round(src.c_sb*1e4);
29        if isempty(profile)
30                hfn = sprintf('%s/src_%s_%d_%d_%d_%d.h', ...
31                        hdir, ctype, src.L, src.M, pbi, sbi);
32        else
33                hfn = sprintf('%s/src_%s_%s_%d_%d_%d_%d.h', ...
34                        hdir, profile, ctype, src.L, src.M, pbi, sbi);
35        end
36        vfn = sprintf('src_%s_%d_%d_%d_%d_fir', ctype, src.L, src.M, pbi, sbi);
37        sfn = sprintf('src_%s_%d_%d_%d_%d', ctype, src.L, src.M, pbi, sbi);
38
39        fprintf('Exporting %s ...\n', hfn);
40        fh = fopen(hfn, 'w');
41	y = datestr(now(), 'yyyy');
42
43	fprintf(fh, '/* SPDX-License-Identifier: BSD-3-Clause\n');
44	fprintf(fh, ' *\n');
45	fprintf(fh, ' * Copyright(c) %s Intel Corporation. All rights reserved.\n', y);
46	fprintf(fh, ' *\n');
47	fprintf(fh, ' */\n');
48	fprintf(fh, '\n');
49	fprintf(fh, '#include <sof/audio/src/src.h>\n');
50	fprintf(fh, '#include <stdint.h>\n');
51	fprintf(fh, '\n');
52
53        switch ctype
54                case 'float'
55                        fprintf(fh, 'const %s %s[%d] = {\n', ...
56                                vtype, vfn, src.filter_length);
57                        fprintf(fh,'\t%16.9e', src.coefs(1));
58                        for n=2:src.filter_length
59                                fprintf(fh, ',\n');
60                                fprintf(fh,'\t%16.9e', src.coefs(n));
61                        end
62                        fprintf(fh,'\n\n};');
63                case 'int32'
64			print_int_coef(src, fh, vtype, vfn, 32);
65                case 'int24'
66			print_int_coef(src, fh, vtype, vfn, 24);
67                case 'int16'
68			print_int_coef(src, fh, vtype, vfn, 16);
69                otherwise
70                        error('Unknown type %s !!!', ctype);
71        end
72
73        fprintf(fh, '\n');
74        switch ctype
75                case 'float'
76                        fprintf(fh, 'struct src_stage %s = {\n', sfn);
77			fprintf(fh, '\t%d, %d, %d, %d, %d, %d, %d, %d, %f,\n\t%s};\n', ...
78                                src.idm, src.odm, src.num_of_subfilters, ...
79                                src.subfilter_length, src.filter_length, ...
80                                src.blk_in, src.blk_out, src.halfband, ...
81                                src.gain, vfn);
82                case { 'int16' 'int24' 'int32' }
83                        fprintf(fh, 'struct src_stage %s = {\n', sfn);
84			fprintf(fh, '\t%d, %d, %d, %d, %d, %d, %d, %d, %d,\n\t%s};\n', ...
85                                src.idm, src.odm, src.num_of_subfilters, ...
86                                src.subfilter_length, src.filter_length, ...
87                                src.blk_in, src.blk_out, src.halfband, ...
88                                src.shift, vfn);
89                otherwise
90                        error('Unknown type %s !!!', ctype);
91        end
92        %fprintf(fh, '\n');
93        fclose(fh);
94        success = 1;
95end
96
97end
98
99function print_int_coef(src, fh, vtype, vfn, nbits)
100        fprintf(fh, 'const %s %s[%d] = {\n', ...
101                vtype, vfn, src.filter_length);
102
103        cint = coef_quant(src, nbits);
104        fprintf(fh,'\t%d', cint(1));
105        for n=2:src.filter_length
106                fprintf(fh, ',\n');
107                fprintf(fh,'\t%d', cint(n));
108        end
109        fprintf(fh,'\n\n};\n');
110end
111
112function cint = coef_quant(src, nbits)
113
114	sref = 2^(nbits-1);
115	pmax = sref-1;
116	nmin = -sref;
117
118        if nbits > 16
119                %% Round() is OK
120                cint0 = round(sref*src.coefs);
121        else
122                %% Prepare to optimize coefficient quantization
123                fs = max(src.fs1, src.fs2);
124                f = linspace(0, fs/2, 1000);
125
126                %% Test sensitivity for stopband and find the most sensitive
127                %  coefficients
128                sbf = linspace(src.f_sb,fs/2, 500);
129                n = src.filter_length;
130                psens = zeros(1,n);
131                bq0 = round(sref*src.coefs);
132                h = freqz(bq0/sref/src.L, 1, sbf, fs);
133                sb1 = 20*log10(sqrt(sum(h.*conj(h))));
134                for i=1:n
135                        bq = src.coefs;
136                        bq(i) = round(sref*bq(i))/sref;
137                        %tbq = bq; %tbq(i) = bq(i)+1;
138                        h = freqz(bq, 1, sbf, fs);
139                        psens(i) = sum(h.*conj(h));
140                end
141                [spsens, pidx] = sort(psens, 'descend');
142
143                %% Test quantization in the found order
144                %  The impact to passband is minimal so it's not tested
145                bi = round(sref*src.coefs);
146                bi0 = bi;
147                dl = -1:1;
148                nd = length(dl);
149                msb = zeros(1,nd);
150                for i=pidx
151                        bit = bi;
152                        for j=1:nd
153                                bit(i) = bi(i) + dl(j);
154                                h = freqz(bit, 1, sbf, fs);
155                                msb(j) = sum(h.*conj(h));
156                        end
157                        idx = find(msb == min(msb), 1, 'first');
158                        bi(i) = bi(i) + dl(idx);
159                end
160                h = freqz(bi/sref/src.L, 1, sbf, fs);
161                sb2 = 20*log10(sqrt(sum(h.*conj(h))));
162
163                %% Plot to compare
164                if 0
165                        f = linspace(0, fs/2, 1000);
166                        h1 = freqz(src.coefs/src.L, 1, f, fs);
167                        h2 = freqz(bq0/sref/src.L, 1, f, fs);
168                        h3 = freqz(bi/sref/src.L, 1, f, fs);
169                        figure;
170                        plot(f, 20*log10(abs(h1)), f, 20*log10(abs(h2)), f, 20*log10(abs(h3)));
171                        grid on;
172                        fprintf('Original = %4.1f dB, optimized = %4.1f dB, delta = %4.1f dB\n', ...
173                                sb1, sb2, sb1-sb2);
174                end
175                cint0 = bi;
176        end
177
178
179        %% Re-order coefficients for filter implementation
180        cint = zeros(src.filter_length,1);
181        for n = 1:src.num_of_subfilters
182                i11 = (n-1)*src.subfilter_length+1;
183                i12 = i11+src.subfilter_length-1;
184                cint(i11:i12) = cint0(n:src.num_of_subfilters:end);
185        end
186
187        %% Done check for no overflow
188        max_fix = max(cint);
189	min_fix = min(cint);
190	if (max_fix > pmax)
191		printf('Fixed point coefficient %d exceeded %d\n.', max_fix, pmax);
192		error('Something went wrong!');
193	end
194	if (min_fix < nmin)
195		printf('Fixed point coefficient %d exceeded %d\n.', min_fix, nmax);
196		error('Something went wrong!');
197        end
198
199
200end
201