# Is it possible to compute the PSD of s(t)^2 given the PSD of s(t)?

Given a well-behaved signal $$s(t)$$ given by,

$$s(t) = \sum_{i}^{N} \cos( \omega_it + \phi_i)$$

with a large $$N$$ and a well-known power spectral density (PSD), $$\tilde{S_1}(\omega)$$. Is it possible to compute the PSD for $$s(t)^2$$, $$\tilde{S_2}(\omega)$$?

I am pretty sure (checked numerically and it makes sense) there is a relationship between $$\tilde{S_2}(f)$$ and the convolutions:

$$\tilde{S_2}(\omega) \propto \bigg( \int_{-\infty}^{+\infty} \tilde{S_1}(\nu)\tilde{S_1}(\omega+\nu)d\nu \bigg) + \bigg( \int_{-\infty}^{+\infty} \tilde{S_1}(\nu)\tilde{S_1}(\omega-\nu)d\nu \bigg)$$

What I've tried so far:

I started by expressing $$s(t)^2$$ as a sum:

$$s(t)^2 = \sum \sum_{i,j}^{N,N} \cos( \omega_it + \phi_i) \cos( \omega_jt + \phi_j)$$

Then invoking the identity $$\cos(a) \cos(b) = \frac{1}{2}(\cos(a-b)+ \cos(a+b))$$ we obtain:

$$s(t)^2 = \frac{1}{2} \sum \sum_{i,j}^{N,N} \bigg( \cos( (\omega_i - \omega_j )t + (\phi_i - \phi_j))+ \cos( (\omega_i + \omega_j)t + (\phi_j-\phi_i)) \bigg)$$

Applying the Fourier transform we obtain:\begin{align*} \mathcal{F}[s^2(t)](\omega) = \frac{1}{2} \sqrt{\frac{\pi}{2}} \bigg( \sum_{i} \sum_{j} \bigg[ e^{-i(\phi_i + \phi_j)}& \delta(\omega_i + \omega_j - \omega) +\\ e^{i(\phi_i - \phi_j)}& \delta(\omega_i - \omega_j + \omega) +\\ e^{-i(\phi_i - \phi_j)}& \delta(-\omega_i + \omega_j + \omega) +\\ e^{i(\phi_i + \phi_j)}& \delta(\omega_i + \omega_j + \omega) \bigg] \bigg) \end{align*}

and then squaring:

\begin{align*} \mathcal{F}[s^2(t)](\omega)^2 = \frac{1}{4} \frac{\pi}{2} \times ...\\ \bigg( \sum_{i} \sum_{j} \bigg[ e^{-i(\phi_i + \phi_j)} \delta(\omega_i + \omega_j - \omega) +& e^{i(\phi_i - \phi_j)} \delta(\omega_i - \omega_j + \omega) +\\ e^{-i(\phi_i - \phi_j)} \delta(-\omega_i + \omega_j + \omega) +& e^{i(\phi_i + \phi_j)} \delta(\omega_i + \omega_j + \omega) \bigg] \bigg)\times ...\\ \bigg( \sum_{k} \sum_{l} \bigg[ e^{-i(\phi_k + \phi_l)} \delta(\omega_k + \omega_l - \omega) +& e^{i(\phi_k - \phi_l)} \delta(\omega_k - \omega_l + \omega) +\\ e^{-i(\phi_k - \phi_l)} \delta(-\omega_k + \omega_l + \omega) +& e^{i(\phi_k + \phi_l)} \delta(\omega_k + \omega_l + \omega) \bigg] \bigg) \end{align*}

Only square terms survive (to be demonstrated). A factor of two comes out because there are two relevant cases $$k = i, l = j$$ and $$k = j, l = i$$ (not demonstrated, but I intuit it should be so), and so:

\begin{align*} \mathcal{F}[s^2(t)](\omega)^2 = \frac{\pi}{4} \times ...\\ \bigg( \sum_{i} \sum_{j} \bigg[ e^{-2i(\phi_i + \phi_j)} \delta(\omega_i + \omega_j - \omega)^2 +& e^{2i(\phi_i - \phi_j)} \delta(\omega_i - \omega_j + \omega)^2 +\\ e^{-2i(\phi_i - \phi_j)} \delta(-\omega_i + \omega_j + \omega)^2 +& e^{2i(\phi_i + \phi_j)} \delta(\omega_i + \omega_j + \omega)^2 \bigg] \bigg) \end{align*}

So every pair $$i,j$$ contributes a unit of power at $$\omega = \omega_i - \omega_j$$ (near DC) and one $$\omega = \omega_i+\omega_j$$ (near 2x the main "carrier" frequency).

Finding out the power spectral density for $$s^2(t)$$, $$\tilde{S_2}(\omega)$$, then becomes a counting problem; how many pairs of contribute power for a given $$\omega$$? Luckily, we know the probability density function for an oscillation unit to be within $$\omega$$ and $$\omega + \delta \omega$$, given by $$\tilde{S_1}(\omega)$$.

The total probability for a pair that satisfies $$\omega = \omega_i - \omega_j$$ to be drawn from $$\tilde{S_1}(\omega)$$ is given by:

$$$$\tilde{S}_{2,DC}(\omega) = \int_{-\infty}^{\infty} \tilde{S_1}(\omega + x) \tilde{S_1}(x) dx,$$$$

The total probability for a pair that satisfies $$\omega = \omega_i + \omega_j$$ to be drawn from $$\tilde{S_1}(\omega)$$ is given by:

$$$$\tilde{S}_{2,2X/Fast}(\omega) = \int_{-\infty}^{\infty} \tilde{S_1}(\omega - x) \tilde{S_1}(x) dx,$$$$

