The trick is to evaluate the retarded Green's function of the d'Alembert operator. So let
$$\Box G_n(x)=\delta^{(1+n)}(x), \quad G_n(x)=0 \quad \text{for} \quad t<0,$$
where ##x=(t,\vec{x})## is the (1+d)-dimensional spacetime vector of the Minkowski space with metric ##\eta_{\mu \nu}=\mathrm{diag}(1,-1,-1,\ldots,-1)##. We set the phase speed of the wave, ##c=1##. The d'Alembert operator is ##\eta^{\mu \nu} \partial_{\mu} \partial_{\nu} = \partial_t^2-\Delta##, where ##\Delta=\vec{\nabla}^2## is the ##n##-dimensional Laplace operator.
It's most easy to work with the socalled Mills representation, i.e.,
$$G_n(x)=\int_{\mathbb{R}^n} \frac{\mathrm{d}^n k}{(2 \pi)^n} \exp(\mathrm{i} \vec{x} \cdot \vec{x}) \tilde{G}_n(t,\vec{k}).$$
Plugging this "Fourier ansatz" into the equation of motion after some calculation one gets the harmonic-oscillator equation
$$(\partial_t^2+\vec{k}^2) \tilde{G}_n(t,\vec{k})=\delta(t).$$
Solving with the condition for the retarded propagator leading to the initial-value problem ##\tilde{G}_n(0^+,\vec{k})=0##, ##\partial_t \tilde{G}_n(0^+,\vec{k})=1## leads to
$$\tilde{G}_n(t,\vec{x}) = \frac{\Theta(t)}{||vec{k}|} \sin(|\vec{k}|t).$$
Now let's consider the cases ##n=3,2,1## (in that order). For ##n=3## we can evaluate the Fourier integral, using polar coordinates for ##\vec{k}##. Since ##G_n(t,\vec{x})## is rotation invariant, it depends only on ##r=|\vec{x}|##, i.e., we can choose the polar axis in direction of ##\vec{x}## without loss of generality. Then we get (with ##u=\cos \vartheta##)
$$G_3(t,\vec{x}) = \frac{\Theta(t)}{(2 \pi)^3} \int_0^{\infty} \mathrm{d} K \int_0^{2 \pi} \mathrm{d} \varphi \int_{-1}^1 \mathrm{d}u K \sin(K t) \exp(\mathrm{i} K r u).$$
Integration over ##\varphi## and ##u## yields
$$G_3(t,\vec{x}) = -\frac{\Theta(t)}{4 \pi^2 r} \int_{\mathbb{R}} \mathrm{d} K \sin(K r) \sin(K t),$$
where we have used the fact that the integrand is an even function in ##K## to extend the integration region for ##K## to the entire real axis. Resolving the sine function in terms of exp functions one gets ##\delta## distributions. Together with the unit-step function and using ##r>0## one finally gets the well-known result
$$G_3(t,\vec{x})=\frac{\Theta(t)}{4 \pi r} \delta(t-r).$$
Thus in ##3## spatial dimensions the Huygens principle is valid. This is seen by using the Green's function for solving the problem of diffraction for scalar waves on obstacles in the famous Kirchhoff approximate solution of the problem, which at the end boils down to the diffracted wave field given by the superposition of the spherical waves originating from each point in the openings of the diffractor (like a slit, double-slit, grating etc.).
Now from the Fourier integral of the Mills representation, it is easy to see that
$$G_{n-1}(t,\vec{x})=\int_{\mathbb{R}} \mathrm{d} x_n G_{n}(t,\vec{x}).$$
Using this for ##n=3##, we find
$$G_{2}(t,\vec{x})=\frac{\Theta(t)}{4 \pi} \int_{\mathbb{R}} \mathrm{d} x_3 \frac{\delta(t-r)}{r} = \frac{\Theta(t)}{2 \pi} \int_{-\infty}^{\infty} \mathrm{d} x_3 \delta (t^2-r^2).$$
Setting ##\rho^2=\vec{x}^2-x_3^2=x_1^2+x_2^2##, we get
$$G_{2}(t,\vec{x})=\frac{\Theta(t)}{2 \pi} \Theta(t^2-\rho^2) \int_{-\infty}^{\infty} \frac{1}{2|x_3|} [\delta(x_3-\sqrt{t^2-\rho^2}-\delta(x_3+\sqrt{t^2-\rho^2})] = \frac{\Theta(t-\rho)}{4 \pi \sqrt{t^2-\rho^2}}.$$
As you see, that's no ##\delta## distribution anymore. Of course, there's the singularity along the future light cone ##\rho=t## as it must be to get a Green's function, but it's not a ##\delta## distribution on the light cone. Thus the diffracted wave pattern is not simply a superposition of spherical waves going out from each point in the openings if the diffraction grating but it's somewhat "smeared around the light cone".
For ##n=1## one better uses the Mills-Fourier integral directly,
$$G_1(t,x)=\Theta(t) \int_{\mathbb{R}} \frac{\mathrm{d} k}{2 \pi} \frac{\sin(x t)}{k} \exp(\mathrm{i} k x).$$
To evaluate this integral, we note that ##\sin(k t)/k## for ##t>0## is the Fourier transform of
$$f(x)=\frac{1}{2} \Theta(-t<x<t)=\frac{1}{2} \Theta(t-|x|),$$
because
$$\int_{-t}^t \mathrm{d} x \exp(-\mathrm{i} k x)=\frac{2 \sin(k x)}{k}.$$
Thus we finally get
$$G_1(x)=\frac{1}{2} \Theta(t-|x|).$$
So in 1 dimension Huygens's principle doesn't hold either, but that's a special case.
If I remember right, the general case for odd dimensions ##n \in \{3,5,7,\ldots \}## Huygens's principle is valid. A proof is in
S. Hassani, Mathematical Physics, Springer Verlag, Cham, Heidelberg, New York,
Dordrecht, London, 2 ed. (2013).
https://dx.doi.org/10.1007/978-3-319-01195-0