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?

| improve this question | |
$\endgroup$
  • $\begingroup$ What is your sample rate? $\endgroup$ – Robert L. Nov 18 '16 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 '16 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 '16 at 14:47
  • $\begingroup$ I edited the question to add the sample rate, 20000. $\endgroup$ – pscheidler Nov 18 '16 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 '16 at 23:30
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).

| improve this answer | |
$\endgroup$
  • $\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 '16 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$ – Maximilian Matthé Nov 18 '16 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 '16 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$ – Maximilian Matthé Nov 18 '16 at 17:34

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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