Polar Form of the 2d DFT

Question

I know that there is a complex and a polar formulation of the Fourier series. The polar form can be transformed to the real valued Fourier series for the 1D case: $$f(x)=a_{0}+\sum_{k}^{\infty} a_{k}\cos(k\omega_{0}x)+b_{k} \sin(k\omega_{0}x)= d_{0}+\sum_{k=1}^{\infty}d_{k}cos(k\omega_{0}x-\phi_{k})$$ with $\omega_{0} := 2\pi /T$ and $$d_{0}=a_{0}, d_k:=\sqrt{a^2_k +b^2_k}\quad\text{as magnitude and the phase}\quad\phi_{k}:=\arctan{\frac{b_{k}}{a_{k}}}$$ , which can be calculated as follows for a dft: $$d_k \cos\phi_k = \sqrt{2/N}\sum^{N}_{j=1}z_j \cos(2\pi k(j-1)/N)$$ $$d_k \sin\phi_k = \sqrt{2/N}\sum^{N}_{j=1}z_j \sin(2\pi k(j-1)/N)$$ It's essentially taking the real and the imaginary parts separately and using trigonometric relationships to calculate the coefficients. Now for the 2D case, I have seen this approach(in a paper): $$d_k \cos\phi_k = 1/N\sum^{N}_{j=1}\left(x_j\cos\left(\frac{2\pi (k+1)(j-1)}{N}\right)-y_j\sin\left(\frac{2\pi (k+1)(j-1)}{N}\right)\right)$$ $$d_k \sin\phi_k = 1/N\sum^{N}_{j=1}\left(x_j\sin\left(\frac{2\pi (k+1)(j-1)}{N}\right)+y_j\cos\left(\frac{2\pi (k+1)(j-1)}{N}\right)\right)$$ I just do not understand, where it is coming from. I also didn't encounter any book on DFT, which explains this. There were rare ones which mentioned at least the polar formulation. I know how the complex case works. Does someone know how the above formulas for the 2d and 1d case are derived or any literature about the polar form?
And maybe how the treatment with just sin and cos work? I am not afraid of complex numbers, but I feel like I want to understand this way of doing it too.
I would really appreciate any information.