Home > m > sigproc > frequency_domain > ltpda_psd.m

ltpda_psd

PURPOSE ^

LTPDA_PSD a psd estimation based on pwelch from MATLAB. The difference is that

SYNOPSIS ^

function varargout = ltpda_psd(x, win, nolap, nfft, fs, scale, inunits)

DESCRIPTION ^

 LTPDA_PSD a psd estimation based on pwelch from MATLAB. The difference is that
 linear detrending is done on each windowed data segment.
 
 function [f, pxx, info] = ltpda_psd(x, win, nolap, nfft, fs, scale, inunits)
 
 
 Inputs:
     x      - data vector
     win    - a window function
              Specify a vector of window samples or a length. If a length
              is specified, a Hamming window is used.
     nfft   - number of samples per fft
     nolap  - number of samples to overlap each segment by
     fs     - sample rate of input time-series
     scale  - specify one of 'ASD', 'AS', 'PSD', 'PS'
   inunits  - Amplitude units of input time-series data.
 
 Outputs:
     f      - a frequency vector 
     Sxx    - spectral estimate in units specified by 'scale'
     info   - a structure of information
 
 NOTE: segment_info() is copied directly from MATLAB's welshparse.m.
 
 M Hewitson 25-09-06
 
 $Id: ltpda_psd.html,v 1.1 2007/06/08 14:15:10 hewitson Exp $

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SUBFUNCTIONS ^

SOURCE CODE ^

