1
$\begingroup$

Assuming the signal ($x$) is in volts.

The output units of pwelch are stated to be $V^2/\text{Hz}$. I wanted to confirm this so I ran:

fs = 48000;
t = [0:1/fs:5-1/fs];
x = 0.5*cos((2*pi*1000)*t);
[pxx f] = pwelch(x, 1024, 512, 1024, fs, 'onesided');
[m,i] = max(pxx);
sqrt(m*f(i))

Which results in $1.2697$ not $0.5$. There must be something I am missing! Anyone know what it is?

$\endgroup$
  • $\begingroup$ The normalization of a spectrum estimator is kind of subject to what the designer of that algorithm wants to achieve, so your test isn't appropriate. Instead, you should take the formula (not the code), and if you really want to, track the units. But yeah, a PSD estimate should be power per bandwidth, and in the digitized sense, that would be amplitude square per frequency. But, whether that is "V² per Hz" or much more "1² per 1/f_sample" is really just an interpretation. $\endgroup$ – Marcus Müller Jul 27 '17 at 21:39
  • $\begingroup$ Pxx = pwelch(X,WINDOW,...,SPECTRUMTYPE) uses the window scaling algorithm specified by SPECTRUMTYPE when computing the power spectrum. SPECTRUMTYPE can be set to 'psd' or 'power': 'psd' - returns the power spectral density 'power' - scales each estimate of the PSD by the equivalent noise bandwidth of the window (in hertz). Use this option to obtain an estimate of the power at each frequency. The default value for SPECTRUMTYPE is 'psd'. $\endgroup$ – user28715 Jul 27 '17 at 22:03
1
$\begingroup$

@Yatekii The linked you shared helped solve the problem. I needed to take into account (1) the gain of the window used and (2) summate over a short range of frequencies. Below is the updated code:

% Create signal
    fs = 48000;t = 0:1/fs:5-1/fs;
    A = 0.5; 
    x = A*cos((2*pi*1000)*t);
    % Run FFT
    [pxx, f] = pwelch(x, 1024, 512, 1024, fs, 'one-sided');

    % Spectral parameters
    df = f(2) - f(1);
    CG = 0.54;
    NG = 0.3974;

    % Extract peak value
    isig = round(1000/df)+1;
    Aout = sqrt(sum(pxx(isig-5:isig+5)*df)) * sqrt(2);
$\endgroup$
0
$\begingroup$

Sadly I cannot comment yet but I also don't have a real answer. I can suggest you this very nice paper about matlabs pwelch: http://www.schmid-werren.ch/hanspeter/publications/2012fftnoise.pdf . Have alook at section 2.2 ;)

Good Luck!

$\endgroup$
  • $\begingroup$ The paper you suggested helped solve my riddle, more or less I missed two steps: (1) taking the windowing into effect and (2) summing over a range. Here is the updated code: % Create signal fs = 48000;t = 0:1/fs:5-1/fs; A = 0.5; x = A*cos((2*pi*1000)*t); % Run FFT [pxx, f] = pwelch(x, 1024, 512, 1024, fs, 'one-sided'); % Spectral parameters df = f(2) - f(1); CG = 0.54; NG = 0.3974; % Extract peak value isig = round(1000/df)+1; Aout = sqrt(sum(pxx(isig-5:isig+5)*df)) * sqrt(2);} $\endgroup$ – Jamie Sep 15 '17 at 20:49
  • $\begingroup$ awesome! also please be aware that the gain factors for the windows are not only dependant on window type but also window size as far as I know. I believe Hanspeter Schmid just calculated it for a single window size. $\endgroup$ – Yatekii Sep 15 '17 at 22:32

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.