[Octave-bug-tracker] [bug #62095] /usr/local/share/octave/4.2.1/m/signal/periodogram.m (version 1.4.1) gives wrong result (namely 0)

[anonymous]

URL:
  <https://savannah.gnu.org/bugs/?62095>

                 Summary: /usr/local/share/octave/4.2.1/m/signal/periodogram.m
(version 1.4.1) gives wrong result (namely 0)
                 Project: GNU Octave
            Submitted by: None
            Submitted on: Tue 22 Feb 2022 12:31:18 PM UTC
                Category: Octave Forge Package
                Severity: 3 - Normal
                Priority: 5 - Normal
              Item Group: Incorrect Result
                  Status: None
             Assigned to: None
         Originator Name: Tim van Zon
        Originator Email: t.vanzon@prolira.com
             Open/Closed: Open
                 Release: 4.2.1
         Discussion Lock: Any
        Operating System: GNU/Linux

    _______________________________________________________

Details:

I constructed an example where periodogram gives the wrong answer. The bug is
that periodogram averages the segments before taking the fft, instead of
averaging the fft output.

++ verbatim ++


N = 8;  % windowsize
Fs_Hz = 10; % sample rate in Hz
dt_sec = 1./Fs_Hz; 
a = randn(N,1);

% make a signal by replicating 'a' with a minus sign
signal = [a;-1.*a];

pxx = periodogram(signal,[],N,Fs_Hz)
% this will result in sum(pxx)==0, while it should not be zero, as there is
amplitude there.
% this is because periodogram averages the segments BEFORE taking the fft.
% this example is constructed to show this, as the average of the 2 segments
is now zero.
 
% And there is a second bug, namely that the documentation says it gives the 
% POWER spectral density (PSD), while it actually gives the ENERGY spectral
density.
% (please recall that  power * time = energy) 
energy_a = sum(a.*conj(a) * dt_sec);  % Not necessary to use conj here. so
this is the Es from https://en.wikipedia.org/wiki/Energy_(signal_processing)
pA = periodogram(a,[],N,Fs_Hz);
% summing over all frequencies in pxx should give the same number as
energy_a,
assert( sum(pA), energy_a, 1e-9) % so this is as expected

% Note: do not compare against pwelch, as that will give the energy per
windowlength, which can 
% be different from energy and from power! 



-- verbatim -- 




    _______________________________________________________

Reply to this item at:

  <https://savannah.gnu.org/bugs/?62095>

_______________________________________________
  Message sent via Savannah
  https://savannah.gnu.org/