67

I have no familiarity with the Multitaper method. That said, you've asked quite a question. In pursuit of my MSEE degree, I took an entire course that covered PSD estimation. The course covered all of what you listed (with exception to the Multitaper method), and also subspace methods. Even this only covers some of the main ideas, and there are many methods ...


26

The power spectral density describes the density of power in a stationary random process $X(t)$ per unit of frequency. By the Wiener-Khinchin theorem, it can be calculated as follows for a wide-sense stationary random process: $$ S_{xx}(f) = \int_{-\infty}^{\infty} r_{xx}(\tau) e^{-j2\pi f \tau} d\tau $$ where $r_{xx}(\tau)$ is the autocorrelation function ...


22

White Gaussian noise in the continuous-time case is not what is called a second-order process (meaning $E[X^2(t)]$ is finite) and so, yes, the variance is infinite. Fortunately, we can never observe a white noise process (whether Gaussian or not) in nature; it is only observable through some kind of device, e.g. a (BIBO-stable) linear filter with transfer ...


17

You are right, PSD has to do with calculating the Fourier Transform of the power of the signal and guess what.....it does. But first let's look at the mathematical relationship between the PSD and the autocorrelation function. Notations: Fourier Transform: $$ \mathcal{F}[ x(t)] = X(\omega) = \int_{-\infty}^{\infty} x(t)e^{-j\omega t}dt $$ (Time) Auto-...


13

I have no idea what your calculation of power spectral density gives since I cannot understand it. If a signal $x(t)$ has Fourier transform $X(f)$, its power spectral density is $|X(f)|^2 = S_X(f)$. The absolute spectral power in the band of frequencies from $f_0$ Hz to $f_1$ Hz is the total power in that band of frequencies, that is, the total power ...


12

I wanted to add in to the only category that the first post didn't cover. The multitaper method is a non-parametric method for computing a power spectrum similar to the periodogram approach. In this method a power spectrum is computed by windowing the data and computing a Fourier transform, taking the magnitude of the result and squaring it. The ...


9

Technically this is not a MATLAB-esque forum, but I can explain the steps in more detail for you: Suppose your input signal is $x[n]$, and its DFT is $X(f)$. For real signals you may use the one-sided DFT, since the other half would be redundant when you look at its Power Spectral Density. (PSD). Once you compute the DFT of your signal, the PSD is simply $|...


8

There is only one correct way of scaling DFT when calculating PSD with RMS values. Given input signal $x$ and its DFT $X$, the exact formula is: $$\mathrm{PSD}=\frac{2\cdot \hat{X}}{f_s\cdot S} $$ where: $\hat{X}=|X|^2=X\cdot X^*$ - squared spectrum magnitude $S=\sum_{i=1}^{N} w_{i}^{2}$ - scaling factor defined as sum of squared samples of window ...


7

