Approximate the magnitude response of an analog filter with a digital filter, starting from its gain expression I have an analog filter with its frequency response curve in dB described by the following expression:
$$
N_{dB}=20log_{10}\omega t_1 \sqrt{\frac{1+(\omega t_2)^2}{1+(\omega t_1)^2}}
$$
This expression is derived from the series connection of two lowpass filters each associated with the following RC circuit:

where, for each circuit, the time constant $t_i=RC$ and of course $\omega = 2\pi f$, where $f$ is the frequency (this is actually the equalization curve for magnetic tape recording/playback, see Annex B, page 14 of this document).

I would like to obtain an approximation of this frequency response using a digital filter. I don't know if there is a method to exploit our knowledge of the analog frequency response or if I should design the filter myself from scratch.
The end goal is to obtain the impulse response and save it as a a .wav file (I know how to do this last part). I just took a basic DSP course at my Uni but we didn't work with analog filters so I am a little bit lost.
Did you try to plot this curve in matlab.
I am not sure whether it is any standard high,low or band pass filter.
So, to calculate the impulse response there might be two ways.
1. take the inversefft of the NdB and then do the convolution between your impulse input and the inversefft signal of filter usinc conv(x,y) command.
2. Or convert the impulse input into frequency domain by fft transform in matlab and then just multiply the gain and fft of input impulse and then convert them in time domain.
Can you do it using above method? 
Do you want a matlab code for this method?