0001 function varargout = ltpda_psd(x, win, nolap, nfft, fs, scale, inunits)
0002 
0003 % LTPDA_PSD a psd estimation based on pwelch from MATLAB. The difference is that
0004 % linear detrending is done on each windowed data segment.
0005 %
0006 % function [f, pxx, info] = ltpda_psd(x, win, nolap, nfft, fs, scale, inunits)
0007 %
0008 %
0009 % Inputs:
0010 %     x      - data vector
0011 %     win    - a window function
0012 %              Specify a vector of window samples or a length. If a length
0013 %              is specified, a Hamming window is used.
0014 %     nfft   - number of samples per fft
0015 %     nolap  - number of samples to overlap each segment by
0016 %     fs     - sample rate of input time-series
0017 %     scale  - specify one of 'ASD', 'AS', 'PSD', 'PS'
0018 %   inunits  - Amplitude units of input time-series data.
0019 %
0020 % Outputs:
0021 %     f      - a frequency vector
0022 %     Sxx    - spectral estimate in units specified by 'scale'
0023 %     info   - a structure of information
0024 %
0025 % NOTE: segment_info() is copied directly from MATLAB's welshparse.m.
0026 %
0027 % M Hewitson 25-09-06
0028 %
0029 % $Id: ltpda_psd.html,v 1.1 2007/06/08 14:15:10 hewitson Exp $
0030 %
0031 
0032 % get segment information
0033 Nx = length(x);
0034 [L,noverlap,win,msg] = segment_info(Nx,win,nolap);
0035 % Compute the number of segments
0036 k = (Nx-noverlap)./(L-noverlap);
0037 
0038 % Uncomment the following line to produce a warning each time the data
0039 % segmentation does not produce an integer number of segements.
0040 if fix(k) ~= k
0041    warning('The number of segments is not an integer, truncating data.');
0042 end
0043 k = fix(k);
0044 
0045 % fix vectors
0046 x   = x(:);
0047 win = win(:);
0048 
0049 % Initialise the spectral estimates
0050 Sxx = zeros(nfft,1);
0051 
0052 % loop over segments
0053 LminusOverlap = L-noverlap;
0054 xStart = 1:LminusOverlap:k*LminusOverlap;
0055 xEnd   = xStart+L-1;
0056 
0057 for i = 1:k,
0058   %disp(sprintf('* computing segment %d', i));
0059   % data segment, detrend and window
0060   dseg = detrend(x(xStart(i):xEnd(i))) .* win;
0061   %dseg = (x(xStart(i):xEnd(i))) .* win;
0062   
0063   % compute fft
0064   xx   = fft(dseg,nfft);  
0065   Sxx = ((Sxx .* (i-1)) +  xx.*conj(xx)) ./ i;
0066 end
0067 % take only one-sided
0068 Sxx = Sxx(1:nfft/2+1);
0069 
0070 % scale the data as required
0071 [Sxx, info] = scale_xx(Sxx, win, fs, scale, inunits);
0072 
0073 % grow frequency vector
0074 f = linspace(0, fs/2, length(Sxx));
0075 
0076 % complete info
0077 info.nsecs = Nx/fs;
0078 info.navs  = k;
0079 info.nolap = noverlap/fs;
0080 
0081 if nargout == 3
0082   
0083   varargout{1} = f;
0084   varargout{2} = Sxx;
0085   varargout{3} = info;
0086   
0087   
0088 else
0089   disp('usage: [f, pxx, info] = ltpda_psd(x, win, nfft, nolap, fs)');
0090   error('# incorrect usage');
0091 end  
0092 
0093 
0094 %-----------------------------------------------------------------------------------------------
0095 function [yy, info] = scale_xx(xx, win, fs, norm, inunits)
0096 
0097 nfft = length(win);
0098 S1   = sum(win);
0099 S2   = sum(win.^2);
0100 enbw = fs * S2 / (S1*S1);
0101 
0102 switch norm
0103   case 'ASD'
0104     yy = sqrt(xx) * sqrt(2 / (fs*S2));
0105     info.units = [sprintf('%s', inunits) '/\surdHz'];
0106   case 'PSD'
0107     yy = xx * 2 /(fs*S2);
0108     info.units = sprintf('%s^2/Hz', inunits);
0109   case 'AS'
0110     yy = sqrt(xx) * sqrt(2) / S1;
0111     info.units = sprintf('%s', inunits);
0112   case 'PS'
0113     yy = xx * 2 / (S1^2);
0114     info.units = sprintf('%s^2', inunits);
0115   case 'none'
0116     yy = xx;
0117     info.units = '';
0118   otherwise
0119     error('Unknown normalisation');
0120 end
0121 
0122 info.nfft = nfft;
0123 info.enbw = enbw;
0124 info.norm = norm;
0125 
0126 %-----------------------------------------------------------------------------------------------
0127 function [L,noverlap,win,msg] = segment_info(M,win,noverlap)
0128 %SEGMENT_INFO   Determine the information necessary to segment the input data.
0129 %
0130 %   Inputs:
0131 %      M        - An integer containing the length of the data to be segmented
0132 %      WIN      - A scalar or vector containing the length of the window or the window respectively
0133 %                 (Note that the length of the window determines the length of the segments)
0134 %      NOVERLAP - An integer containing the number of samples to overlap (may be empty)
0135 %
0136 %   Outputs:
0137 %      L        - An integer containing the length of the segments
0138 %      NOVERLAP - An integer containing the number of samples to overlap
0139 %      WIN      - A vector containing the window to be applied to each section
0140 %      MSG      - A string containing possible error messages
0141 %
0142 %
0143 %   The key to this function is the following equation:
0144 %
0145 %      K = (M-NOVERLAP)/(L-NOVERLAP)
0146 %
0147 %   where
0148 %
0149 %      K        - Number of segments
0150 %      M        - Length of the input data X
0151 %      NOVERLAP - Desired overlap
0152 %      L        - Length of the segments
0153 %
0154 %   The segmentation of X is based on the fact that we always know M and two of the set
0155 %   {K,NOVERLAP,L}, hence determining the unknown quantity is trivial from the above
0156 %   formula.
0157 
0158 % Initialize outputs
0159 L = [];
0160 msg = '';
0161 
0162 % Check that noverlap is a scalar
0163 if any(size(noverlap) > 1),
0164     msg.identifier = generatemsgid('invalidNoverlap');
0165     msg.message = 'You must specify an integer number of samples to overlap.';
0166     return;
0167 end
0168 
0169 if isempty(win),
0170     % Use 8 sections, determine their length
0171     if isempty(noverlap),
0172         % Use 50% overlap
0173         L = fix(M./4.5);
0174         noverlap = fix(0.5.*L);
0175     else
0176         L = fix((M+7.*noverlap)./8);
0177     end
0178     % Use a default window
0179     win = hamming(L);
0180 else
0181     % Determine the window and its length (equal to the length of the segments)
0182     if ~any(size(win) <= 1) | ischar(win),
0183         msg.identifier = generatemsgid('invalidWindow');
0184         msg.message = 'The WINDOW argument must be a vector or a scalar.';
0185         return
0186     elseif length(win) > 1,
0187         % WIN is a vector
0188         L = length(win);
0189     elseif length(win) == 1,
0190         L = win;
0191         win = hamming(win);
0192     end
0193     if isempty(noverlap),
0194         % Use 50% overlap
0195         noverlap = fix(0.5.*L);
0196     end
0197 end
0198 
0199 % Do some argument validation
0200 if L > M,
0201     msg.identifier = generatemsgid('invalidSegmentLength');
0202     msg.message = 'The length of the segments cannot be greater than the length of the input signal.';
0203     return;
0204 end
0205 
0206 if noverlap >= L,
0207     msg.identifier = generatemsgid('invalidNoverlap');
0208     msg.message = 'The number of samples to overlap must be less than the length of the segments.';
0209     return;
0210 end
0211 
0212 %-----------------------------------------------------------------------------------------------
0213 
0214 % END

Generated on Fri 08-Jun-2007 16:09:11 by m2html © 2003