I have $2^{nd}$ order bessel filter transfer function. I want to model this analog bessel filter in FPGA. For this I need its digital or discrete form. I am trying to get digital filter's coefficients in MATLAB using $c2d$ function. But I see there are many options for this function, e.g. by-default it uses zero-order hold. When $\textbf{tustin}$ is specified it uses bilinear-transformation for discretization.
I don't know when to use which option while using $c2d$ ?
Before I detail the mapping approaches further let me comment since the OP mentions using a Bessel filter, that a Bessel filter is an analog filter that has a maximally-flat group delay and very nice transient characteristics (no overshoot/ringing), and is often selected primarily due to these features. The maximally flat-group delay doesn't carry over well when mapped from s to z using the common mapping approaches offered. Instead of mapping at all, I recommend implementing the all-pole Thiran filter, which is the digital equivalent to the Bessel filter that can be implemented directly. (The Thiran is also commonly used for fractional delays as an all-pass filter when implemented with poles and zeros). The Thiran filter has the maximally-flat group delay characteristic and nice transient response and is the popular choice for implementation when the digital equivalent of a Bessel filter is needed. For further implementation details of the Thiran, see ConcernedCitizen's answer to DSP.SE #82411. Also related to this, see the 6th paragraph down in this excellent response by MattL that has influenced me on this, where he states that the Bessel filter is the only analog filter that he would consider obsolete as a prototype for digital design, and references the originating 1971 paper by Thiran on digital maximally flat filter implementation.
Regarding general mapping from s to z; there is no single transformation that can provide an exact match both in the time domain and in the frequency domain. The two most popular methods when it is desired to match either the time response or match the frequency response are the “Method of Impulse Invariance” (using 'impulse'), and the “Bilinear Transform” (using 'tustin'). If the sampling rate is sufficiently higher than any of the finite poles and zeros in the filter, and there isn't a motivation to reduce it, then either approach will provide an excellent match in both time and frequency. However typically we are motivated for consideration of resources to run the filter at the lowest sampling rate possible, and it would be in this condition where careful selection of the mapping used is important.
The Method of Impulse Invariance is used when an exact time domain match is needed (the resulting impulse response, formally referred to as the unit sample response, of the resulting digital system will be exact samples of the impulse response of the continuous time system). The result is an exact match in the time domain, but not in the frequency domain and thus would not maintain the maximally flat group-delay of a Bessel filter.
The Bilinear Transform is used when an exact match in the frequency domain is needed, although the response will be “warped” in that the response in the frequency range from DC to infinity from the continuous-time case will be logarithmically mapped into the DC to Nyquist range. The tool supports prewarping to get an exact result for any particular frequency, using chosen where critical poles exist (but an exact match at all frequencies is not possible, thus the method is ill suited for maintaining the maximally flat group-delay of a Bessel filter).
Those are the primary considerations and the Bilinear Transform is often my go-to unless I do need an exact match of the time domain impulse response. There are other advantages and reasons to use the other methods (as well as further details on the two I mentioned) detailed in these other DSP.SE posts, as well as the posts listed in the "Related" column listed to the right on DSP.SE:
https://dsp.stackexchange.com/a/90065/21048
https://dsp.stackexchange.com/a/115/21048
https://dsp.stackexchange.com/a/62284/21048
https://dsp.stackexchange.com/a/54033/21048
https://dsp.stackexchange.com/a/113/21048
https://dsp.stackexchange.com/a/62565/21048
Here is the complete summary of all supported methods from Matlab’s documentation:
'zoh' — Zero-order hold (default). Assumes the control inputs are piecewise constant over the sample time Ts.
'foh' — Triangle approximation (modified first-order hold). Assumes the control inputs are piecewise linear over the sample time Ts.
'impulse' — Impulse invariant discretization
'tustin' — Bilinear (Tustin) method. To specify this method with frequency prewarping (formerly known as the 'prewarp' method), use the PrewarpFrequency option of c2dOptions.
'matched' — Zero-pole matching method
'least-squares' — Least-squares method
'damped' — Damped Tustin approximation based on the TRBDF2 formula for sparse models only.
