1 %% Copyright (C) 2006 Peter V. Lanspeary <pvl@mecheng.adelaide.edu.au>
3 %% This program is free software; you can redistribute it and/or modify it under
4 %% the terms of the GNU General Public License as published by the Free Software
5 %% Foundation; either version 3 of the License, or (at your option) any later
8 %% This program is distributed in the hope that it will be useful, but WITHOUT
9 %% ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
10 %% FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
13 %% You should have received a copy of the GNU General Public License along with
14 %% this program; if not, see <http://www.gnu.org/licenses/>.
17 %% [spectra,freq] = pwelch(x,window,overlap,Nfft,Fs,
18 %% range,plot_type,detrend,sloppy)
19 %% Estimate power spectral density of data "x" by the Welch (1967)
20 %% periodogram/FFT method. All arguments except "x" are optional.
21 %% The data is divided into segments. If "window" is a vector, each
22 %% segment has the same length as "window" and is multiplied by "window"
23 %% before (optional) zero-padding and calculation of its periodogram. If
24 %% "window" is a scalar, each segment has a length of "window" and a
25 %% Hamming window is used.
26 %% The spectral density is the mean of the periodograms, scaled so that
27 %% area under the spectrum is the same as the mean square of the
28 %% data. This equivalence is supposed to be exact, but in practice there
29 %% is a mismatch of up to 0.5% when comparing area under a periodogram
30 %% with the mean square of the data.
32 %% [spectra,freq] = pwelch(x,y,window,overlap,Nfft,Fs,
33 %% range,plot_type,detrend,sloppy,results)
34 %% Two-channel spectrum analyser. Estimate power spectral density, cross-
35 %% spectral density, transfer function and/or coherence functions of time-
36 %% series input data "x" and output data "y" by the Welch (1967)
37 %% periodogram/FFT method.
38 %% pwelch treats the second argument as "y" if there is a control-string
39 %% argument "cross", "trans", "coher" or "ypower"; "power" does not force
40 %% the 2nd argument to be treated as "y". All other arguments are
41 %% optional. All spectra are returned in matrix "spectra".
43 %% [spectra,Pxx_ci,freq] = pwelch(x,window,overlap,Nfft,Fs,conf,
44 %% range,plot_type,detrend,sloppy)
45 %% [spectra,Pxx_ci,freq] = pwelch(x,y,window,overlap,Nfft,Fs,conf,
46 %% range,plot_type,detrend,sloppy,results)
47 %% Estimates confidence intervals for the spectral density.
48 %% See Hint (7) below for compatibility options. Confidence level "conf"
49 %% is the 6th or 7th numeric argument. If "results" control-string
50 %% arguments are used, one of them must be "power" when the "conf"
51 %% argument is present; pwelch can estimate confidence intervals only for
52 %% the power spectrum of the "x" data. It does not know how to estimate
53 %% confidence intervals of the cross-power spectrum, transfer function or
54 %% coherence; if you can suggest a good method, please send a bug report.
57 %% All but the first argument are optional and may be empty, except that
58 %% the "results" argument may require the second argument to be "y".
60 %% x %% [non-empty vector] system-input time-series data
61 %% y %% [non-empty vector] system-output time-series data
63 %% window %% [real vector] of window-function values between 0 and 1; the
64 %% %% data segment has the same length as the window.
65 %% %% Default window shape is Hamming.
66 %% %% [integer scalar] length of each data segment. The default
67 %% %% value is window=sqrt(length(x)) rounded up to the
68 %% %% nearest integer power of 2; see 'sloppy' argument.
70 %% overlap %% [real scalar] segment overlap expressed as a multiple of
71 %% %% window or segment length. 0 <= overlap < 1,
72 %% %% The default is overlap=0.5 .
74 %% Nfft %% [integer scalar] Length of FFT. The default is the length
75 %% %% of the "window" vector or has the same value as the
76 %% %% scalar "window" argument. If Nfft is larger than the
77 %% %% segment length, "seg_len", the data segment is padded
78 %% %% with "Nfft-seg_len" zeros. The default is no padding.
79 %% %% Nfft values smaller than the length of the data
80 %% %% segment (or window) are ignored silently.
82 %% Fs %% [real scalar] sampling frequency (Hertz); default=1.0
84 %% conf %% [real scalar] confidence level between 0 and 1. Confidence
85 %% %% intervals of the spectral density are estimated from
86 %% %% scatter in the periodograms and are returned as Pxx_ci.
87 %% %% Pxx_ci(:,1) is the lower bound of the confidence
88 %% %% interval and Pxx_ci(:,2) is the upper bound. If there
89 %% %% are three return values, or conf is an empty matrix,
90 %% %% confidence intervals are calculated for conf=0.95 .
91 %% %% If conf is zero or is not given, confidence intervals
92 %% %% are not calculated. Confidence intervals can be
93 %% %% obtained only for the power spectral density of x;
96 %% CONTROL-STRING ARGUMENTS -- each of these arguments is a character string.
97 %% Control-string arguments must be after the other arguments but can be in
100 %% range %% 'half', 'onesided' : frequency range of the spectrum is
101 %% %% zero up to but not including Fs/2. Power from
102 %% %% negative frequencies is added to the positive side of
103 %% %% the spectrum, but not at zero or Nyquist (Fs/2)
104 %% %% frequencies. This keeps power equal in time and
105 %% %% spectral domains. See reference [2].
106 %% %% 'whole', 'twosided' : frequency range of the spectrum is
107 %% %% -Fs/2 to Fs/2, with negative frequencies
108 %% %% stored in "wrap around" order after the positive
109 %% %% frequencies; e.g. frequencies for a 10-point 'twosided'
110 %% %% spectrum are 0 0.1 0.2 0.3 0.4 0.5 -0.4 -0.3 -0.2 -0.1
111 %% %% 'shift', 'centerdc' : same as 'whole' but with the first half
112 %% %% of the spectrum swapped with second half to put the
113 %% %% zero-frequency value in the middle. (See "help
115 %% %% If data (x and y) are real, the default range is 'half',
116 %% %% otherwise default range is 'whole'.
118 %% plot_type %% 'plot', 'semilogx', 'semilogy', 'loglog', 'squared' or 'db':
119 %% %% specifies the type of plot. The default is 'plot', which
120 %% %% means linear-linear axes. 'squared' is the same as 'plot'.
121 %% %% 'dB' plots "10*log10(psd)". This argument is ignored and a
122 %% %% spectrum is not plotted if the caller requires a returned
125 %% detrend %% 'no-strip', 'none' -- do NOT remove mean value from the data
126 %% %% 'short', 'mean' -- remove the mean value of each segment from
127 %% %% each segment of the data.
128 %% %% 'linear', -- remove linear trend from each segment of
130 %% %% 'long-mean' -- remove the mean value from the data before
131 %% %% splitting it into segments. This is the default.
133 %% sloppy %% 'sloppy': FFT length is rounded up to the nearest integer
134 %% %% power of 2 by zero padding. FFT length is adjusted
135 %% %% after addition of padding by explicit Nfft argument.
136 %% %% The default is to use exactly the FFT and window/
138 %% %% segment lengths specified in argument list.
140 %% results %% specifies what results to return (in the order specified
141 %% %% and as many as desired).
142 %% %% 'power' calculate power spectral density of "x"
143 %% %% 'cross' calculate cross spectral density of "x" and "y"
144 %% %% 'trans' calculate transfer function of a system with
145 %% %% input "x" and output "y"
146 %% %% 'coher' calculate coherence function of "x" and "y"
147 %% %% 'ypower' calculate power spectral density of "y"
148 %% %% The default is 'power', with argument "y" omitted.
151 %% If return values are not required by the caller, the results are
152 %% plotted and nothing is returned.
154 %% spectra %% [real-or-complex matrix] columns of the matrix contain results
155 %% %% in the same order as specified by "results" arguments.
156 %% %% Each column contains one of the result vectors.
158 %% Pxx_ci %% [real matrix] estimate of confidence interval for power
159 %% %% spectral density of x. First column is the lower
160 %% %% bound. Second column is the upper bound.
162 %% freq %% [real column vector] frequency values
166 %% if you don't want to use an optional argument you can leave it empty
167 %% by writing its value as [].
169 %% The profusion of arguments may make pwelch difficult to use, and an
170 %% unskilled user can easily produce a meaningless result or can easily
171 %% mis-interpret the result. With real data "x" and sampling frequency
172 %% "Fs", the easiest and best way for a beginner to use pwelch is
173 %% probably "pwelch(x,[],[],[],Fs)". Use the "window" argument to
174 %% control the length of the spectrum vector. For real data and integer
175 %% scalar M, "pwelch(x,2*M,[],[],Fs)" gives an M+1 point spectrum.
176 %% Run "demo pwelch" (octave only).
177 %% 3) WINDOWING FUNCTIONS:
178 %% Without a window function, sharp spectral peaks can have strong
179 %% sidelobes because the FFT of a data in a segment is in effect convolved
180 %% with a rectangular window. A window function which tapers off
181 %% (gradually) at the ends produces much weaker sidelobes in the FFT.
182 %% Hann (hanning), hamming, bartlett, blackman, flattopwin etc are
183 %% available as separate Matlab/sigproc or Octave functions. The sidelobes
184 %% of the Hann window have a roll-off rate of 60dB/decade of frequency.
185 %% The first sidelobe of the Hamming window is suppressed and is about 12dB
186 %% lower than the first Hann sidelobe, but the roll-off rate is only
187 %% 20dB/decade. You can inspect the FFT of a Hann window by plotting
188 %% "abs(fft(postpad(hanning(256),4096,0)))".
189 %% The default window is Hamming.
191 %% Zero-padding reduces the frequency step in the
192 %% spectrum, and produces an artificially smoothed spectrum. For example,
193 %% "Nfft=2*length(window)" gives twice as many frequency values, but
194 %% adjacent PSD (power spectral density) values are not independent;
195 %% adjacent PSD values are independent if "Nfft=length(window)", which is
196 %% the default value of Nfft.
197 %% 5) REMOVING MEAN FROM SIGNAL:
198 %% If the mean is not removed from the signal there is a large spectral
199 %% peak at zero frequency and the sidelobes of this peak are likely to
200 %% swamp the rest of the spectrum. For this reason, the default behaviour
201 %% is to remove the mean. However, the matlab pwelch does not do this.
202 %% 6) WARNING ON CONFIDENCE INTERVALS
203 %% Confidence intervals are obtained by measuring the sample variance of
204 %% the periodograms and assuming that the periodograms have a Gaussian
205 %% probability distribution. This assumption is not accurate. If, for
206 %% example, the data (x) is Gaussian, the periodogram has a Rayleigh
207 %% distribution. However, the confidence intervals may still be useful.
209 %% 7) COMPATIBILITY WITH Matlab R11, R12, etc
210 %% When used without the second data (y) argument, arguments are compatible
211 %% with the pwelch of Matlab R12, R13, R14, 2006a and 2006b except that
212 %% 1) overlap is expressed as a multiple of window length ---
213 %% effect of overlap scales with window length
214 %% 2) default values of length(window), Nfft and Fs are more sensible, and
215 %% 3) Goertzel algorithm is not available so Nfft cannot be an array of
216 %% frequencies as in Matlab 2006b.
217 %% Pwelch has four persistent Matlab-compatibility levels. Calling pwelch
218 %% with an empty first argument sets the order of arguments and defaults
219 %% specified above in the USAGE and ARGUMENTS section of this documentation.
220 %% prev_compat=pwelch([]);
221 %% [Pxx,f]=pwelch(x,window,overlap,Nfft,Fs,conf,...);
222 %% Calling pwelch with a single string argument (as described below) gives
223 %% compatibility with Matlab R11 or R12, or the R14 spectrum.welch
224 %% defaults. The returned value is the PREVIOUS compatibility string.
226 %% Matlab R11: For compatibility with the Matlab R11 pwelch:
227 %% prev_compat=pwelch('R11-');
228 %% [Pxx,f]=pwelch(x,Nfft,Fs,window,overlap,conf,range,units);
229 %% %% units of overlap are "number of samples"
230 %% %% defaults: Nfft=min(length(x),256), Fs=2*pi, length(window)=Nfft,
231 %% %% window=Hanning, do not detrend,
232 %% %% N.B. "Sloppy" is not available.
234 %% Matlab R12: For compatibility with Matlab R12 to 2006a pwelch:
235 %% prev_compat=pwelch('R12+');
236 %% [Pxx,f]=pwelch(x,window,overlap,nfft,Fs,...);
237 %% %% units of overlap are "number of samples"
238 %% %% defaults: length(window)==length(x)/8, window=Hamming,
239 %% %% Nfft=max(256,NextPow2), Fs=2*pi, do not detrend
240 %% %% NextPow2 is the next power of 2 greater than or equal to the
241 %% %% window length. "Sloppy", "conf" are not available. Default
242 %% %% window length gives very poor amplitude resolution.
244 %% To adopt defaults of the Matlab R14 "spectrum.welch" spectrum object
245 %% associated "psd" method.
246 %% prev_compat=pwelch('psd');
247 %% [Pxx,f] = pwelch(x,window,overlap,Nfft,Fs,conf,...);
248 %% %% overlap is expressed as a percentage of window length,
249 %% %% defaults: length(window)==64, Nfft=max(256,NextPow2), Fs=2*pi
251 %% %% NextPow2 is the next power of 2 greater than or equal to the
252 %% %% window length. "Sloppy" is not available.
253 %% %% Default window length gives coarse frequency resolution.
257 %% [1] Peter D. Welch (June 1967):
258 %% "The use of fast Fourier transform for the estimation of power spectra:
259 %% a method based on time averaging over short, modified periodograms."
260 %% IEEE Transactions on Audio Electroacoustics, Vol AU-15(6), pp 70-73
262 %% [2] William H. Press and Saul A. Teukolsky and William T. Vetterling and
263 %% Brian P. Flannery",
264 %% "Numerical recipes in C, The art of scientific computing", 2nd edition,
265 %% Cambridge University Press, 2002 --- Section 13.7.
266 %% [3] Paul Kienzle (1999-2001): "pwelch", http://octave.sourceforge.net/
268 function [varargout] = pwelch(x,varargin)
270 %% COMPATIBILITY LEVEL
271 %% Argument positions and defaults depend on compatibility level selected
272 %% by calling pwelch without arguments or with a single string argument.
273 %% native: compatib=1; prev_compat=pwelch(); prev_compat=pwelch([]);
274 %% matlab R11: compatib=2; prev_compat=pwelch('R11-');
275 %% matlab R12: compatib=3; prev_compat=pwelch('R12+');
276 %% spectrum.welch defaults: compatib=4; prev_compat=pwelch('psd');
277 %% In each case, the returned value is the PREVIOUS compatibility string.
279 compat_str = {[]; 'R11-'; 'R12+'; 'psd'};
281 if ( isempty(compatib) || compatib<=0 || compatib>4 )
282 %% legal values are 1, 2, 3, 4
287 error( 'pwelch: Need at least 1 arg. Use "help pwelch".' );
288 elseif ( nargin==1 && (ischar(x) || isempty(x)) )
289 varargout{1} = compat_str{compatib};
290 if ( isempty(x) ) % native
292 elseif ( strcmp(x,'R11-') )
294 elseif ( strcmp(x,'R12+') )
296 elseif ( strcmp(x,'psd') )
299 error( 'pwelch: compatibility arg must be empty, R11-, R12+ or psd' );
303 %% Check fixed argument
304 elseif ( isempty(x) || ~isvector(x) )
305 error( 'pwelch: arg 1 (x) must be vector.' );
307 %% force x to be COLUMN vector
312 %% Look through all args to check if cross PSD, transfer function or
313 %% coherence is required. If yes, the second arg is data vector "y".
316 nvarargin = length(varargin);
318 arg = varargin{iarg};
319 if ( ~isempty(arg) && ischar(arg) && ...
320 ( strcmp(arg,'cross') || strcmp(arg,'trans') || ...
321 strcmp(arg,'coher') || strcmp(arg,'ypower') ))
322 %% OK. Need "y". Grab it from 2nd arg.
324 if ( nargin<2 || isempty(arg) || ~isvector(arg) || length(arg)~=x_len )
325 error( 'pwelch: arg 2 (y) must be vector, same length as x.' );
327 %% force COLUMN vector
335 %% To select default argument values, "compatib" is used as an array index.
336 %% Index values are 1=native, 2=R11, 3=R12, 4=spectrum.welch
338 %% argument positions:
339 %% arg_posn = varargin index of window, overlap, Nfft, Fs and conf
340 %% args respectively, a value of zero ==>> arg does not exist
341 arg_posn = [1 2 3 4 5; %% native
342 3 4 1 2 5; %% Matlab R11- pwelch
343 1 2 3 4 0; %% Matlab R12+ pwelch
344 1 2 3 4 5]; %% spectrum.welch defaults
345 arg_posn = arg_posn(compatib,:) + arg2_is_y;
347 %% SPECIFY SOME DEFAULT VALUES for (not all) optional arguments
348 %% Use compatib as array index.
349 %% Fs = sampling frequency
350 Fs = [ 1.0 2*pi 2*pi 2*pi ];
352 %% plot_type: 1='plot'|'squared'; 5='db'|'dB'
353 plot_type = [ 1 5 5 5 ];
354 plot_type = plot_type(compatib);
355 %% rm_mean: 3='long-mean'; 0='no-strip'|'none'
356 rm_mean = [ 3 0 0 0 ];
357 rm_mean = rm_mean(compatib);
358 %% use max_overlap=x_len-1 because seg_len is not available yet
359 %% units of overlap are different for each version:
360 %% fraction, samples, or percent
361 max_overlap = [ 0.95 x_len-1 x_len-1 95];
362 max_overlap = max_overlap(compatib);
363 %% default confidence interval
364 %% if there are more than 2 return values and if there is a "conf" arg
365 conf = 0.95 * (nargout>2) * (arg_posn(5)>0);
367 is_win = 0; % =0 means valid window arg is not provided yet
368 Nfft = []; % default depends on segment length
369 overlap = []; % WARNING: units can be #samples, fraction or percentage
370 range = ~isreal(x) || ( arg2_is_y && ~isreal(y) );
379 %% DECODE AND CHECK OPTIONAL ARGUMENTS
380 end_numeric_args = 0;
381 for iarg = 1+arg2_is_y:nvarargin
382 arg = varargin{iarg};
384 %% first string arg ==> no more numeric args
385 %% non-string args cannot follow a string arg
386 end_numeric_args = 1;
388 %% decode control-string arguments
389 if ( strcmp(arg,'sloppy') )
390 is_sloppy = ~is_win || is_win==1;
391 elseif ( strcmp(arg,'plot') || strcmp(arg,'squared') )
393 elseif ( strcmp(arg,'semilogx') )
395 elseif ( strcmp(arg,'semilogy') )
397 elseif ( strcmp(arg,'loglog') )
399 elseif ( strcmp(arg,'db') || strcmp(arg,'dB') )
401 elseif ( strcmp(arg,'half') || strcmp(arg,'onesided') )
403 elseif ( strcmp(arg,'whole') || strcmp(arg,'twosided') )
405 elseif ( strcmp(arg,'shift') || strcmp(arg,'centerdc') )
407 elseif ( strcmp(arg,'long-mean') )
409 elseif ( strcmp(arg,'linear') )
411 elseif ( strcmp(arg,'short') || strcmp(arg,'mean') )
413 elseif ( strcmp(arg,'no-strip') || strcmp(arg,'none') )
415 elseif ( strcmp(arg, 'power' ) )
417 n_results = n_results+1;
418 do_power = n_results;
420 elseif ( strcmp(arg, 'cross' ) )
422 n_results = n_results+1;
423 do_cross = n_results;
425 elseif ( strcmp(arg, 'trans' ) )
427 n_results = n_results+1;
428 do_trans = n_results;
430 elseif ( strcmp(arg, 'coher' ) )
432 n_results = n_results+1;
433 do_coher = n_results;
435 elseif ( strcmp(arg, 'ypower' ) )
437 n_results = n_results+1;
438 do_ypower = n_results;
441 error( 'pwelch: string arg %d illegal value: %s', iarg+1, arg );
443 %% end of processing string args
445 elseif ( end_numeric_args )
447 %% found non-string arg after a string arg ... oops
448 error( 'pwelch: control arg must be string' );
451 %% first 4 optional arguments are numeric -- in fixed order
453 %% deal with "Fs" and "conf" first because empty arg is a special default
454 %% -- "Fs" arg -- sampling frequency
455 elseif ( iarg == arg_posn(4) )
458 elseif ( ~isscalar(arg) || ~isreal(arg) || arg<0 )
459 error( 'pwelch: arg %d (Fs) must be real scalar >0', iarg+1 );
464 %% -- "conf" arg -- confidence level
465 %% guard against the "it cannot happen" iarg==0
466 elseif ( arg_posn(5) && iarg == arg_posn(5) )
469 elseif ( ~isscalar(arg) || ~isreal(arg) || arg < 0.0 || arg >= 1.0 )
470 error( 'pwelch: arg %d (conf) must be real scalar, >=0, <1',iarg+1 );
475 %% skip all empty args from this point onward
476 elseif ( isempty(arg) )
479 %% -- "window" arg -- window function
480 elseif ( iarg == arg_posn(1) )
483 elseif ( isvector(arg) )
484 is_win = length(arg);
485 if ( size(arg,2)>1 ) %% vector must be COLUMN vector
492 error( 'pwelch: arg %d (window) must be scalar or vector', iarg+1 );
493 elseif ( is_win==1 && ( ~isreal(arg) || fix(arg)~=arg || arg<=3 ) )
494 error( 'pwelch: arg %d (window) must be integer >3', iarg+1 );
495 elseif ( is_win>1 && ( ~isreal(arg) || any(arg<0) ) )
496 error( 'pwelch: arg %d (window) vector must be real and >=0',iarg+1);
501 %% -- "overlap" arg -- segment overlap
502 elseif ( iarg == arg_posn(2) )
503 if (~isscalar(arg) || ~isreal(arg) || arg<0 || arg>max_overlap )
504 error( 'pwelch: arg %d (overlap) must be real from 0 to %f', ...
505 iarg+1, max_overlap );
509 %% -- "Nfft" arg -- FFT length
510 elseif ( iarg == arg_posn(3) )
511 if ( ~isscalar(arg) || ~isreal(arg) || fix(arg)~=arg || arg<0 )
512 error( 'pwelch: arg %d (Nfft) must be integer >=0', iarg+1 );
517 error( 'pwelch: arg %d must be string', iarg+1 );
520 if ( conf>0 && (n_results && ~do_power ) )
521 error('pwelch: can give confidence interval for x power spectrum only' );
524 %% end DECODE AND CHECK OPTIONAL ARGUMENTS.
526 %% SETUP REMAINING PARAMETERS
527 %% default action is to calculate power spectrum only
532 need_Pxx = do_power || do_trans || do_coher;
533 need_Pxy = do_cross || do_trans || do_coher;
534 need_Pyy = do_coher || do_ypower;
536 nearly_one = 0.99999999999;
538 %% compatibility-options
539 %% provides exact compatibility with Matlab R11 or R12
541 %% Matlab R11 compatibility
544 Nfft = min( 256, x_len );
547 seg_len = min( length(window), Nfft );
548 window = window(1:seg_len);
551 %% window arg is scalar
556 %% make Hann window (don't depend on sigproc)
558 window = 0.5 - 0.5 * cos( (2*pi/xx)*[0:xx].' );
561 %% Matlab R12 compatibility
562 elseif ( compatib==3 )
564 %% window arg provides window function
565 seg_len = length(window);
567 %% window arg does not provide window function; use Hamming
569 %% window arg is scalar
572 %% window arg not available; use R12 default, 8 windows
573 %% ignore overlap arg; use overlap=50% -- only choice that makes sense
574 %% this is the magic formula for 8 segments with 50% overlap
575 seg_len = fix( (x_len-3)*2/9 );
577 %% make Hamming window (don't depend on sigproc)
579 window = 0.54 - 0.46 * cos( (2*pi/xx)*[0:xx].' );
582 Nfft = max( 256, 2^ceil(log(seg_len)*nearly_one/log_two) );
585 %% Matlab R14 psd(spectrum.welch) defaults
586 elseif ( compatib==4 )
588 %% window arg provides window function
589 seg_len = length(window);
591 %% window arg does not provide window function; use Hamming
593 %% window arg is scalar
596 %% window arg not available; use default seg_len = 64
599 %% make Hamming window (don't depend on sigproc)
601 window = 0.54 - 0.46 * cos( (2*pi/xx)*[0:xx].' );
603 %% Now we know segment length,
604 %% so we can set default overlap as number of samples
605 if ( ~isempty(overlap) )
606 overlap = fix(seg_len * overlap / 100 );
609 Nfft = max( 256, 2^ceil(log(seg_len)*nearly_one/log_two) );
612 %% default compatibility level
613 else %if ( compatib==1 )
614 %% calculate/adjust segment length, window function
616 %% window arg provides window function
617 seg_len = length(window);
619 %% window arg does not provide window function; use Hamming
620 if ( is_win ) %% window arg is scalar
623 %% window arg not available; use default length:
624 %% = sqrt(length(x)) rounded up to nearest integer power of 2
625 if ( isempty(overlap) )
628 seg_len = 2 ^ ceil( log(sqrt(x_len/(1-overlap)))*nearly_one/log_two );
630 %% make Hamming window (don't depend on sigproc)
632 window = 0.54 - 0.46 * cos( (2*pi/xx)*[0:xx].' );
634 %% Now we know segment length,
635 %% so we can set default overlap as number of samples
636 if ( ~isempty(overlap) )
637 overlap = fix(seg_len * overlap);
640 %% calculate FFT length
645 Nfft = 2 ^ ceil( log(Nfft) * nearly_one / log_two );
648 %% end of compatibility options
650 %% minimum FFT length is seg_len
651 Nfft = max( Nfft, seg_len );
652 %% Mean square of window is required for normalising PSD amplitude.
653 win_meansq = (window.' * window) / seg_len;
655 %% Set default or check overlap.
656 if ( isempty(overlap) )
657 overlap = fix(seg_len /2);
658 elseif ( overlap >= seg_len )
659 error( 'pwelch: arg (overlap=%d) too big. Must be <length(window)=%d',...
663 %% Pad data with zeros if shorter than segment. This should not happen.
664 if ( x_len < seg_len )
665 x = [x; zeros(seg_len-x_len,1)];
667 y = [y; zeros(seg_len-x_len,1)];
671 %% end SETUP REMAINING PARAMETERS
675 %% Remove mean from the data
677 n_ffts = max( 0, fix( (x_len-seg_len)/(seg_len-overlap) ) ) + 1;
678 x_len = min( x_len, (seg_len-overlap)*(n_ffts-1)+seg_len );
679 if ( need_Pxx || need_Pxy )
680 x = x - sum( x(1:x_len) ) / x_len;
682 if ( arg2_is_y || need_Pxy)
683 y = y - sum( y(1:x_len) ) / x_len;
687 %% Calculate and accumulate periodograms
688 %% xx and yy are padded data segments
689 %% Pxx, Pyy, Pyy are periodogram sums, Vxx is for confidence interval
701 for start_seg = [1:seg_len-overlap:x_len-seg_len+1]
702 end_seg = start_seg+seg_len-1;
703 %% Don't truncate/remove the zero padding in xx and yy
704 if ( need_Pxx || need_Pxy )
705 if ( rm_mean==1 ) % remove mean from segment
706 xx(1:seg_len) = window .* ( ...
707 x(start_seg:end_seg) - sum(x(start_seg:end_seg)) / seg_len);
708 elseif ( rm_mean == 2 ) % remove linear trend from segment
709 xx(1:seg_len) = window .* detrend( x(start_seg:end_seg) );
710 else % rm_mean==0 or 3
711 xx(1:seg_len) = window .* x(start_seg:end_seg);
715 if ( need_Pxy || need_Pyy )
716 if ( rm_mean==1 ) % remove mean from segment
717 yy(1:seg_len) = window .* ( ...
718 y(start_seg:end_seg) - sum(y(start_seg:end_seg)) / seg_len);
719 elseif ( rm_mean == 2 ) % remove linear trend from segment
720 yy(1:seg_len) = window .* detrend( y(start_seg:end_seg) );
721 else % rm_mean==0 or 3
722 yy(1:seg_len) = window .* y(start_seg:end_seg);
727 %% force Pxx to be real; pgram = periodogram
728 pgram = real(fft_x .* conj(fft_x));
730 %% sum of squared periodograms is required for confidence interval
732 Vxx = Vxx + pgram .^2;
736 %% Pxy (cross power spectrum) is complex. Do not force to be real.
737 Pxy = Pxy + fft_y .* conj(fft_x);
740 %% force Pyy to be real
741 Pyy = Pyy + real(fft_y .* conj(fft_y));
746 %% Calculate confidence interval
747 %% -- incorrectly assumes that the periodogram has Gaussian probability
748 %% distribution (actually, it has a single-sided (e.g. exponential)
750 %% Sample variance of periodograms is (Vxx-Pxx.^2/n_ffts)/(n_ffts-1).
751 %% This method of calculating variance is more susceptible to round-off
752 %% error, but is quicker, and for double-precision arithmetic and the
753 %% inherently noisy periodogram (variance==mean^2), it should be OK.
754 if ( conf>0 && need_Pxx )
758 %% Should use student distribution here (for unknown variance), but tinv
759 %% is not a core Matlab function (is in statistics toolbox. Grrr)
760 Vxx = (erfinv(conf)*sqrt(2*n_ffts/(n_ffts-1))) * sqrt(Vxx-Pxx.^2/n_ffts);
764 %% Convert two-sided spectra to one-sided spectra (if range == 0).
765 %% For one-sided spectra, contributions from negative frequencies are added
766 %% to the positive side of the spectrum -- but not at zero or Nyquist
767 %% (half sampling) frequencies. This keeps power equal in time and spectral
768 %% domains, as required by Parseval theorem.
771 if ( ~ rem(Nfft,2) ) %% one-sided, Nfft is even
774 Pxx = Pxx(1:psd_len) + [0; Pxx(Nfft:-1:psd_len+1); 0];
776 Vxx = Vxx(1:psd_len) + [0; Vxx(Nfft:-1:psd_len+1); 0];
780 Pxy = Pxy(1:psd_len) + conj([0; Pxy(Nfft:-1:psd_len+1); 0]);
783 Pyy = Pyy(1:psd_len) + [0; Pyy(Nfft:-1:psd_len+1); 0];
785 else %% one-sided, Nfft is odd
786 psd_len = (Nfft+1)/2;
788 Pxx = Pxx(1:psd_len) + [0; Pxx(Nfft:-1:psd_len+1)];
790 Vxx = Vxx(1:psd_len) + [0; Vxx(Nfft:-1:psd_len+1)];
794 Pxy = Pxy(1:psd_len) + conj([0; Pxy(Nfft:-1:psd_len+1)]);
797 Pyy = Pyy(1:psd_len) + [0; Pyy(Nfft:-1:psd_len+1)];
800 else %% two-sided (and shifted)
803 %% end MAIN CALCULATIONS
805 %% SCALING AND OUTPUT
806 %% Put all results in matrix, one row per spectrum
807 %% Pxx, Pxy, Pyy are sums of periodograms, so "n_ffts"
808 %% in the scale factor converts them into averages
809 spectra = zeros(psd_len,n_results);
810 spect_type = zeros(n_results,1);
811 scale = n_ffts * seg_len * Fs * win_meansq;
813 spectra(:,do_power) = Pxx / scale;
814 spect_type(do_power) = 1;
816 Vxx = [Pxx-Vxx Pxx+Vxx]/scale;
820 spectra(:,do_cross) = Pxy / scale;
821 spect_type(do_cross) = 2;
824 spectra(:,do_trans) = Pxy ./ Pxx;
825 spect_type(do_trans) = 3;
828 %% force coherence to be real
829 spectra(:,do_coher) = real(Pxy .* conj(Pxy)) ./ Pxx ./ Pyy;
830 spect_type(do_coher) = 4;
833 spectra(:,do_ypower) = Pyy / scale;
834 spect_type(do_ypower) = 5;
836 freq = [0:psd_len-1].' * ( Fs / Nfft );
838 %% range='shift': Shift zero-frequency to the middle
840 len2 = fix((Nfft+1)/2);
841 spectra = [ spectra(len2+1:Nfft,:); spectra(1:len2,:)];
842 freq = [ freq(len2+1:Nfft)-Fs; freq(1:len2)];
844 Vxx = [ Vxx(len2+1:Nfft,:); Vxx(1:len2,:)];
848 %% RETURN RESULTS or PLOT
849 if ( nargout>=2 && conf>0 )
852 if ( nargout>=(2+(conf>0)) )
853 %% frequency is 2nd or 3rd return value,
854 %% depends on if 2nd is confidence interval
855 varargout{2+(conf>0)} = freq;
858 varargout{1} = spectra;
861 %% Plot the spectra if there are no return variables.
862 plot_title=['power spectrum x ';
866 'power spectrum y ' ];
867 for ii = 1: n_results
868 if ( conf>0 && spect_type(ii)==1 )
876 if ( plot_type == 1 )
877 plot(freq,[abs(spectra(:,ii)) Vxxxx]);
878 elseif ( plot_type == 2 )
879 semilogx(freq,[abs(spectra(:,ii)) Vxxxx]);
880 elseif ( plot_type == 3 )
881 semilogy(freq,[abs(spectra(:,ii)) Vxxxx]);
882 elseif ( plot_type == 4 )
883 loglog(freq,[abs(spectra(:,ii)) Vxxxx]);
884 elseif ( plot_type == 5 ) % db
885 ylabel( 'amplitude (dB)' );
886 plot(freq,[10*log10(abs(spectra(:,ii))) 10*log10(abs(Vxxxx))]);
888 title( char(plot_title(spect_type(ii),:)) );
889 ylabel( 'amplitude' );
890 %% Plot phase of cross spectrum and transfer function
891 if ( spect_type(ii)==2 || spect_type(ii)==3 )
893 if ( plot_type==2 || plot_type==4 )
894 semilogx(freq,180/pi*angle(spectra(:,ii)));
896 plot(freq,180/pi*angle(spectra(:,ii)));
898 title( char(plot_title(spect_type(ii),:)) );
908 %! rand('seed',2038014164);
909 %! a = [ 1.0 -1.6216505 1.1102795 -0.4621741 0.2075552 -0.018756746 ];
910 %! white = rand(1,16384);
911 %! signal = detrend(filter(0.70181,a,white));
912 %! % frequency shift by modulating with exp(j.omega.t)
913 %! skewed = signal.*exp(2*pi*i*2/25*[1:16384]);
914 %! Fs = 25; % sampling frequency
918 %! disp('Default settings: Fs=1Hz, overlap=0.5, no padding' )
919 %! input('Onesided power spectral density (real data). Press ENTER', 's' );
922 %! disp('Frequency-shifted complex data. Twosided wrap-around spectrum.' );
923 %! input('Area is same as one-sided spectrum. Press ENTER', 's' );
924 %! pwelch(signal,'shift','semilogy');
925 %! input('Twosided, centred zero-frequency, lin-log plot. Press ENTER', 's' );
928 %! pwelch(skewed,[],[],[],Fs,'shift','semilogy');
929 %! input('Actual Fs=25 Hz. Note change of scales. Press ENTER', 's' );
930 %! pwelch(skewed,[],[],[],Fs,0.95,'shift','semilogy');
931 %! input('Spectral density with 95% confidence interval. Press ENTER', 's' );
933 %! pwelch(signal,'squared');
934 %! input('Spectral density with Matlab R12 defaults. Press ENTER', 's' );
937 %! pwelch(signal,3640,[],4096,2*pi,[],'no-strip');
938 %! input('Same spectrum with 95% confidence interval. Press ENTER', 's' );
940 %! pwelch(signal,[],[],[],2*pi,0.95,'no-strip');
941 %! input('95% confidence interval with native defaults. Press ENTER', 's' );
942 %! pwelch(signal,64,[],[],2*pi,'no-strip');
943 %! input('Only 32 frequency values in this spectrum. Press ENTER', 's' );
945 %! pwelch(signal,64,[],256,2*pi,'no-strip');
946 %! input('4:1 zero padding gives artificial smoothing. Press ENTER', 's' );
949 %! pwelch(signal,'squared');
950 %! input('Just like Matlab spectrum.welch(...) defaults. Press ENTER', 's' );
953 %! pwelch(white,signal,'trans','coher','short')
954 %! input('Transfer and coherence functions. Press ENTER', 's' );
955 %! disp('Use "close all" to remove plotting windows.' );