If the definition of the flatness dictates that you use a power spectrum, then yes, you should square the magnitudes as the reference from the SciPy documentation indicates. In the equation that you referenced where you didn't see a squaring, I don't think you can read much into it; it says that $$ S_{flatness} = \frac{\exp\left(\frac{1}{N} \sum_k \log (...


7

As the documentation states, periodogram provides a power spectral density estimate pxx: [pxx, w] = periodogram(x); meaning that it shows how the total variance of the signal, var(x), is distributed over the frequency w. This implies that the integral of the estimated spectral density over the frequency axis trapz(w, pxx) should correspond to that total ...


7

The problem is not the spectrogram parameters, these are correct since they only depend on what resolution you want in time and frequency domain. Also, the spectrogram interpretation is correct, there are multiple frequency peaks. The problem may be: I expected to see one high power frequency after pressure rise, instead of multiple frequencies Why? If ...


6

Nice derivation but I think you can do this even easier Auto correlation $r(t) = x(t)*x(-t)$, it's the convolution of the signal with it's time flipped self. Convolution in the time domain is multiplication in the frequency domain. Time flip in the time domain is "complex conjugate" in the frequency domain. Hence we get $$ R(\omega) = \mathcal{F}\{r(t)\} ...


6

Interesting question. I juts happen to live right outside a waterfall so I can relate some hands on experience. The sound depends A LOT and water flow and weather conditions. Water flow is the main engine but wind can generate significant modulations Depending on conditions you can get strong harmonics in the signal. Sometimes you can actually see standing ...


6

$$\begin{align} R_x(\tau) &=\int_{-\infty}^\infty x(t)x^*(t-\tau)\,\mathrm dt\\ &= \int_{-\infty}^\infty \left[\int_{-\infty}^\infty X(f)e^{j2\pi ft}\,\mathrm df\right]x^*(t-\tau)\,\mathrm dt\\ &= \int_{-\infty}^\infty X(f) \left[\int_{-\infty}^\infty x^*(t-\tau) e^{j2\pi ft}\,\mathrm dt\right]\,\mathrm df\\ &= \int_{-\infty}^\infty X(f) \...


6

I don't really understand what do you mean by multiply them in the time domain and multiply them with window function. I think that you are trying to implement the Welch's PSD calculation. If so, steps should be: Split your data into possibly overlapping segments of length $N$ Preferably remove the mean of each segment. Window each of the segments by ...


6

A simple, "non-technical" way of thinking of it is the fact that the Doppler frequency is proportional to $\cos\theta$. The amplitudes of cosine, however, are not uniformly distributed, but are heavily weighted towards $\pm 1$. Example plot to demonstrate, using Python/Pylab code: theta = linspace(0, 2*pi, 1001) x = cos(theta) hist(x) More rigor can be ...


5

As Jason was saying in the comments, the power spectral density of white noise is flat. This is equivalent to saying that the autocorrelation of white noise is a delta dirac function (i.e. that there is no correlation, positive or negative, between one noise sample and another), not that the noise itself is a delta dirac function.


5

I Just do here My source code: [x, Fs, nbits] = wavread('ederwander.wav'); winSize = 2048; n_samples = length(x); %50% overlap or 0 to not use overlap OverlapStep = 50; if OverlapStep > 0 Overlap = floor((OverlapStep*winSize) / 100); nFrames=floor(n_samples/Overlap)-1; else Overlap= winSize; ...


5

The issue (beside the scaling) here is that your selected frequencies do not align with your FFT grid and that the energy of each complex exponential is smeared out over multiple FFT bins. Your frequency resolution is given by sample rate divide by FFT length. In your case that's 3.906 Hz. Only for frequencies that are an integer multiple of this will all ...


5

There are several issues with your code: The "Gibbs" issue is a non-issue, so you shouldn't need to window. You are better off generating the complete $-\pi$ to $+\pi$ spectrum that you want, rather than relying on getting the symmetry correct. The code FFT = magnitude-np.mean(magnitude)+phase is changing the phase of the PSD magnitude by addition. It ...


5

In principle you are doing the right thing. Step 4 should produce a real result unless there is a coding error. Sometimes you add up some residual imaginary part due to numerical noise but, if any, that should be very small. Here is an example %% random signal of length n n = 128; x = rand(n,1); % zero pad x = [x; zeros(n,1)]; % fft fx = fft(x); % mag ...


5

Well, because in a lot of real world problems this If I have a deterministic signal with a fixed number of samples, shouldn't I be able to directly determine its spectral information? is just not the case. Very often, measured signals are more of a random process. A simple and common case would be to have the desired signal and some additive noise, very ...


5

You can compute the power of the process from its power spectrum as well as from its PDF. Equating the two gives you a relation between the constants $A$ and $B$. More specifically you get $$\int_{-B}^BG_x(f)df=\frac{1}{2A}\int_{-A}^Ax^2dx$$ If I'm not mistaken this should give $A=\sqrt{3B}$.


5

As mentioned by @MBaz in comments, you could use the fact that filtering white noise with unitary power spectral density (for example from independent Gaussian samples zero mean and unitary standard deviation) with a filter having a frequency response $H(f)$ gives you colored noise with power spectral density $S_{xx}(f) = |H(f)|^2$. Then comes the question ...


5

In addition to Carlos answer, I want to correct your general understanding: What I understand of Doppler spread is that the relative motion between Transmitter (TX) and Receiver (RX) change the exposing time of signal. In rapport to a constant-distance TX-RX, a moving toward each other TX-RX "compresses" signal in time (signal takes less time to ...


5

That equation is indeed simple, but it's also wrong. So your doubts are completely justified. What is probably meant is something like $$S_X(f)=\lim_{T\rightarrow\infty}\frac{1}{T}\left|X_T(f)\right|^2\tag{1}$$ with $$X_T(f)=\int_{-T/2}^{T/2}x(t)e^{-j2\pi f t}dt\tag{2}$$ This definition of the power spectrum is used for functions $x(t)$ with finite power ...


4

Suppose we have a discrete-time sequence $x[t]$ which is stationary, zero mean, white noise with variance $\sigma^2$. Then the autocorrelation of $x$ is: $$ \begin{array} RR_{xx}[\tau] &=& E\left[ x[t] x[t+\tau] \right]\\ &=& \left \{ \begin{array} EE \left[ x[t]^2 \right], {\rm if\ }\tau=0 \\ 0, {\rm otherwise} \end{array} \right. \\ &=&...


4

The most authoritative reference I can come up with is from Jayant & Noll, Digital Coding Of Waveforms, (c) Bell Telephone Laboratories, Incorporated 1984, published by Prentice-Hall, Inc. On page 57, they define the spectral flatness: and, previously, on page 55 they define $S_{xx}$: So the FFT-squared version is the one you want. It looks like ...


4

Definitions vary, don't they? The first thing that has to be settled is whether we agree that the power spectral density is equivalent to the power spectrum, or else define what we mean by both. Proakis and Salehi use them synonymously. Moving on, I think the discrepancies are due to differing definitions, for signals that have one, of the power spectrum. ...


4

I am not sure I understand your question, nevertheless, I would like to point out that matlab computes fft considering this expression: $$ D_{n}=\sum_{k=0}^{K-1}x(k)\cdot e^{-j\cdot n\cdot2\cdot\pi\cdot\frac{k}{K}} $$ Which as you know it misses a factor of $$\frac{1}{K} $$ That is: $$ D_{n}=\frac{1}{K}\sum_{k=0}^{K-1}x(k)\cdot e^{-j\cdot n\cdot2\cdot\...


Only top voted, non community-wiki answers of a minimum length are eligible