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 $
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