For finite differences, as suggested by Andy, you can use MATLAB/Octave's freqz, but with some tricks to get the phase frequency response right. I'm using Wikipedia's definitions of higher-order finite differences. In Octave:
graphics_toolkit("gnuplot");
n = 2; #finite differences derivative order
s = 0; #-1 for forward, 0 for central, 1 for backward
b = [1];
a = [1];
for k=1:n,
b = conv(b, [1, 0, -1]);
a = conv(a, horzcat(zeros(1, 1-s), [1]));
end;
w = [(0:512)/512*pi-pi/2];
h = freqz(b, a, w);
freqz_plot(w*2/pi, h, true);
Feel free to change n and s. b will contain the filter coefficients, diluted by zeros to trick freqz to calculate frequency responses when taps may be at half-integer times, and a will have a 1 to mark the tap at time zero in the impulse response. For the above:
b =
1 0 -2 0 1
a =
0 0 1
These are coefficients of the numerator and denominator polynomials of $z^{-1}$ of the transfer function, respectively. For forward and central finite differences the filter is acausal. The difference to the causal backward finite differences filter is that the impulse response has been shifted backwards $m$ samples by multiplying the transfer function by $z^m$ with positive $m.$ This introduces negative powers of $z^{-1}$ such as $(z^{-1})^m$ in the numerator. This can be made valid input to freqz by multiplying both the numerator and denominator by $z^{-m}$ to ensure both are polynomials of $z^{-1}$. Effectively, the script above does that.
Figure 1. freqz_plot output for central 2nd derivative finite differences.
The magnitude frequency response will be equal for all: forward, backward, central.