Repeating the question: Is there a relationship between the PSD for $$s(t)$$ and the PSD for $$s(t)^2$$?

Is what I have presented correct? Does any one know about this? Is there a known theorem or something?

Remember that the Fourier transform and the inverse Fourier transform are the conjugate of each other (possibly scaled by a constant), and that the convolution theorem holds for the inverse transform as well as for the direct transform.

If the spectrum of $$s(t)$$ is $$S(\omega)$$, the PSD of $$s(t)$$ is $$|S(\omega)|^2$$. By the convolution theorem you know that the spectrum of $$s(t)^2$$ is the convolution of $$S(\omega)$$ with itself (possibly scaled). Then the PSD of $$s(t)^2$$ must be proportional to

$$\left|\int S(\omega) S(\omega-x) dx\right|^2$$

Try this with $$s(t)$$ as a sine or a cosine function, you must get impulses at $$\pm 2\omega$$ and $$0$$ for the PSD of $$s(t)^2$$ with this formula, consistently with what you expect from squaring $$s(t)^2 = (1\pm cos^2(2\omega))$$, then taking the PSD.

EDIT: Proof of the proposition of TimWescott

Proposition: For a stochastic signal, Convolving $$\sqrt{\mathbb{E}\left\{|S(\omega)|^2\right\}}$$ with itself gives an upper bound to (in my notation, square root of) the expected PSD of $$s(t)^2$$.

Results I will use without prooving

1. Jensen's inequality $$\mathbb{E}[|X|^2] \ge \mathbb{E}[|X|]^2$$, for any complex variable $$X$$.
2. From the triangle inequality $$\int |f(x)|dx \ge |\int f(x)dx|$$
3. $$\mathbb{E}\left\{X+Y\right\} = \mathbb{E}\left\{X \right\} + \mathbb{E}\left\{X \right\}$$
4. $$\mathbb{E}\left\{ A B \right\} = \mathbb{E}\left\{A\right\} \mathbb{E}\left\{ B \right\}$$

From the Jensen's inequality $$\sqrt{\mathbb{E}\left\{|S(\omega)|^2\right\}} \ge \sqrt{\mathbb{E}\left\{|S(\omega)|\right\}^2} = \mathbb{E}\left\{|S(\omega)|\right\}$$ So, the convolution of $$\sqrt{\mathbb{E}\left\{|S(\omega)|^2\right\}}$$ with itself is greater then $$\mathbb{E}\left\{|S(\omega)|\right\}$$ with itself. And I will prove that the second quantity is a (tighter) upperbound for the PSD of $$s(t)^2$$

Teorem: For a stochastic signal, convolving $$|S(\omega)|$$ with itself gives an upper bound to the expected PSD of $$s(t)^2$$.

Proof: Applying Jensen's inequaltiy, then triangle inequality

$$\begin{eqnarray}\mathbb{E} \left\{ \textrm{PSD}(s(t)^2) \right\} &=& \mathbb{E}\left\{\left|\int S(\omega) S(\omega-x) dx\right|^2\right\} \\ &\le& \mathbb{E}\left\{\int \left|S(\omega)\right| \left|S(\omega-x)\right| dx\right\}^2 \\ &=& \left(\int \mathbb{E}\left\{\left|S(\omega)\right|\left|S(\omega-x)\right|\right\} dx\right)^2 \\ &\le&\left(\int \mathbb{E}\left\{\left|S(\omega)\right|\right\} \mathbb{E}\left\{\left|S(\omega-x)\right|\right\} dx\right)^2 \end{eqnarray}$$

QED

• This assumes that you know $S(\omega)$. If $s(t)$ is a stochastic signal, then you can know the PSD (properly, you can know the expected PSD: $E \left \lbrace \left | S(\omega) \right | ^2 \right \rbrace$ But that's not the same as knowing $S(\omega)$. I'm pretty sure that if you just convolved $\sqrt{E \left \lbrace \left | S(\omega) \right | ^2 \right \rbrace}$ with itself you'd end up with an upper bound on the actual expected PSD of $s^2(t)$, not the PSD itself. Sep 24, 2021 at 20:30
• @TimWescott, I added a proof to your bound and a tighter one. I typed in a rush tough, if you find any mistakes feel free to edit. Thanks
– Bob
Sep 26, 2021 at 8:58
• For a (wide-sense-stationary) stochastic signal (which the sinusoids in the OP's question very definitely are not), the PSD is usually defined as the Fourier transform of the autocorrelation function of the signal, and not as the expected value of the $|S(\omega)|^2$. Note that $S(\omega)$ is itself a (complex-valued) stochastic process with the usual time variable $t$ being replaced by $\omega$. Sep 26, 2021 at 18:46
• @DilipSarwate: I was under the impression that $E \left \lbrace \left | S(\omega) \right |^2 \right \rbrace$ and the Fourier transform of the autocorrelation function of the signal were one and the same. If the difference can't be stated in less than one line, I may need to ask a question on here! Sep 27, 2021 at 4:03
• @TimWescott Sample paths of WSS processes are power signals, and we don't have a simple means of defining the Fourier transform of a power signal except in special cases such as sinusoidal signals, and even there, it was necessary to extend classical Fourier theory to include Dirac deltas so as to be able to include sinusoidal power signals in Fourier theory. For an arbitrary power signal $s(t)$, the Fourier transform $S(\omega)$ is might not even be definable in the Fourier theory thus extended. Sorry for extending past the one-liner requested, but at least I stayed within 500 characters! Sep 27, 2021 at 19:20