There are two parts to your question. For the second part, assuming X(t) and Y(t') are independent, requires the cross spectral density to be zero. The cross spectral density is the Fourier transform of the cross correlation function. The cross correlation is the ensemble average of the time-shifted product of X(t) and Y(t'), and if these are independent zero-mean processes than the ensemble average is the product of the two means is zero, thus making the cross spectral density zero.
To find an expression for the cross spectral density start by showing the relationship of $X(t)$ to $S_X(f)$. A WSS random function $X(t)$ can be represented by a Fourier series over the interval $T$ $$X(t) = \sum_{m=-\infty}^{+\infty}A_me^{i2\pi mt/T + i\zeta _m}$$ where $A_m$ is the Fourier coefficient related to the spectrum $S_X(f_m)$, the frequency is given by $f_m = m/T$, and the random behavior of $X(t)$ is carried in the random phase $\zeta_m$. The phase $\zeta_m$ is a uniformly distributed random variable between $0$ and $2\pi$; a random value of $\zeta_m$ is assigned to each term in the Fourier series. An infinite number of ensembles of $X(t)$ can be created by choosing new values of the phases $\zeta_m$ from the distribution. The mean value of $X(t)$ is given by $A_0$ with $\zeta_0 = 0$, and for a zero-mean process $A_0 = 0$. (Note: most of the literature that I have seen erroneously treats $\eta_0$ as a random variable; this is wrong because this makes the mean value a member of a distribution instead of a constant as required by stationarity)
The autocorrelation $R(\tau)$ is given by the temporal average over $t$ or the ensemble averages over $\zeta_m$; here we use the ensemble averages $$R(\tau) = \langle X(t)X^*(t+\tau)\rangle = \sum_{m=-\infty}^{+\infty}\sum_{n=-\infty}^{+\infty}A_m A_n e^{i2\pi (m-n)t/T -i2\pi n\tau/T} \langle e^{i(\zeta _m-\zeta_n)}\rangle $$ The ensemble averages are zero everywhere except when $m = n$ so that the autocorrelation reduces to $$R(\tau) = \sum_{m=-\infty}^{+\infty}A_m^2 e^{-i2\pi m\tau/T} $$
An estimate of the $m^{th}$ component of the spectrum of $X(t)$ on the finite interval $T$ is given by taking the Fourier transform of the autocorrelation estimate$$S_m = \int_{-T/2}^{T/2}R(\tau) e^{i2\pi p\tau/T}d\tau = \sum_{m=-\infty}^{+\infty}A_m^2 \int_{-T/2}^{T/2} e^{-i2\pi (m-p)\tau/T}d\tau = A_m^2T $$ where the integral is zero except when $p = m$. If we think of $S_m$ as a sample of $S_X(f_m)$ then the Fourier coefficient $A_m$ can be related to the power spectral density$$A_m = \sqrt\frac{S_m}T = \sqrt\frac{S_X(f_m)}T$$ Substituting this value for $A_m$ into the autocorrelation function yields $$R(\tau) = \sum_{m=-\infty}^{+\infty}\frac{S_m}T e^{-i2\pi m\tau/T} $$ Define $2\pi /T = \Delta\omega$ and after substituting this into the expression for autocorrelation take the limit as $T$ goes to infinity $$R(\tau) = \lim_{T \to \infty}\sum_{m=-\infty}^{+\infty} S_m e^{-i\Delta\omega \tau} \Delta\omega = \frac{1}{2\pi}\int_{-\infty}^{+\infty}S(\omega) e^{-i\omega\tau}d\omega $$ where the Reimann sum has been converted to an integral. Assume all continuity and convergence requirements are met.
All of the foregoing was just to establish that $X(t)$ and $Y(t+\tau)$ can be written in terms of their spectra as $$X(t) = \sum_{m=-\infty}^{+\infty}\sqrt\frac{S_X(f_m)}T e^{i2\pi mt/T + i\zeta _m}$$ $$Y(t+\tau) = \sum_{n=-\infty}^{+\infty}\sqrt\frac{S_Y(f_n)}T e^{i2\pi n(t+\tau)/T + i\mu_n}$$ The cross correlation is the ensemble average of the product of $X(t)$ and the complex conjugate of $Y(t+\tau)$ $$R_{XY}(\tau) = \sum_{m=-\infty}^{+\infty}\sum_{n=-\infty}^{+\infty}\frac{\sqrt{S_X(f_m)S_Y(f_n)}}T e^{i2\pi (m-n)t/T -i2\pi n\tau/T} \langle e^{i(\zeta _m-\mu_n)}\rangle $$
To achieve a non-zero cross spectral density the ensemble average of the cross correlation function cannot be zero; this means that phases between the spectral components must be related. Otherwise, if the phases $\zeta_m$ and $\mu_n$ are independent the ensemble average will always be zero. So if the phases between components with the same frequency are assumed to be related by $\zeta_m - \mu_m = \delta_m$ and $\zeta_m - \mu_n = 0$ the cross correlation can be written as $$R_{XY}(\tau) = \sum_{m=-\infty}^{+\infty}\frac{\sqrt{S_X(f_m)S_Y(f_m)}}T e^{-i2\pi m\tau/T+i\delta_m} $$
Again, define $2\pi /T = \Delta\omega$ and after substituting this into the expression for autocorrelation take the limit as $T$ goes to infinity $$R_{XY}(\tau) = \lim_{T \to \infty}\sum_{m=-\infty}^{+\infty} \sqrt{S_X(f_m)S_Y(f_m)} e^{-i\Delta\omega \tau+i\delta(f_m)} \Delta\omega $$ $$ = \frac{1}{2\pi}\int_{-\infty}^{+\infty} \sqrt{S_X(\omega)S_Y(\omega)} e^{-i\omega\tau+i\delta(\omega)}d\omega $$ where the Reimann sum has been converted to an integral.
The cross power spectral density is just the Fourier transform of the cross correlation and it can be obtained by visual inspection of the above expression $$ S_{XY}(\omega) = \sqrt{S_X(\omega)S_Y(\omega)} e^{i\delta(\omega)} $$
A few observations on the above result. Although the question pertained to the spectra of the two processes, the more crucial parameter is the phase relationship $\delta(\omega)$. Also, consider relaxing the WSS requirement. One common mechanism that produces correlated processes is for one of them to be the time evolved version of the other; this only happens when WSS does not apply.
From the point of view of radiolocation, the general task can be divided into two parts: the consideration of jamming trends and the allocation of signal information. As far as I understand the essence of the issue, in meteorology the same information can be perceived as an interfering trend (with a short-term forecast), and as signaling information (with a long-term one).
Time trends (including polynomial one) can be used to equalize the signal level at the edges of the samples, which eliminates the use of weight windows that reduce the resolution of the DFT.
The same task can be performed by methods based on an autoregressive smoothing model, for which fast algorithms also exist.
Trends computed in the spectral region are used to normalize the density of the power spectrum. The alignment of the interference spectrum brings the consistent signal processing to the optimum.
The signal spectrum is chaotic and poorly predictable. Non-parametric processing methods are better suited for it — for example, median filtering in a sliding window.
Waiting for detalization in comments.
$\textbf{Comments to the comments}$
Note that the AR model is significantly marked from polynomial.
The polynomial model is not fundamentally cyclical, while the AR model of small order ideally approximates a harmonic signal.
Fundamentally different structure of matrices. The matrix structure in polynomial regression problems leads to the Vandermonde determinant and is characterized by computational instability, which is manifested in models of order 10 and higher.
The matrix for the AR model has a special Toeplitz structure, and its solution is built on recurrent formulas, which make it possible to consistently increase the order of the model. The considered methods to a significant extent rely on the autocorrelation function, which is not typical of polynomial models.
- The signal maxima in the framework of the polynomial model are also calculated by well-known algorithms.
The signal maxima within the AR model correspond to the poles of the signal spectrum.
In my opinion, these differences are enough to abandon the idea of linking these models. If the search continues, they should take into account these differences.
Best Answer
The last equality in the first line of your equation is wrong. Note that you have
$$u(t)u(t+\tau)=\begin{cases}u(t),&\tau>0\\u(t+\tau),&\tau<0\end{cases}$$
Now compute the integral for both cases ($\tau>0$ and $\tau<0$), and you'll get the desired result.