Determining power spectral density from a sample

Question

To start, I have absolutely no background in electrical engineering, I am a computer scientist by trade.

I'm currently working on a project where I am attempting to code a function that uses the power spectrum density of a particular signal at a given frequency.

Specifically, I want to calculate the following function

$$F(x,f) = \int [P(f) * X(x,f)] df$$

where $F$ is the function I want to calculate, $f$ is the frequency, $X$ is an function that doesn't pertain to my question, $x$ is a variable that also doesn't pertain to the question, and $P$ is the power spectral density of a signal.

I don't have a function describing the signal, instead I have a vector of samples of the signal. What I'm interested in calculating is $P$ using my sample vector where $P$ has been normalized such that

$$ \int [P(f)] df = 1$$

Questions:

1. How might I go about this in MATLAB? I have looked at the pwelch function (Welch's power spectrum density estimate) and I think it does what I want. However, the values seem to be scaled to $[0,\pi]$ which makes me a bit wary. In addition, just calling the function will not scale the values so that $ \int [P(f)] df = 1$. Is using pwelch an acceptable approach for what I'm trying to do and how can I get things scaled correctly?

2. In addition to Computer Science, I also have a background in math. However, I'm a bit thrown off by the integrals that I have been provided since the result has been "evaluated" despite no upper and lower limits ($F$ should evaluate to an actual value). Is it common notation to integrate from negative to positive infinity in signal processing using this notation? I understand that this question is a bit vague since I haven't really explained much of what I am trying to do. This question is meant more as a general notation question regarding convention in the signal processing / electrical engineering community. Specifically, is there a common convention as to what limits to integrate over when the limits have been omitted?

Answer 1

In Matlab, you are dealing with sampled signals, so your frequency variable goes from $0$ to Nyquist. This is why the x-axis is labeled in normalized frequency that goes from $0$ to $\pi$. If you like you can use [P, f] = pwelch(___, fs) where fs is your sampling rate, and then the function will return the power spectral density (PSD) estimate P as a function of analog frequencies in the vector f. See example here.

pwelch is not the only method for estimating PSD. The algorithm choice will depend on your application. A nice summary of other competing algorithms can be found in the answer to this question.

As for these integrals, the limits should go from $-\infty$ to $\infty$ for continuous signals. But if you know that the signal is bandlimited to some range of frequencies $(f_1, f_2)$ then you can get away with integrating from $f_1$ to $f_2$. Recall the definition of PSD is inspired by continuous time signals $x(t)$ with finite energy: $$ E = \int_{-\infty}^{\infty} |x(t)|^2 dt < \infty $$ and denoting $X(f)$ to be the Fourier transform of $x(t)$, by Parseval's theorem we have: $$ \int_{-\infty}^{\infty} |x(t)|^2 dt = \int_{-\infty}^{\infty} |X(f)|^2 df $$ This leads to the intuition that $|X(f)|^2$ is the "density" of "spectral energy" around the frequency $f.$

Requiring that the PSD integrate to 1 implies that your signals of interest are unit energy. For a non-unit-energy finite energy signal, you can simply divide the signal $x(t)$ by $\sqrt{E}$ to transform it into a unit energy signal.

For implementation purposes, all your signals and transforms are discrete. So you will have to approximate these integrals with sums.