週末やったFFTによるフィルタの実験プログラム

とりあえず、プログラムだけ。

function yn=fftfilter(xn, blocksize)
  index = 1;
  zeropad = [1:blocksize] - [1:blocksize];
  w = window( 'hn', blocksize )
  [h, wfm, wfr] = wfir( 'lp', blocksize, [0.1 0.2], 'hn', [0 0]);
  hspectrum = fft([h zeropad]);
  eoblock = index+blocksize-1
  eodoubleblock = index+blocksize*2-1
  
  yn = xn - xn; // make yn as same figure with xn, but zero padding
  while ( eodoubleblock <= size( xn,2 ) ) do
    block = [xn( index:eoblock ) .* w zeropad];  // get block
    spectrum = fft( block ) .* hspectrum;       // convolution
    block = ifft(spectrum);
    yn( index:eodoubleblock ) = yn( index:eodoubleblock ) + block;
    index = index + blocksize/2; // slide block window half
    eoblock = eoblock + blocksize/2;
    eodoubleblock = eodoubleblock + blocksize/2;
  end
endfunction