Mercurial > hg > ltpda
comparison m-toolbox/classes/+utils/@math/computepsd.m @ 0:f0afece42f48
Import.
| author | Daniele Nicolodi <nicolodi@science.unitn.it> |
|---|---|
| date | Wed, 23 Nov 2011 19:22:13 +0100 |
| parents | |
| children |
comparison
equal
deleted
inserted
replaced
| -1:000000000000 | 0:f0afece42f48 |
|---|---|
| 1 function varargout = computepsd(Sxx,Svxx,w,range,nfft,Fs,esttype) | |
| 2 % | |
| 3 % Slight modification of original MATLAB's computepsd to include correct | |
| 4 % scaling for the variance, i.e var(a*x) = a^2*var(x) | |
| 5 % | |
| 6 % VERSION: $Id: computepsd.m,v 1.1 2009/07/22 15:29:59 miquel Exp $ | |
| 7 % | |
| 8 % HISTORY: 21-07-09 M Nofrarias | |
| 9 % Creation | |
| 10 % | |
| 11 %%%%%%%%%%%%%%%%% | |
| 12 % | |
| 13 % COMPUTEPSD Compute the one-sided or two-sided PSD or Mean-Square. | |
| 14 % [Pxx,Pvxx,W,UNITS] = COMPUTEPSD(Sxx,Svxx,W,RANGE,NFFT,Fs,ESTTYPE) where the | |
| 15 % inputs and outputs are: | |
| 16 % | |
| 17 % Inputs: | |
| 18 % Sxx - Whole power spectrum [Power]; it can be a vector or a matrix. | |
| 19 % For matrices the operation is applied to each column. | |
| 20 % W - Frequency vector in rad/sample or in Hz. | |
| 21 % RANGE - Determines if a 'onesided' or a 'twosided' Pxx and Sxx are | |
| 22 % returned. | |
| 23 % NFFT - Number of frequency points. | |
| 24 % Fs - Sampling Frequency. | |
| 25 % ESTTYPE - A string indicating the estimate type: 'psd', or 'ms' value. | |
| 26 % | |
| 27 % Outputs: | |
| 28 % Pxx - One-sided or two-sided PSD or MEAN-SQUARE (not scaled by Fs) | |
| 29 % depending on the input arguments RANGE and TYPE. | |
| 30 % W - Frequency vector 0 to 2*Nyquist or 0 to Nyquist depending on | |
| 31 % range, units will be either rad/sample (if Fs is empty) or Hz | |
| 32 % (otherwise). | |
| 33 % UNITS - Either 'rad/sample' or 'Hz' | |
| 34 % Author(s): R. Losada | |
| 35 % Copyright 1988-2005 The MathWorks, Inc. | |
| 36 % $Revision: 1.1 $ $Date: 2009/07/22 15:29:59 $ | |
| 37 | |
| 38 if nargin < 7 | |
| 39 iswinfs = false; | |
| 40 if nargin < 6, | |
| 41 esttype = 'psd'; | |
| 42 end | |
| 43 end | |
| 44 | |
| 45 w = w(:); % Make sure we always returns a column vector for frequency | |
| 46 | |
| 47 % Generate the one-sided spectrum [Power] if so wanted | |
| 48 if strcmp(range,'onesided'), | |
| 49 if rem(nfft,2), | |
| 50 select = 1:(nfft+1)/2; % ODD | |
| 51 Sxx_unscaled = Sxx(select,:); % Take only [0,pi] or [0,pi) | |
| 52 Svxx_unscaled = Svxx(select,:); % Take only [0,pi] or [0,pi) | |
| 53 Sxx = [Sxx_unscaled(1,:); 2*Sxx_unscaled(2:end,:)]; % Only DC is a unique point and doesn't get doubled | |
| 54 Svxx = [Svxx_unscaled(1,:); 4*Svxx_unscaled(2:end,:)]; | |
| 55 else | |
| 56 select = 1:nfft/2+1; % EVEN | |
| 57 Sxx_unscaled = Sxx(select,:); % Take only [0,pi] or [0,pi) | |
| 58 Svxx_unscaled = Svxx(select,:); % Take only [0,pi] or [0,pi) | |
| 59 Sxx = [Sxx_unscaled(1,:); 2*Sxx_unscaled(2:end-1,:); Sxx_unscaled(end,:)]; % Don't double unique Nyquist point | |
| 60 Svxx = [Svxx_unscaled(1,:); 4*Svxx_unscaled(2:end-1,:); Svxx_unscaled(end,:)]; % Don't double unique Nyquist point | |
| 61 end | |
| 62 w = w(select); | |
| 63 end | |
| 64 | |
| 65 % Compute the PSD [Power/freq] | |
| 66 if ~isempty(Fs), | |
| 67 Pxx = Sxx./Fs; % Scale by the sampling frequency to obtain the psd | |
| 68 Pvxx = Svxx./Fs^2; | |
| 69 units = 'Hz'; | |
| 70 else | |
| 71 Pxx = Sxx./(2.*pi); % Scale the power spectrum by 2*pi to obtain the psd | |
| 72 Pvxx = Svxx./(2.*pi)^2; | |
| 73 units = 'rad/sample'; | |
| 74 end | |
| 75 | |
| 76 if strcmpi(esttype,'ms'), | |
| 77 varargout = {Sxx,Svxx,w,units}; % Mean-square | |
| 78 else | |
| 79 varargout = {Pxx,Pvxx,w,units}; % PSD | |
| 80 end | |
| 81 | |
| 82 % [EOF] computepsd.m |