diff options
Diffstat (limited to 'macros/periodogram.sci')
| -rwxr-xr-x | macros/periodogram.sci | 225 |
1 files changed, 183 insertions, 42 deletions
diff --git a/macros/periodogram.sci b/macros/periodogram.sci index ca8995e..11971d6 100755 --- a/macros/periodogram.sci +++ b/macros/periodogram.sci @@ -1,46 +1,187 @@ -function [d,n]=periodogram(a,b,c,d,e) -//Return the periodogram (Power Spectral Density) of X -//Calling Sequence -// [PXX, W] = periodogram (X) -// [PXX, W] = periodogram (X, WIN) -// [PXX, W] = periodogram (X, WIN, NFFT) -// [PXX, W] = periodogram (X, WIN, NFFT, FS) -// [PXX, W] = periodogram (..., "RANGE") -//Parameters -// 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 th 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]. -//Description -//The optional second output W are the normalized angular frequencies. For a one-sided calculation W is in the range [0, pi]. If NFFT is even and [0, pi) if NFFT is odd. Similarly, for a two-sided calculation W is in the range [0, 2*pi] or [0, 2*pi)depending on NFFT. -// -//If a sampling frequency is specified, FS, then the output frequencies F will be in the range [0, FS/2] or [0, FS/2) for one-sided calculations. For two-sided calculations the range will be [0, FS). -// -//When called with no outputs the periodogram is immediately plotted in the current figure window. - funcprot(0); - lhs= argn(1); - rhs= argn(2); - if(rhs<1 | rhs>5) - error("Wrong number of input arguments"); +// 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 +// Author:[insert name] +// 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]. +// +// + + +//Test cases: +////1. +//n=0:319; +//x=cos(%pi/4*n)+rand(size(n,"r"),"normal"); +//[pxx,w]=periodogram(x,ones(1,320),256,2000,"onesided"); +//plot2d(w,10*log10(pxx)) +//xtitle('periodogram','frequency','magnitude(db)') +//xgrid() +// +// + +[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 - if(lhs>2 | lhs< 2) - error("Wrong number of output arguments"); + 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,"r"); + + if (~isempty (win)) + if (~ isvector (win) | length (win) ~= n) + error ("periodogram: WIN must be a vector of the same length as X"); end - select(rhs) - case 1 then - [d,n]= callOctave('periodogram',a); - case 2 then - [d,n]= callOctave('periodogram',a,b); - - case 3 then - [d,n]= callOctave('periodogram',a,b,c); - - case 4 then - [d,n]= callOctave('periodogram',a,b,c,d); - - case 5 then - [d,n]= callOctave('periodogram',a,b,c,d,e); + 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 (~strcmpi (ran, "onesided")) + ran = 1; + elseif (~strcmpi (ran, "twosided")) + ran = 2; + elseif (~strcmpi (ran, "centered")) + error ('periodogram: centered ran type is not implemented'); + else + ran = 2-isreal (x); + end + + // compute periodogram + + if (n > nfft) + Pxx = 0; + rr = modulo(length (x), nfft); + if (rr) + x = [x(:); zeros(nfft-rr, 1)]; end -endfunction + x = sum (matrix (x, nfft,length(x)/nfft), 2); + end + + if (~ isempty (win)) + n = sum(win.*conj(win)); + end; + Pxx = (abs (fft (x))) .^ 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) + 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 + + + + pxx = Pxx; + + + + +endfunction |