I have limited experience with optics and imaging, so was unfamiliar with the term Modulation Transfer Function. This is what I have learned from this article on Wikipedia combined with my experience with Fourier Transforms which is hopefully sufficient to provided a helpful answer. This is an interesting application demonstrating navigating all the Fourier Transform options and which to apply and how, especially when working in other domains other than time and frequency (here space and frequency).
The more generalized Optical Transfer Function (OTF) is the Fourier Transform of the Point Spread Function (PSF) for an optical system. The PSF is the system's impulse response. The OTF will be complex with magnitude and phase, consistent with any causal linear time invariant system. Thus the complex result can be represented in magnitude and phase form, specifically $A(f)e^{j\phi(f)}$, with $A(f)$ and $\phi(f)$ both real functions where $A(f)$ is the Magnitude Transfer Function and $\phi(f)$ is the Phase Transfer Function.
The frequency units are spatial frequency in cycles/mm, and the resulting transform can be in multiple dimensions such as 1D, 2D, or 3D depending on the dimension of the Point Spread Function. Here the OP has given a 1D point spread function as 3 impulses. For purposes of providing units in this answer I will assume each of those impulses are separated by 1 mm. This would physically be realized by a point source going through a lens that converges to three points each separated by 1 mm on the image plane due to lens distortions.
Which Transform to use??
More generally the Point Spread Function would be a continuous function in space (such as a 2D ellipse for example, consistent with "point spread"), so for that the Continuous Space Fourier Transform would be used (CSFT). Given the Point Spread Function in this case by the OP is discrete, we can use here the Discrete Space Fourier Transform (DSFT). I provide the equations for both forms below along with that for the Discrete Fourier Transform (DFT), since samples of the DSFT can be computed from the DFT. This is an important point since the Fast Fourier Transform (FFT) is a very efficient algorithm for computing the DFT. The CSFT and DSFT are continuous functions in frequency, while the DFT is discrete in frequency. By zero padding the spatial domain waveform prior to computing the DFT, we can interpolate more samples of the desired DSFT thus providing an efficient algorithmic solution for general applications.
CSFT:
$$Y(f) = \int_{x=-\infty}^{\infty}y(x)e^{-j2\pi f x}dx \tag{1}\label{1}$$
Where
$f$ is frequency in cycles/mm as a continuous function
$x$ is distance in mm as a continuous function
$y(x)$ is a spatial domain function as a function of x
$Y(f)$ is a frequency domain function as a function of f
DSFT:
$$Y(f) =\sum_{n=-\infty}^{\infty} y[n]e^{-j 2 \pi f n} \tag{2}\label{2}$$
Where
$f$ is frequency in cycles/mm as a continuous function
$n$ is distance in mm as a discrete function
DFT:
$$Y(k) =\sum_{n=0}^{N-1} y[n]e^{-j 2 \pi n k/N} \tag{3}\label{3}$$
Where
$k$ is a frequency index from $0$ to $N-1$
$n$ is a spatial index from $0$ to $N-1$
$N$ is the total number of samples in space and frequency
Solution
Given that the Point Spread Function (PSF) is discrete, as $y[n] = [1,1,1]$ we would use the DSFT (Equation \ref{2}).
We can reasonably assume the spatial center is in the middle of the PSF (the PSF as given extends from $n=-1$ to $n=1$), the DSFT result would be:
$$Y(f) =\sum_{n=-\infty}^{\infty} y[n]e^{-j2 \pi f n} = e^{-j2 \pi f} + 1 + e^{j2 \pi f}$$
Using this relationship from our friend Euler:
$$2\cos(\omega) = e^{j \omega} + e^{-j \omega}$$
The above simplifies to the solution provided by the OP as:
$$Y(f) = 1 + 2\cos(2\pi f)$$
This is the OTF, but since the result is completely real and always positive, it happens to be the MTF as well, but that wouldn't necessarily always be the case. Whenever the PSF is symmetric, the result will be real, but if it is signed, the MTF would be the magnitude of that result.
I prefer using the units $f$ in the result and not $k$ since the resulting response in frequency is not discrete, and $k$ is very commonly used as the discrete frequency index in the DFT. We note if the units in the Point Spread Function (PSF) are mm, the units of frequency $f$ in the solution is cycles/mm.
Also to comment if the PSF was not centered on $n=0$ as I applied here, this would have no effect on the MTF since any such shift would effect the phase only. For example if the PSF started at the $n=0$, the resulting OTF would be:
$$Y(f) =\sum_{n=-\infty}^{\infty} y[n]e^{-j2 \pi f n} = 1 + e^{-j2 \pi f} + + e^{-j4 \pi f}$$
$$ = e^{-j2 \pi f}(e^{-j2 \pi f} + 1 + e^{j2 \pi f}) =e^{-j2 \pi f}( 1 + 2\cos(2\pi f)) $$
Thus giving us the same MTF and we see the PTF would be $\phi(f) = 2 \pi f$.
