// FIR filter desing by inverse DFT(FFT) and transition samples

clear;

FS=10000;              // sampling frequency
FIRLEN=255;            // FIRLEN should be odd number
speclen=FIRLEN+1;      // filter specification length is even number

// set filter specification
freq_response=[zeros(1:20), 0.9, ones(1:20), 0.7, 0.1, zeros(1:30), 0.1, 0.7, ones(1:20), 0.7, 0.3, zeros(1:speclen)];
//freq_response=[zeros(1:20), 0.0, ones(1:20), 0.0, 0.0, zeros(1:30), 0.0, 0.0, ones(1:20), 0.0, 0.0, zeros(1:speclen)];
freq_response=freq_response(1:(speclen/2+1));
freq_response=[freq_response, mtlb_fliplr(freq_response(2:(speclen/2)))];

// plot filter specification
scf(0);
clf;
plot2d(freq_response(1:speclen/2), axesflag=1, rect=[1, 0, speclen/2, 2.0], style=5);
plot2d(freq_response(1:speclen/2), axesflag=1, rect=[1, 0, speclen/2, 2.0], style=-1);

h=real(fft(freq_response, 1));                       // inverse FFT
h=[h((speclen/2+2):speclen), h((1:speclen/2))];      // sort

// plot filter coefficient
scf(1);
clf;
plot(h);

hspec=20*log10(abs(fft([h, zeros(1:(FS-FIRLEN))], -1))); 
hspec=hspec(1:(FS/2));

// plot filter frequency response
scf(2);
clf;
plot2d(hspec, axesflag=1, rect=[0, -80, FS/2, 10]);

// print filter coefficient
for i=1:FIRLEN
  mprintf("%3i %15.12f\n", i-1, h(i));
end