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