# PSD subtraction

I am interested seeing the difference between two power spectral densities (PSD) as a noise reduction exercise

The blue line is the psd of my signal, and the the orange line is the psd of the signal without the part I am interested in. So to isolate the part I am interested in I would like to 'subtract' these to see the difference.

However it occurs to me that this would be dangerous since if we consider the DFT's of the 'blue and orange signals' for a particular frequency $$k$$, to be $$F_k$$ and $$G_k$$ respectively. Then the PSD of one minus the other is not the same as the PSD of one minus the PSD of the other since by taking the PSD you lose the phase information. Or less confusingly the problem boils down to this:

$$|F_k|^2 - |G_k|^2 \neq |F_k - G_k|^2$$

The obvious solution to this is to simply subtract the output of the DFT step, then PSD. However in python this will be difficult since I am using scipys version of welches method to calculate the PSDs which means I would have to rewrite a lot of my code.

So my question is this:

Is there a way to get a good approximation to the difference between the two PSDs using the PSDs rather than doing the earlier subtraction? Or again less confusingly:

Can I get an approximation to $$|F_k - G_k|^2$$ from $$|F_k|^2$$ and $$|G_k|^2$$?

• Maybe I just don't understand, but you have the original time-domain signals, right? Why not just subtract those, and then calculate the PSD (with Welch or whatever method of your choosing)? As you mention, the PSD is nonlinear, but the DFT is, so just subtracting in the time domain before should work. – mateC Sep 9 at 8:49
• Yeah, you're right I think except I can't guarantee my time signals are exactly in phase, so I would really be doing $F_k[y(t) + x(t + \tau)] \neq F_k[y(t) + x(t)]$ for unknown $\tau$. I think I will have to just change the code to allow for cross correlations – Stephen Jackson Sep 9 at 10:34

In general, if you have a signal $$x$$ corrupted by noise $$n$$ that is known (or assumed to be) independent of the signal, then it is true that the PSD $$S_{x+n}(f)$$ of $$x+n$$ is the sum of the PSDs $$S_x(f)$$ and $$S_n(f)$$ of $$x$$ and $$n$$ respectively: $$S_{x+n}(f) = S_x(f) + S_n(f).\tag{1}$$ This is because the cross-power spectral densities $$S_{x,n}(f)$$ and $$S_{n,x}(f)$$ that should appear on the right side of $$(1)$$ have value $$0$$ because $$x$$ and $$n$$ are independent. So, if you have measured $$S_{x+n}(f)$$ and have some assumptions about $$S_n(f)$$, it is perfectly legitimate to take $$S_x(f)$$ to be the difference between $$S_{x+n}(f)$$ and $$S_{n}(f)$$. Be aware, though, that this procedure can be fraught with peril when your measurements of $$S_{x+n}(f)$$ aren't quite on target and your assumptions about $$S_{n}(f)$$ are flaky and so the purported $$S_x(f) = S_{x+n}(f) - S_n(f)$$ works out to be negative for some values of $$f$$. When in doubt, blame the instruments!
• Also can I do $S_x(f) = S_{x+n}(f) - S_n(f) - S_{x,n}(f) - S_{n,x}(f)$ If the two signals are dependent? – Stephen Jackson Apr 11 at 10:24
• Power spectral densities are Fourier transforms of autocorrelation functions $R$, and so, we have that $$R_{x+n}(\tau)=E[(x(t)+n(t))(x(t+\tau)+n(t+\tau))]=R_x(\tau)+R_n(\tau)+R_{x,n}(\tau)+R_{n,x}(\tau)$$ using FOIL and so the power spectral density result follows. When $x$ and $n$ are independent (and zero-mean), their cross-correlation function is $0$ leading to $(1)$ above. – Dilip Sarwate Apr 12 at 1:35