Mercurial > hg > ltpda

comparison m-toolbox/classes/@ao/welchparse.m @ 0:f0afece42f48

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
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 % WELCHPARSE Parser for the PWELCH & SPECTROGRAM functions
2 %
3 % Outputs:
4 % X - First input signal (used when esttype = MS & PSD)
5 % M - An integer containing the length of the data to be segmented
6 % isreal_x - Boolean for input complexity
7 % Y - Second input signal (used when esttype = CPSD, TFE, MSCOHERE)
8 % Ly - Length of second input signal (used when esttype = CPSD, TFE, MSCOHERE)
9 % WIN - A scalar or vector containing the length of the window or the
10 % window respectively (Note that the length of the window determines
11 % the length of the segments)
12 % WINNAME - String with the window name.
13 % WINPARAM - Window parameter.
14 % NOVERLAP - An integer containing the number of samples to overlap (may
15 % be empty)
16 % K - Number of segments
17 % OPTIONS - A structure with the following fields:
18 % OPTIONS.nfft - number of freq. points at which the psd is estimated
19 % OPTIONS.Fs - sampling freq. if any
20 % OPTIONS.range - 'onesided' or 'twosided' psd
21 %
22 % Direct copy from MATLAB
23 %
24 % M Hewitson 08-05-08
25 %
26 % $Id: welchparse.m,v 1.3 2010/07/23 14:53:40 mauro Exp $
27 %
28
29 % Author(s): P. Costa
30 % Copyright 1988-2006 The MathWorks, Inc.
31 % $Revision: 1.3 $ $Date: 2010/07/23 14:53:40 $
32
33 function [x,M,isreal_x,y,Ly,win,winName,winParam,noverlap,k,L,options] = ...
34 welchparse(x,esttype,varargin)
35
36 % Parse input arguments.
37 [x,M,isreal_x,y,Ly,win,winName,winParam,noverlap,opts,errid,errmsg] = ...
38 parse_inputs(x,esttype,varargin{:});
39 if ~isempty(errmsg), error(errid,errmsg); end
40
41 % Obtain the necessary information to segment x and y.
42 [L,noverlap,win,errid,errmsg] = segment_info(M,win,noverlap);
43 if ~isempty(errmsg), error(errid,errmsg); end
44
45 % Parse optional args nfft, fs, and spectrumType.
46 [options,msg] = welch_options(isreal_x,L,opts{:});
47 if ~isempty(msg), error(generatemsgid('InvalidParam'),msg); end
48
49 % Compute the number of segments
50 k = (M-noverlap)./(L-noverlap);
51
52 % Uncomment the following line to produce a warning each time the data
53 % segmentation does not produce an integer number of segements.
54 %if fix(k) ~= k,
55 % warning(generatemsgid('MustBeInteger'),'The number of segments is not an integer, truncating data.');
56 %end
57
58 k = fix(k);
59 end
60
61 %-----------------------------------------------------------------------------------------------
62 function [x,Lx,isreal_x,y,Ly,win,winName,winParam,noverlap,opts,errid,errmsg] = ...
63 parse_inputs(x,esttype,varargin)
64 % Parse the inputs to the welch function.
65
66 % Assign defaults in case of early return.
67 isreal_x = 1;
68 y = [];
69 Ly = 0;
70 is2sig = false;
71 win = [];
72 winName = 'User Defined';
73 winParam = '';
74 noverlap = [];
75 opts = {};
76 errid = '';
77 errmsg = '';
78
79 % Determine if one or two signal vectors was specified.
80 Lx = length(x);
81 if iscell(x),
82 if Lx > 1, % Cell array.
83 y = x{2};
84 is2sig = true;
85 end
86 x = x{1};
87 Lx = length(x);
88 else
89 if ~any(strcmpi(esttype,{'psd','ms'})),
90 errid = generatemsgid('invalidSignalVectors');
91 errmsg = 'You must specify a cell array with two signal vectors to estimate either the cross power spectral density or the transfer function.';
92 return;
93 end
94 end
95
96 x = x(:);
97 isreal_x = isreal(x);
98
99 % Parse 2nd input signal vector.
100 if is2sig,
101 y = y(:);
102 isreal_x = isreal(y) && isreal_x;
103 Ly = length(y);
104 if Ly ~= Lx,
105 errid = generatemsgid('invalidSignalVectors');
106 errmsg = 'The length of the two input vectors must be equal to calculate the cross spectral density.';
107 return;
108 end
109 end
110
111 % Parse window and overlap, and cache remaining inputs.
112 lenargin = length(varargin);
113 if lenargin >= 1,
114 win = varargin{1};
115 if lenargin >= 2,
116 noverlap = varargin{2};
117
118 % Cache optional args nfft, fs, and spectrumType.
119 if lenargin >= 3, opts = varargin(3:end); end
120 end
121 end
122
123 if isempty(win) | isscalar(win),
124 winName = 'hamming';
125 winParam = 'symmetric';
126 end
127
128 end
129
130 %-----------------------------------------------------------------------------------------------
131 %SEGMENT_INFO Determine the information necessary to segment the input data.
132 %
133 % Inputs:
134 % M - An integer containing the length of the data to be segmented
135 % WIN - A scalar or vector containing the length of the window or the window respectively
136 % (Note that the length of the window determines the length of the segments)
137 % NOVERLAP - An integer containing the number of samples to overlap (may be empty)
138 %
139 % Outputs:
140 % L - An integer containing the length of the segments
141 % NOVERLAP - An integer containing the number of samples to overlap
142 % WIN - A vector containing the window to be applied to each section
143 % ERRID - A string containing possible error identifier
144 % ERRMSG - A string containing possible error messages
145 %
146 %
147 % The key to this function is the following equation:
148 %
149 % K = (M-NOVERLAP)/(L-NOVERLAP)
150 %
151 % where
152 %
153 % K - Number of segments
154 % M - Length of the input data X
155 % NOVERLAP - Desired overlap
156 % L - Length of the segments
157 %
158 % The segmentation of X is based on the fact that we always know M and two of the set
159 % {K,NOVERLAP,L}, hence determining the unknown quantity is trivial from the above
160 % formula.
161
162 function [L,noverlap,win,errid,errmsg] = segment_info(M,win,noverlap)
163
164 % Initialize outputs
165 L = [];
166 errid = '';
167 errmsg = '';
168
169 % Check that noverlap is a scalar
170 if any(size(noverlap) > 1),
171 errid = generatemsgid('invalidNoverlap');
172 errmsg = 'You must specify an integer number of samples to overlap.';
173 return;
174 end
175
176 if isempty(win),
177 % Use 8 sections, determine their length
178 if isempty(noverlap),
179 % Use 50% overlap
180 L = fix(M./4.5);
181 noverlap = fix(0.5.*L);
182 else
183 L = fix((M+7.*noverlap)./8);
184 end
185 % Use a default window
186 win = hamming(L);
187 else
188 % Determine the window and its length (equal to the length of the segments)
189 if ~any(size(win) <= 1) | ischar(win),
190 errid = generatemsgid('invalidWindow');
191 errmsg = 'The WINDOW argument must be a vector or a scalar.';
192 return
193 elseif length(win) > 1,
194 % WIN is a vector
195 L = length(win);
196 elseif length(win) == 1,
197 L = win;
198 win = hamming(win);
199 end
200 if isempty(noverlap),
201 % Use 50% overlap
202 noverlap = fix(0.5.*L);
203 end
204 end
205
206 % Do some argument validation
207 if L > M,
208 errid = generatemsgid('invalidSegmentLength');
209 errmsg = 'The length of the segments cannot be greater than the length of the input signal.';
210 return;
211 end
212
213 if noverlap >= L,
214 errid = generatemsgid('invalidNoverlap');
215 errmsg = 'The number of samples to overlap must be less than the length of the segments.';
216 return;
217 end
218 end
219
220 %------------------------------------------------------------------------------
221 %WELCH_OPTIONS Parse the optional inputs to the PWELCH function.
222 % WELCH_OPTIONS returns a structure, OPTIONS, with following fields:
223 %
224 % options.nfft - number of freq. points at which the psd is estimated
225 % options.Fs - sampling freq. if any
226 % options.range - 'onesided' or 'twosided' psd
227
228 function [options,msg] = welch_options(isreal_x,N,varargin)
229
230 % Generate defaults
231 options.nfft = max(256,2^nextpow2(N));
232 options.Fs = []; % Work in rad/sample
233
234 % Determine if frequency vector specified
235 freqVecSpec = false;
236 if (length(varargin) > 0 && length(varargin{1}) > 1)
237 freqVecSpec = true;
238 end
239
240 if isreal_x && ~freqVecSpec,
241 options.range = 'onesided';
242 else
243 options.range = 'twosided';
244 end
245 msg = '';
246
247 if any(strcmp(varargin, 'whole'))
248 warning(generatemsgid('invalidRange'), '''whole'' is not a valid range, use ''twosided'' instead.');
249 elseif any(strcmp(varargin, 'half'))
250 warning(generatemsgid('invalidRange'), '''half'' is not a valid range, use ''onesided'' instead.');
251 end
252
253 [options,msg] = psdoptions(isreal_x,options,varargin{:});
254
255 end
256