%% Octave/Matlab implementation of Michael Mann's lowpass "filter" from
%% http://www.meteo.psu.edu/holocene/public_html/Mann/tools/MTM-SVD/mtm-svd-recon.f
%% Relevant parts of the code reproduced for the reference (lines starting with '%c')
%%
%% Jean S/17.1.2013
%%
%% usage: y = lowmann(x, Wc [, window])
%% 	lowpass filter with cutoff pi*Wc radians
%%      window type:
%%       1) rectangular
%%       2) tapered rectangular
%%       3) triangular
%%       4) Hann(ing)
%%       5) Hamming (default)
%%       6) Blackman
function [y,b] = lowmann(x, W, varargin)
  
  if nargin < 2 || nargin > 3
    usage('y = lowmann(x, W [, window])');
  end
  
  %% Assign a window
  if length(varargin)==1
    iwndo = varargin{1};
  else 
    iwndo = 5; %% Hamming window (default in Mann's code)
  end

  x = x(:);                   % make it a column vector  
  lx = size(x,1);
  hx = fix(lx/2); 
  
  %% Impulse response (of length lx !!!) by windowing method
  b=spfird_low(lx-1,W,iwndo);
  %%
  %c      call spfird(N-1,iband,f1,f2,iwndo,b,ierror)
  
  %% Now convolve the input signal with the impulse response
  %% and cut off half of the filter length from the both ends.
  y = conv(x, b); 
  y = y(hx+1:hx+lx);
  %%
  %c      do i=0,2*N
  %c         c(i)=zero
  %c      end do
  %c      do i=0,N-1
  %c         do j=0,N-1
  %c            c(i+j)=c(i+j)+x(i)*b(j)
  %c        end do
  %c      end do
  %c      do i=0,N-1
  %c         yy(i+1)=c(i+N/2)
  %c      end do
  %c      return
end

function b = spfird_low(n,W,iwndo)
  %% 
  %% Design a lowpass FIR filter of order n (length=n+1) by windowing method
  %% 
  %% Mann is using code from 
  %% Samuel D. Stearns and Ruth A. David: Signal Processing Algorithms in Fortran and C, Prentice Hall, 1992.
  %%
  %% Original FORTRAN code from SBLIB.FOR -library available at
  %% http://www.ece.unm.edu/faculty/stearns/downloads/Website%20Fortran/Signal_Processing_Algorithms/
  %% 
  %% Lowpass Octave implementation of spfird follows.
  hn = fix(n/2); 
  evenorder = hn==n/2;
  hn=hn-evenorder;
  wc1=pi*W;

  b=ones(n+1,1)*NaN;
  s=[0:hn]'-n/2;
  window=spwndo(iwndo,n+1); 
  b(1:hn+1)=sin(wc1*s)./(pi*s).*window(1:hn+1);
  b(end-hn:end)=b(hn+1:-1:1);
  if evenorder, b(hn+2)=wc1/pi; end % if odd length add the middle coefficient
  %%
  %c       dly=l/2.
  %c       lim=int(l/2)
  %c       mid=0
  %c       if(dly.eq.lim) then
  %c         lim=lim-1
  %c         mid=1
  %c       endif
  %c       wc1=2.*pi*fln
  %c       if(iband.ge.3) wc2=2.*pi*fhn
  %c       go to (5,10,15,20) iband
  %c c-lowpass design
  %c     5 continue
  %c       do 6 i=0,lim
  %c         s=i-dly
  %c         b(i)=((sin(wc1*s))/(pi*s))*spwndo(iwndo,l+1,i)
  %c         b(l-i)=b(i)
  %c     6 continue
  %c       if(mid.eq.1) b(l/2)=wc1/pi
  %c       return
end

function window = spwndo(itype,lw)
  %c c-itype=1(rectangular), 2(tapered rectangular), 3(triangular),
  %c c-      4(hanning), 5(hamming), or 6(blackman).
  if itype==1, % Rectangular
      window = ones(lw,1);
  elseif itype==2, % Tapered rectangular
      l=fix((lw-2)/10);
      window=ones(lw,1);
      window(1:l)=0.5*(1-cos([0:l-1]*pi/(l+1)));
      window(end-l+1:end)=flipud(window(1:l));
      %c  2 l=(n-2)/10
      %c    if(k.le.l) spwndo=0.5*(1.0-cos(k*pi/(l+1)))
      %c    if(k.gt.n-l-2) spwndo=0.5*(1.0-cos((n-k-1)*pi/(l+1)))
  elseif itype==3, % Triangular
      window = 1-abs(1-2*[0:lw-1]/(lw-1));
      %c  3 spwndo=1.0-abs(1.0-2*k/(n-1.0))
  elseif itype==4, % Hann(ing)
      window = 0.5*(1-cos(2*pi*[0:(lw-1)]'/(lw-1)));
      %c  4 spwndo=0.5*(1.0-cos(2*k*pi/(n-1)))
  elseif itype==5, % Hamming
      window = 0.54-0.46*cos(2*pi*[0:(lw-1)]'/(lw-1));
      %c  5 spwndo=0.54-0.46*cos(2*k*pi/(n-1))
  elseif itype==6, % Blackman
      window = 0.42-0.5*cos(2*pi*[0:(lw-1)]'/(lw-1))+0.08*cos(4*pi*[0:(lw-1)]'/(lw-1));
      %c  6 spwndo=0.42-0.5*cos(2*k*pi/(n-1))+0.08*cos(4*k*pi/(n-1))
  else
     error('Window type is not defined.');
  end
end