Hilbert transfrorm via FIR filter procedure

function xH=hilbfir(x)
x=x(:);
% A length N=231 Remez Hilbert transformer with filtering procedure
% for long length data > (N-1)*3+1
% fp – Filter cutoff frequency (0.01>fp>0.5)
N=231-1;
fp=0.01;
h=firpm(N,[fp 0.5-fp]*2,[1 1],’Hilbert’);

xH0=filter(h,1,x);

% 90 degree – phase forward and reverse digital filtering
n=ceil((N+1)/2);
l=length(x); a=1;
b=h(:);
xH=filter(b,a,x);
xH=xH(n:l-n-1);
xb=x(l:-1:l-3*N);
xHb=-filter(b,a,xb);
xHb=xHb(3*N:-1:n);
xH=[xH; xHb(3*N+2-3*n:3*N-n+1)];

return
N=231-1;
fp=0.01;
h=firpm(N,[fp 0.5-fp]*2,[0.9998 0.9999],’Hilbert’);

% Filter spectral characteristics
[H,F]=freqz(h);
figure(1);
subplot(121)
plot(F/pi,abs(H),’.-‘);
axis([0 0.1 0.995 1.005])
xlabel(‘Frequency’)
ylabel(‘|FRF|’)
subplot(122)
plot(F/pi,abs(H),’.-‘);
axis([0.9 1 0.995 1.005])
xlabel(‘Frequency’);
ylabel(‘|FRF|’)

% Compare signal before and after filtering
figure(2)
x=sin(0.1*(1:1024))’;
xH=hilbfir(x);
plot(x,’-.’), grid on, hold on
plot(xH,’k-‘); hold off
legend(‘before’, ‘after filtering’)