Code: Select all
#include <aucommon.h>Code: Select all
#include <audiofs.h>Code: Select all
function bIntOut=fir2vsdsp(b,bits,upsampleBy,downsampleBy,varargin)
% fir2vsdsp(B,BITS,UPSAMPLEBY,DOWNSAMPLEBY,...);
% Creates C code for a BITS-bit filter NAME that first upsamples
% data by UPSAMPLEBY, then filters with Fir filter B, then
% downsamples by DOWNSAMPLEBY. Prints the filter structure as C code.
% VS1005 header file <fir.h> contains definitions required for running
% the filters.
%
% Returns:
% The vector B rounded to the accuracy that could be represented by its
% internal integer accuracy.
%
% Parameters:
% B - Fir filter
% BITS - Number of bits. Currently only 16 is supported.
% UPSAMPLEBY - How many times to upsample before filtering
% DOWNSAMPLEBY - How many times to downsample after filtering
% ... - Varargs options may contain some of the following fields:
% NAME - Name of the filter. If more than 1 channels are required,
% many names can be specified. There must be at least one name.
% '-split' - If defined in the argument list, make a split filter
% that outputs twice the input data (e.g. low-pass and high-pass path).
% '-rate',FS - If defined, followed by FS, the samplerate of the source
% signal.
% Example: fir2vsdsp(b,16,2,3,'myFilter','-rate',48000). In this case
% a frequency response of the filter is plotted. If a split filter,
% frequency response of both paths is plotted.
% '-magnify' - Show Y scale of -2 to 2 dB instead of the default.
% '-file',NAME - Write source code to file NAME.c and NAME.h.
% '-append' - If writing to a file, append the code the files instead of
% creating the files.
%
% 2015-02-18 VLSI Solution Oy
canBeHalfBand = 0;
makeMod1 = 0;
makeMod2 = 0;
flags = 0;
centerOffset=0;
stdout = 1;
fsIn=-1;
magnify=0;
if nargin < 5,
error 'Required argument(s) missing';
end
legendStrings = {};
outFileNameBase='-';
coeffName='-';
plotName='-';
append=0;
splitSignal=0;
for i=1:length(varargin),
if fsIn == -2,
%if ~strcmp(typeinfo(varargin{i}), 'scalar'),
if ~isnumeric(varargin{i}) || length(varargin{i}) ~= 1
error 'Parameter for \"-rate\" must be scalar';
end
fsIn = varargin{i};
elseif strcmp(outFileNameBase, '+'),
if ~ischar(varargin{i})
error 'Name parameter must be string';
end
outFileNameBase = varargin{i};
elseif strcmp(varargin{i}, '-rate'),
fsIn = -2;
elseif strcmp(varargin{i}, '-split'),
splitSignal=1;
elseif strcmp(varargin{i}, '-magnify'),
magnify=1;
elseif strcmp(varargin{i}, '-append'),
append=1;
elseif strcmp(varargin{i}, '-file'),
outFileNameBase='+';
elseif strcmp(coeffName,'-'),
%if ~strcmp(typeinfo(varargin{i}), 'string'),
if ~ischar(varargin{i})
error 'Name parameter must be string';
end
coeffName=varargin{i};
fileIndex(1) = i;
else
if ~ischar(varargin{i})
error 'Name parameter must be string';
end
fileIndex = [fileIndex i];
end
end
if strcmp(outFileNameBase, '+'),
error 'No parameter given for \"-file\"';
end
if fsIn == -2,
error 'No parameter given for \"-rate\"';
end
if strcmp(coeffName, '-')
error 'No filter name specified';
end
if bits ~= 16,
error 'Number of bits must be 16';
end
% Make vector always be in known direction
[bx,by]=size(b);
if bx < by,
b=b';
end
% Is filter too short?
len=length(b);
if len < 3,
b = [b(1:end)' zeros(1,3-length(b))]';
len=3;
end
shiftLeft = -1;
% Take care of scale. This intentionally allows max(bInt) to be
% one larger than what can be represented. That is taken care of
% later.
b=b*upsampleBy;
scale = 1;
while scale,
bInt=round(b*(2^(bits-1-shiftLeft)));
bIntOut=bInt/(2^(bits-1-shiftLeft));
if (max(abs(bInt))) > 2^(bits-1),
shiftLeft = shiftLeft + 1;
elseif max(abs(bInt)) < 2^(bits-2),
shiftLeft = shiftLeft - 1;
else
scale = 0;
end
end
% Strip zero elements from beginning and end
while bInt(1) == 0.0 && length(bInt) > 3,
bInt=bInt(2:end);
b=b(2:end);
end
while bInt(length(bInt)) == 0.0 && length(bInt) > 3,
bInt=bInt(1:end-1);
b=b(1:end-1);
end
len=length(bInt);
fsHigh = fsIn*upsampleBy;
fsOut = fsHigh/downsampleBy;
bwIn = floor(fsIn/2);
bwHigh = floor(fsHigh/2);
bwOut = floor(fsOut/2);
memLen=len;
% MAKEMODF bufSize-1
makeMod1 = hex2dec('8000') + (len-1);
if mod(length(bInt),2) == 1,
canBeHalfBand = 1;
% Check zeros
zerosum=sum(abs(bInt(2:2:end)))-abs(bInt((len+1)/2));
if zerosum > 0.5,
canBeHalfBand = 0;
end
% Check upsample and downsample ratios
if upsampleBy == 1 && downsampleBy == 1,
% up/down ok
elseif upsampleBy == 1 && downsampleBy == 2,
% up/down ok
elseif upsampleBy == 2 && downsampleBy == 1,
% up/down ok
else,
canBeHalfBand = 0;
end
% Check center value
if bInt((len+1)/2) ~= 2^(bits-2-shiftLeft)*upsampleBy,
canBeHalfBand = 0;
end
% Check symmetry
if sum(abs(bInt(1:(len-1)/2) - bInt(end:-1:(len+3)/2))) >= 0.5,
canBeHalfBand = 0;
end
% Check maximum size. This has to do with the VSDSP modifier register
% limitations. Check also minimum size.
if memLen > 256 || memLen < 7,
canBeHalfBand = 0;
end
% Split filters have only downsampling by 2 supported
if splitSignal && (upsampleBy ~= 1 || downsampleBy ~= 2),
canBeHalfBand = 0;
end
% If everything is ok, make it half-band
if canBeHalfBand,
flags = bitor(flags,1);
bInt=bInt(1:2:len);
if upsampleBy == 2,
len = ceil(len/2);
memLen = len;
end
if len < 64,
if upsampleBy == 1 && downsampleBy == 1,
memLen = memLen+1; % Required by the quick HalfBand filter
end
% Make modifiers: one for normal samples, one for center
% MAKEMOD step bufSize-1
makeMod1 = hex2dec('2000') + 1*(2^6) + (memLen-1);
makeMod2 = hex2dec('2000') + 2*(2^6) + (memLen-1);
else
memLen = ceil(memLen/64)*64;
% MAKEMOD64 step bufSize/64-1
makeMod1 = hex2dec('4000') + 1*(2^6) + (memLen/64-1);
makeMod2 = hex2dec('4000') + 2*(2^6) + (memLen/64-1);
end
end
end
% Check that not a split filter for which there is no function
if splitSignal,
if canBeHalfBand,
if upsampleBy~=1 || downsampleBy~=2,
error 'Cannot make half-band split filter that is not decimation-by-2';
end
else
if upsampleBy~=1 || downsampleBy~=1,
error 'Cannot make normal split filter that has interpolation/decimation';
elseif mod(len,2) == 0
error 'Split filter length must be odd';
end
end
flags = bitor(flags,2);
end
bInt = -bInt; % Solves the case where max = 32768
% Finally, scale down if necessary
if max(bInt) >= 2^(bits-1),
bInt = bInt/2;
shiftLeft = shiftLeft + 1;
end
funcName='***ERROR***';
centerOffset = ceil(len/2);
plotName = 'Freq. resp. of ';
if bits==16,
if upsampleBy==1 && downsampleBy==1,
coeffLen = length(bInt);
if canBeHalfBand,
if memLen < 64,
funcName='FirFilter16x16HB';
plotName = sprintf('%s%s', plotName, 'half-band filter');
else
funcName='FirFilter16x16HBLong';
plotName = sprintf('%s%s', plotName, 'long half-band');
end
else
if splitSignal,
funcName='FirFilter16x16NSplit';
plotName = sprintf('%s%s', plotName, 'basic split');
else
funcName='FirFilter16x16N';
plotName = sprintf('%s%s', plotName, 'basic');
end
end
else
if canBeHalfBand,
if upsampleBy==1 && downsampleBy==2,
coeffLen = (memLen+1)/2;
centerOffset = coeffLen+1;
if splitSignal,
plotName = sprintf('%s%s', plotName, 'half-band split');
funcName='FirFilter16x16HBUp1Dn2Split';
else
plotName = sprintf('%s%s', plotName, 'half-band');
funcName='FirFilter16x16HBUp1Dn2';
end
elseif upsampleBy==2 && downsampleBy==1,
coeffLen = memLen;
plotName = sprintf('%s%s', plotName, 'half-band');
funcName='FirFilter16x16HBUp2Dn1';
else
error 'Internal error #1';
end
else
memLen = ceil(memLen/upsampleBy);
if memLen < 3,
memLen = 3;
end
coeffLen = memLen*upsampleBy;
len = memLen;
% MAKEMODF bufSize-1
makeMod1 = hex2dec('8000') + (memLen-1);
funcName='FirFilter16x16NUpDn';
plotName = sprintf('%s%s', plotName, 'basic');
end
end
end
plotName = sprintf('%s filter, %d Hz', plotName, fsIn);
if upsampleBy > 1,
plotName = sprintf('%s, intep. by %d->%d Hz', plotName, upsampleBy, fsHigh);
end
if downsampleBy > 1,
plotName = sprintf('%s, decim. by %d->%d Hz', plotName, downsampleBy, fsOut);
end
if ~strcmp(outFileNameBase, '-'),
fileName = sprintf('%s.c', outFileNameBase);
if append,
oFpC = fopen(fileName, 'r+');
else
oFpC = fopen(fileName, 'w');
end
if oFpC < 0,
error 'Could not open output C file';
end
fileName = sprintf('%s.h', outFileNameBase);
if append,
oFpH = fopen(fileName, 'r+');
else
oFpH = fopen(fileName, 'w');
end
if oFpH < 0,
fclose(oFpC);
error 'Could not open output H file';
end
else
oFpC = stdout;
end
if oFpC ~= stdout,
if ~append,
fprintf(oFpC, '/*\n This is an automatically generated file by VLSI Solutions\n');
fprintf(oFpC, ' Octave program fir2vsdsp.m.\n\n');
fprintf(oFpC, ' DO NOT MODIFY unless you ABSOLUTELY know what you are doing!\n*/\n\n');
fprintf(oFpC, '#include <vstypes.h>\n#include <fir.h>\n');
fprintf(oFpC, '#include \''%s.h\''\n', outFileNameBase);
fprintf(oFpH, '/*\n This is an automatically generated file by VLSI Solutions\n');
fprintf(oFpH, ' Octave program fir2vsdsp.m.\n\n');
fprintf(oFpH, ' DO NOT MODIFY unless you ABSOLUTELY know what you are doing!\n*/\n\n');
fprintf(oFpH, '#ifndef __FIR_%s_H__\n', upper(outFileNameBase));
fprintf(oFpH, '#define __FIR_%s_H__\n\n', upper(outFileNameBase));
fprintf(oFpH, '#include <vstypes.h>\n#include <fir.h>\n\n');
else
% Seek to end of .h file, except just before the final #endif,
% which will be overwritten. This method will work with both
% Unix and Windows line feed characters.
fseek(oFpC, 0, 'eof');
fseek(oFpH, -10, 'eof');
val = -1;
% Look for the '#' of the final #endif
while val ~= 35,
val = fread(oFpH, 1, 'uchar');
end
% Then take one char back
fseek(oFpH, -1, 'cof');
end
end
fprintf(oFpC, '\n/* Filter structures for FIR filter \''%s\'' */\n', coeffName);
fprintf(oFpC, '\nconst s_int16 %sCoeff[%d] = {', coeffName, coeffLen);
for i=0:coeffLen-1,
if (mod(i,8) == 0),
fprintf(oFpC, '\n ');
end
if i>=length(bInt),
t=0;
else
t=bInt(i+1);
end
fprintf(oFpC, ' %6d,', t);
end
fprintf(oFpC, '\n};\n');
onePerDownsampleBy = ceil(32768/downsampleBy);
for i=1:length(fileIndex),
name = varargin{fileIndex(i)};
fprintf(oFpC, 's_int16 __align __mem_y %sMem[%d] = {0};\n', name, memLen);
fprintf(oFpC, 'struct FirFilter __mem_y %s = {\n', name);
fprintf(oFpC, ' %s, 0x%04x/*flags*/, %d/*len*/,\n', funcName, flags, len);
fprintf(oFpC, ' %sCoeff, %sMem, %sMem+%d,\n', coeffName, name, name, centerOffset);
fprintf(oFpC, ' 0x%04x/*mod1*/, 0x%04x/*mod2*/, %d/*shL*/,\n', makeMod1, makeMod2, shiftLeft);
fprintf(oFpC, ' %d/*up*/, %d/*down*/, %d/*dnPhase*/, 0x%04x/*1/dn*/\n', upsampleBy, downsampleBy, downsampleBy, onePerDownsampleBy);
fprintf(oFpC, '};\n');
if oFpC ~= stdout,
fprintf(oFpH, 'extern struct FirFilter __mem_y %s;\n', name);
end
end
fprintf(oFpC, '\n');
if oFpC ~= stdout,
fprintf(oFpH, '\n#endif\n');
end
if oFpC ~= stdout,
fclose(oFpC);
fclose(oFpH);
end
if fsIn >= 0,
x=(0:floor(fsHigh/2)-1);
bSpl=bIntOut*upsampleBy;
bSpl(floor((length(bIntOut)+1)/2)) = bSpl(floor((length(bIntOut)+1)/2)) - 1;
y=[bIntOut(1:end)'/upsampleBy zeros(1,fsHigh-length(b))];
y=20*log10(abs(fft(y)));
y=y(1:length(x));
ySpl=[bSpl(1:end)'/upsampleBy zeros(1,fsHigh-length(bSpl))];
ySpl=20*log10(abs(fft(ySpl)));
ySpl=ySpl(1:length(x));
plot(x,y);
legendStrings = [legendStrings,'Base filter frequency response'];
hold on
if splitSignal,
plot(x,ySpl);
legendStrings = [legendStrings,'Split filter frequency response'];
end
if bwIn < bwOut,
stepSize=floor((bwOut-bwIn)/16);
plot(x(bwIn:stepSize:bwOut),y(bwIn:stepSize:bwOut),'ro');
legendStrings = [legendStrings,'Base upsampling aliasing;'];
if splitSignal,
plot(x(bwIn:stepSize:bwOut),ySpl(bwIn:stepSize:bwOut));
legendStrings = [legendStrings,'Split upsampling aliasing;'];
end
end
if downsampleBy > 1,
x2=x;
for i=1:floor((bwHigh-20)/bwOut),
startIdx = i*bwOut;
endIdx = min([(i+1)*bwOut-1 length(x2)]);
if mod(i,2),
x2(startIdx:endIdx) = bwOut-1:-1:0;
else,
x2(startIdx:endIdx) = 0:bwOut-1;
end
end
plot(x2(bwOut:endIdx),y(bwOut:endIdx));
legendStrings = [legendStrings,'Base downsampling aliasing;'];
if splitSignal,
plot(x2(bwOut:endIdx),ySpl(bwOut:endIdx));
legendStrings = [legendStrings,'Split downsampling aliasing;'];
end
end
if upsampleBy > 1,
inLimit = sprintf('Input bandwidth %d Hz;', bwIn);
plot(linspace(fsIn/2, fsIn/2, 26),[-100:4:0],'k+-')
legendStrings = [legendStrings,inLimit];
end
if downsampleBy > 1,
outLimit = sprintf('Output bandwidth %d Hz;', bwOut);
plot(linspace(fsOut/2, fsOut/2, 26),[-100:4:0],'kx-');
legendStrings = [legendStrings,outLimit];
end
hold off
if magnify,
axis([0 fsHigh/2 -2 2]);
else
axis([0 fsHigh/2 -120 10+max(y)]);
end
title(plotName);
xlabel('frequency / Hz');
ylabel('gain / dB');
grid on;
legend(legendStrings);
end