Differentiator via FIR filter procedure

function xd=diffir(x,Fs)

% A length N=231 Remez differentiator with filtering procedure
% for long length data > (N-1)*3+1
% Fs – Sampling frequency
N=231-1;

h=firpm(N,[0 0.96],[0 0.96*pi*Fs],’d’);

n=ceil((N+1)/2);
l=length(x); a=1;
b=h(:); x=x(:);
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);
xd=[xH; xHb(3*N+2-3*n:3*N-n+1)];Frequency response function

%return
% Filter spectral characteristics
N=231-1;
Fs=1;
h=firpm(N,[0 0.96],[0 0.96*pi*Fs],’d’);
[H,F]=freqz(h);
figure(1);
subplot(111)
plot(F/pi,abs(H),’.-‘);
xlabel(‘Frequency’)Initial signal and veliocity
ylabel(‘|FRF|’)

% Compare signal before and after filtering
figure(2)
x=sin(1*(1:1024))’;
xd=diffir(x,Fs);
plot(x,’-.’), grid on, hold on
plot(xd,’k-‘); hold off
axis([500 600 -1.5 1.5]);
legend(‘before’, ‘after differentiating’)