Get the frequency response curve from FIR filter coefficients + sampling rate
As suggested by a comment, scipy.signal.freqz gives the solution. Don't forget to use this to do the rad/sample to Hz conversion.

Here is a ready-to-use code, posted here for future reference:

import numpy as np
from scipy import signal
import matplotlib.pyplot as plt
sr = 44100
w, h = signal.freqz(b=[1, -1], a=1)
x = w * sr * 1.0 / (2 * np.pi)
y = 20 * np.log10(abs(h))
plt.figure(figsize=(10,5))
plt.semilogx(x, y)
plt.ylabel('Amplitude [dB]')
plt.xlabel('Frequency [Hz]')
plt.title('Frequency response')
plt.grid(which='both', linestyle='-', color='grey')
plt.xticks([20, 50, 100, 200, 500, 1000, 2000, 5000, 10000, 20000], ["20", "50", "100", "200", "500", "1K", "2K", "5K", "10K", "20K"])
plt.show()


With logarithmic x-axis:



Another example with a comb filter (COEF = [1] + [0] * 200 + [-1], i.e. y[n] = x[n] - x[n-201]):