diff options
Diffstat (limited to 'macros/periodgram.sci')
| -rw-r--r-- | macros/periodgram.sci | 232 |
1 files changed, 232 insertions, 0 deletions
diff --git a/macros/periodgram.sci b/macros/periodgram.sci new file mode 100644 index 0000000..0d406ec --- /dev/null +++ b/macros/periodgram.sci @@ -0,0 +1,232 @@ +// Copyright (C) 2018 - IIT Bombay - FOSSEE +// This file must be used under the terms of the CeCILL. +// This source file is licensed as described in the file COPYING, which +// you should have received as part of this distribution. The terms +// are also available at +// http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt +// Original Source : https://octave.sourceforge.io/signal/ +// Modifieded by: Abinash Singh Under FOSSEE Internship +// Last Modified on : Feb 2024 +// Organization: FOSSEE, IIT Bombay +// Email: toolbox@scilab.in + +function [pxx, f] = periodogram (x, varargin) + //Calling Sequence: + //[PXX, F] = periodogram (X, WIN, NFFT, FS) + //[PXX, F] = periodogram (..., "RANGE") + + // The possible inputs are: + // + // X + // + // data vector. If X is real-valued a one-sided spectrum is + // estimated. If X is complex-valued, or "RANGE" specifies + // "twosided", the full spectrum is estimated. + // + // WIN + // window weight data. If window is empty or unspecified a + // default rectangular window is used. Otherwise, the window is + // applied to the signal ('X .* WIN') before computing the + // periodogram. The window data must be a vector of the same + // length as X. + // + // NFFT + // number of frequency bins. The default is 256 or the next + // higher power of 2 greater than the length of X ('max (256, + // 2.^nextpow2 (length (x)))'). If NFFT is greater than the + // length of the input then X will be zero-padded to the length + // of NFFT. + // + // FS + // sampling rate. The default is 1. + // + // RANGE + // range of spectrum. "onesided" computes spectrum from + // [0..nfft/2+1]. "twosided" computes spectrum from [0..nfft-1]. + // + // + // Dependencies + // hamming fft1 + [nargout,nargin]=argn(); + // check input arguments + if (nargin < 1 | nargin > 5) + error("wrong no. of input arguments") + end + + nfft = []; + fs = []; + ran = ""; + win = []; + j = 2; + for k = 1:length (varargin) + if (type (varargin(k))==10) + ran = varargin(k); + else + select (j) + case 2 + win = varargin(k); + case 3 + nfft = varargin(k); + case 4 + fs = varargin(k); + end + j=j+1; + end + end + + if (~ isvector (x)) + error ("periodogram: X must be a real or complex vector"); + end + x = x(:); // Use column vectors from now on + + n = size(x,1); + + if (~isempty (win)) + if (~ isvector (win) | length (win) ~= n) + error ("periodogram: WIN must be a vector of the same length as X"); + end + win = win(:); + x =x.* win; + else + win=window("re",length(x)); + win=win(:); + x=x.*win; + + end + + if (isempty (nfft)) + nfft = max (256, 2.^nextpow2 (n)); + elseif (~ isscalar (nfft)) + error ("periodogram: NFFT must be a scalar"); + end + + use_w_freq = isempty (fs); + if (~use_w_freq & ~ isscalar (fs)) + error ("periodogram: FS must be a scalar"); + end + + if (~strcmp (ran, "onesided")) + ran = 1; + elseif (~strcmp (ran, "twosided")) + ran = 2; + elseif (~strcmp (ran, "centered")) + error ('periodogram: centered ran type is not implemented'); + else + ran = 2 - double(isreal (x)); + end + + // compute periodogram + + if (n > nfft) + Pxx = 0; + rr = modulo(length(x), nfft); + if (rr) + x = [x(:); zeros(nfft-rr, 1)]; + end + x = sum (matrix (x, nfft,-1), 2); + end + + if (~ isempty (win)) + n = sum(win.*conj(win)); + end; + + Pxx = (abs (fft1 (x,nfft))) .^ 2 / n; + + if (use_w_freq) + Pxx =Pxx/(2*%pi); + else + Pxx =Pxx/fs; + end + + // generate output arguments + + if (ran == 1) // onesided + if (modulo(nfft,2)==0) // nfft is even + psd_len = (nfft/2)+1; + Pxx = Pxx(1:psd_len) + [0; Pxx(nfft:-1:psd_len+1); 0]; + else // nfft is odd + psd_len = (nfft+1)/2; + Pxx = Pxx(1:psd_len) + [0; Pxx(nfft:-1:psd_len+1)]; + end + end + + //if (nargout() ~= 1) FIXME: fix nargout + if (ran == 1) + f = (0:nfft/2)' / nfft; + elseif (ran == 2) + f = (0:nfft-1)' / nfft; + end + if (use_w_freq) + f =f* 2*%pi; // generate w=2*pi*f + else + f =f* fs; + end + //end + if (nargout() ~= 2) + if (use_w_freq) + plot (f/(2*%pi), 10*log10 (Pxx)); + xlabel ("normalized frequency [x pi rad]"); + ylabel ("Power density [dB/rad/sample]"); + else + plot (f, 10*log10 (Pxx)); + xlabel ("frequency [Hz]"); + ylabel ("Power density [dB/Hz]"); + end + title ("Periodogram Power Spectral Density Estimate"); + end + pxx = Pxx; +endfunction + +/* +pi = %pi; // ezecute on scilab only + +t = 0:0.01:1; +x = sin(2*pi*10*t); +periodogram(x); + +x = complex(0:0.01:1, 0:0.01:1); +periodogram(x); + +x = cos(0:0.01:1); +win = hamming(101); +periodogram(x, win); + + + +x = tan(0:0.01:1); +nfft = 512; +periodogram(x, [], nfft); + + +t = 0:0.01:1; +x = sin(2*pi*10*t); +Fs = 100; +periodogram(x, [], [], Fs) + + + +x = sin(0:0.01:1); +periodogram(x, [], [], [], 'onesided'); + +x = sin(0:0.01:1); +periodogram(x, [], [], [], 'twosided') + + +Fs = 1000; +t = 0:1/Fs:1-1/Fs; +f0 = 100; +x = sin(2*pi*f0*t); +[Pxx, f] = periodogram(x, [], [], Fs); +[_, idx] = max(Pxx); +detected_freq = f(idx); + + +// Test error : invalid window length +x = randn(100,1); +win = hamming(50); +periodogram(x, win) + +// Test invalid nfft (negative) +periodogram(x, [], -256) + +*/ |