2
$\begingroup$

I'm having some trouble getting consistent results with PSD calculations. If I just take a basic sine wave

fs = 20000;
len = 16384;
fsig = 1000;
Sig = 1.5 * sin( (1:len) * 2 * pi * fsig / fs);
fsig = 700;
Sig = Sig + 0.9 * sin( (1:len) * 2 * pi * fsig / fs);

and I run it through the follow two Welch's functions with MATLAB's pwelch:

[matWelchD1, matWelchF1] = pwelch(Sig, length(Sig), 0, [], fs);
[matWelchD2, matWelchF2] = pwelch(Sig, ceil(length(Sig/2)), 0, [], fs);

I get the $1\textrm{ kHz}$ peak at approximately $0.65$ and $0.26$. I would love to blame frequency leakage due to how the frequencies line up, and there is some support for this because the $0.65$ peak has a much narrower base, but I see a similar result on real world data with FFTs that have broadband frequencies.

Anyone know why the two implementations would be different?

$\endgroup$
5
  • $\begingroup$ What is your sample rate? $\endgroup$
    – Robert L.
    Nov 18, 2016 at 1:15
  • 1
    $\begingroup$ In the past, i performed the SAME study you already did. I had a lot of problems by using welch and those familes. I would recommend you to drop it, and to use the standard FFT. The welch methods are not properly configured, and it is equivalent to use a frequency filtered FFT instead the welch inconsistencies...... $\endgroup$
    – Brethlosze
    Nov 18, 2016 at 2:01
  • $\begingroup$ @hypfco Thanks, that is what I was afraid of. Can you expand on "not properly configured" and what kind of filtering it applies to the FFT? $\endgroup$
    – pscheidler
    Nov 18, 2016 at 14:47
  • $\begingroup$ I edited the question to add the sample rate, 20000. $\endgroup$
    – pscheidler
    Nov 18, 2016 at 14:49
  • $\begingroup$ I will check the accepted answer and my study. When i didn that, i included noise signals, and i didnt get the pretty results of the answer. $\endgroup$
    – Brethlosze
    Nov 18, 2016 at 23:30

1 Answer 1

1
$\begingroup$

Your value for $f_s$ is missing. Assuming you work above Nyquist sampling rate (i.e. $f_s > 2*f_{sig}$), I obtain the following results:

fs = 20000;
len = 16384;
fsig = 1000;
Sig = 1.5 * sin( (1:len) * 2 * pi * fsig / fs);
fsig = 700;
Sig = Sig + 0.9 * sin( (1:len) * 2 * pi * fsig / fs);

[matWelchD1, matWelchF1] = pwelch(Sig, length(Sig), 0, [], fs);
[matWelchD2, matWelchF2] = pwelch(Sig, ceil(length(Sig)/2), 0, [], fs);

subplot(2,1,1);
plot(matWelchF1, (matWelchD1), '-o'); grid;

subplot(2,1,2);
plot(matWelchF2, (matWelchD2), '-o'); grid;

df1 = matWelchF1(2) - matWelchF1(1);
df2 = matWelchF2(2) - matWelchF2(1);

df1
df2

Output:

df1 =

    1.2207


df2 =

    2.4414

Program output

As you can see, both frequencies are nicely represented in the figures.

One more comments on the window length (i.e. second parameter to pwelch): It determines how big the sections are, into which the signal is divided. Each section is multiplied with a Hamming window and then the FFT is taken. Afterwards, the values of all FFTs are summed together, yielding the PSD estimate.

I.e. putting the signal length as the window length, would result in a single section. For your stationary signal this is fine, but in reality you would want to have a smaller value, or leave the parameter out (or use the empty vector []), such that the PSD can be more accurately estimated.

As you can see, different window length create a different amount of frequency samples: The longer the window, the more frequency samples. NOte that pwelch calculates the Power Spectral Density, i.e. the power per Hertz. In order to get the energy of a frequency you need to multiply with the bandwidth of each bin:

>> sum(matWelchD1) * (matWelchF1(2)-matWelchF1(1))

ans =

    1.5300

>> sum(matWelchD2) * (matWelchF2(2)-matWelchF2(1))

ans =

    1.5300

Both PSDs contain the same overall energy. Since the frequency samples in the first PSD are closer to each other, the PSD has higher peaks for each frequency (to deliver the same overall power).

$\endgroup$
4
  • $\begingroup$ I edited my code to show that Fs was 20 kHz. I understand how the windowing works, but I don't have a good intuition for why one result is over twice the magnitude of the other result. The result is the same in your graphs, you just hide the difference with log plots. I do understand that log plots are the normal way to display this kind of data. $\endgroup$
    – pscheidler
    Nov 18, 2016 at 14:45
  • $\begingroup$ I have edited my answer to explain the difference in the amplitudes of the PSD. It's related to the D in PSD: "Density", i.e. energy/Hertz $\endgroup$ Nov 18, 2016 at 15:28
  • $\begingroup$ That makes sense! Thanks! And I assume that the 1.53 value comes from the square of the RMS energy of the 0.9 and 1.5 amplitude sine waves, which would be 0.5 * (1.5^2 + 0.9^2), since the RMS = Peak/sqrt(2). Is it common to just multiply each point by the bin size? It seems like that would give a better representation. $\endgroup$
    – pscheidler
    Nov 18, 2016 at 16:02
  • $\begingroup$ You can multiply the points with the bin sizes. This would then correspond more to a signal that has a discrete spectrum at the given frequency bins (and zero between them) instead of a continuous spectrum (which has the energy distributed among all real-valued frequencies). $\endgroup$ Nov 18, 2016 at 17:34

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.