function [u, b] = zp_filter(x, b, varargin)
% Zero-phase filter
% u = zp_filter(x, 'hann')
% u = zp_filter(x, 'boxcar', N)
% u = zp_filter(x, 'trapez', N, N2M)
% u = zp_filter(x, b)
%
% N: Filter length
% N2M: cut-off period in units of delta t
% length(b) must be odd
% (first and last filtered data points filled with NaN)
%
% April 2001, Nils Olsen, DSRI
%

[N_row, N_col] = size(x);
u = NaN*ones(size(x));
% x = x(:); % row vector

if ischar(b)
    if strcmp(b, 'boxcar')
        if nargin < 3 error('boxcar requires filter length as 3rd input'); end
        N = varargin{1};
        b = ones(N, 1)/N;
    elseif strcmp(b, 'hann')
        b = [.25 .5 .25];
    elseif strcmp(b, 'trapez')
        if nargin < 4 error('trapez requires 4 input elements'); end
        NFILT = varargin{1};
        if mod((NFILT-1)/2, 1) ~= 0; error('N must be odd'); end
        N2M = varargin{2}/2;
        NFILT = (NFILT-1)/2;
        E = 1.0E-6; JM = 6; FAC = [0.77,0.95,0.97]'; k=0;
        N=abs(NFILT);
        T=N2M;
        IQ=floor(N/T+0.5);
        A=pi/T;
        b=pi/N;
        if IQ <= 3; b=b*FAC(IQ); end
        W0=1/T;
        %     NULLSTELLENSUCHE VON Y(PI*BETA)
        while 1 > 0
            fx=W0-1;
            fxp=0;
            for I=1:N
                U=A*I;
                V=b*I;
                D=sin(U)/U;
                if (V < 1.0E-4)
                    V2=V^2;
                    P=1.-V2/6;
                    Q=-V*(1.-0.1*V2)/3;
                else
                    P=sin(V)/V;
                    Q=(cos(V)-P)/V;
                end
                fx=fx+2*W0*D*P;
                fxp=fxp+2.*W0*I*D*Q;
            end
            if (abs(fx) < E); break; end
            k=k+1;
            if (k > 30)
                error(num2str(k));
                return
            end
            b=b-fx/fxp;
        end
        % Filter weights w
        w = [];
        for I=1:N
            U=A*I;
            V=b*I;
            w(I)=W0*sin(U)*sin(V)/U/V;
        end
        b = [fliplr(w) W0 w]';
    else
        error('unknown filter type');
        return
    end
end
N = (length(b)-1)/2;
if mod(N, 1) ~= 0; warning('b must be of odd length'); end

for k = 1:N_col
    y = filter(b, 1, x(:,k));
    % u = [y(N+1:end); NaN*ones(N, 1)];
    u(:,k) = [NaN*ones(N, 1); y(2*N+1:end); NaN*ones(N, 1)];
end

end
