// FIR filter desing by inverse DFT(FFT) and window function

clear;

FS=48000;             // FS should be even number
FIRLEN=255;           // FIRLEN should be odd number
SCALE16bit=32768;     // (2^16)/2=32768

// freq_response(1:FS/2) is FIR filter frequency response specification
freq_response=[ones(1:2000), zeros(1:4000), ones(1:4000), zeros(1:3000), 0.1*ones(1:3000), ones(1:2000), 0.999^(1:4000), zeros(1:2000)];
//freq_response=[zeros(1:4000), ones(1:20000)];
//freq_response=[ones(1:2000), zeros(1:4000), ones(1:4000), zeros(1:6000), ones(1:2000), zeros(1:6000)];

freq_response=[freq_response, freq_response(FS/2), mtlb_fliplr(freq_response(1:FS/2-1))];

//plot2d(freq_response, rect=[0, -0.2, 48000, 1.2], axesflag=1);
plot2d(20*log10(freq_response(1:FS/2)+1e-19), style=2, axesflag=1, rect=[0, -80, FS/2, 10]);

coef=real(fft(freq_response, 1));          // inverse FFT
coef=[coef(FS/2+1:FS), coef(1:FS/2)];

//plot2d(coef, rect=[0, -0.5, FS, 0.5], axesflag=1);
//plot2d(coef, rect=[FS/2-500, -0.5, FS/2+500, 0.5], axesflag=1);
//plot(20*log10(abs(fft([coef, zeros(1:FS-FIRLEN)], -1))+1e-19));

coef=coef(FS/2-(FIRLEN-1)/2+1:FS/2+(FIRLEN-1)/2+1);

//plot2d(coef, rect=[0, -0.5, FIRLEN, 0.5], axesflag=1);
//plot(20*log10(abs(fft([coef, zeros(1:FS-FIRLEN)], -1))+1e-19));

// generate window function coefficient
//win=sin(%pi/(FIRLEN-1)*(0:FIRLEN-1));                  // sine window
win=0.5-0.5*cos(2.0*%pi/(FIRLEN-1)*(0:FIRLEN-1));        // Hanning window
//win=0.54-0.46*cos(2.0*%pi/(FIRLEN-1)*(0:FIRLEN-1));    // Hamming window
//win=0.85-0.15*cos(2.0*%pi/(FIRLEN-1)*(0:FIRLEN-1));    // unknown trivial window

coef=coef.*win;            // apply window

//plot2d(coef, rect=[0, -0.5, FIRLEN, 0.5], axesflag=1);

designed_response=20*log10(abs(fft([coef, zeros(1:FS-FIRLEN)], -1))+1e-19);
designed_response=designed_response(1:FS/2);

//plot2d(designed_response, style=5, axesflag=1, rect=[0, -80, FS/2, 10]);

for i=1:FIRLEN
  mprintf("%3i %15.12f\n", i-1, coef(i));
end

//-------------------------------------------------------------
// FIR filter coefficient conversion (float --> 16bit fixed point)

coef_gain=1.1*max(abs(coef));
scaled_coef=coef./coef_gain;

//plot(scaled_coef);

coef16bit=fix(scaled_coef.*SCALE16bit);

//plot(coef16bit);

for i=1:FIRLEN
  mprintf("%3i %5i\n", i-1, coef16bit(i));
end
mprintf("coef_gain : %f\n", coef_gain);

filter_16bit_response=20*log10(abs(fft([coef16bit./SCALE16bit*coef_gain, zeros(1:FS-FIRLEN)], -1))+1e-19);
filter_16bit_response=filter_16bit_response(1:FS/2);

plot2d(filter_16bit_response, style=5, axesflag=1, rect=[0, -80, FS/2, 10]);

//-------------------------------------------------------------
// fixed point FIR filter simulation

rand('normal'); 
x=rand(1:FS);
y=zeros(1:FS);
buf=zeros(1:FIRLEN);

for i=1:FS
  buf=[x(i), buf(1:FIRLEN-1)];
  y(i)=(buf*coef16bit')/SCALE16bit;     // FIR filter 
end

y=y.*coef_gain;                         // adjust gain

xspec=20*log10(abs(fft(x, -1))+1e-19);
xspec=xspec(1:FS/2);

yspec=20*log10(abs(fft(y, -1))+1e-19);
yspec=yspec(1:FS/2);

//plot2d(xspec, style=2);
//plot2d(yspec, style=5);