## Copyright (C) 2006 Peter V. Lanspeary ## ## This program is free software: you can redistribute it and/or modify ## it under the terms of the GNU General Public License as published by ## the Free Software Foundation, either version 3 of the License, or ## (at your option) any later version. ## ## This program is distributed in the hope that it will be useful, ## but WITHOUT ANY WARRANTY; without even the implied warranty of ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the ## GNU General Public License for more details. ## ## You should have received a copy of the GNU General Public License ## along with this program; see the file COPYING. If not, see ## . ## -*- texinfo -*- ## @deftypefn {Function File} {[@var{spectra},@var{freq}] =} pwelch(@var{x}, @var{window}, @var{overlap}, @var{Nfft}, @var{Fs}, @var{range}, @var{plot_type}, @var{detrend}, @var{sloppy}) ## Estimate power spectral density of data "x" by the Welch (1967) periodogram/FFT method. ## ## All arguments except "x" are optional. ## ## The data is divided into segments. If "window" is a vector, each ## segment has the same length as "window" and is multiplied by "window" ## before (optional) zero-padding and calculation of its periodogram. If ## "window" is a scalar, each segment has a length of "window" and a ## Hamming window is used. ## ## The spectral density is the mean of the periodograms, scaled so that ## area under the spectrum is the same as the mean square of the ## data. This equivalence is supposed to be exact, but in practice there ## is a mismatch of up to 0.5% when comparing area under a periodogram ## with the mean square of the data. ## ## [spectra,freq] = pwelch(x,y,window,overlap,Nfft,Fs, ## range,plot_type,detrend,sloppy,results) ## ## Two-channel spectrum analyser. Estimate power spectral density, cross- ## spectral density, transfer function and/or coherence functions of time- ## series input data "x" and output data "y" by the Welch (1967) ## periodogram/FFT method. ## ## pwelch treats the second argument as "y" if there is a control-string ## argument "cross", "trans", "coher" or "ypower"; "power" does not force ## the 2nd argument to be treated as "y". All other arguments are ## optional. All spectra are returned in matrix "spectra". ## ## [spectra,Pxx_ci,freq] = pwelch(x,window,overlap,Nfft,Fs,conf, ## range,plot_type,detrend,sloppy) ## ## [spectra,Pxx_ci,freq] = pwelch(x,y,window,overlap,Nfft,Fs,conf, ## range,plot_type,detrend,sloppy,results) ## ## Estimates confidence intervals for the spectral density. ## ## See Hint (7) below for compatibility options. ## ## Confidence level "conf" is the 6th or 7th numeric argument. If "results" control-string ## arguments are used, one of them must be "power" when the "conf" ## argument is present; pwelch can estimate confidence intervals only for ## the power spectrum of the "x" data. It does not know how to estimate ## confidence intervals of the cross-power spectrum, transfer function or ## coherence; if you can suggest a good method, please send a bug report. ## ## ARGUMENTS ## ## All but the first argument are optional and may be empty, except that ## the "results" argument may require the second argument to be "y". ## @table @asis ## @item x ## [non-empty vector] system-input time-series data ## @item y ## [non-empty vector] system-output time-series data ## @item window ## [real vector] of window-function values; the ## data segment has the same length as the window. ## Default window shape is Hamming. ## ## [integer scalar] length of each data segment. The default ## value is window=sqrt(length(x)) rounded up to the ## nearest integer power of 2; see 'sloppy' argument. ## ## @item overlap ## [real scalar] segment overlap expressed as a multiple of ## window or segment length. 0 <= overlap < 1, ## The default is overlap=0.5 . ## ## @item Nfft ## [integer scalar] Length of FFT. The default is the length ## of the "window" vector or has the same value as the ## scalar "window" argument. If Nfft is larger than the ## segment length, "seg_len", the data segment is padded ## with "Nfft-seg_len" zeros. The default is no padding. ## Nfft values smaller than the length of the data ## segment (or window) are ignored silently. ## ## @item Fs ## [real scalar] sampling frequency (Hertz); default=1.0 ## ## @item conf ## [real scalar] confidence level between 0 and 1. Confidence ## intervals of the spectral density are estimated from ## scatter in the periodograms and are returned as Pxx_ci. ## Pxx_ci(:,1) is the lower bound of the confidence ## interval and Pxx_ci(:,2) is the upper bound. If there ## are three return values, or conf is an empty matrix, ## confidence intervals are calculated for conf=0.95 . ## If conf is zero or is not given, confidence intervals ## are not calculated. Confidence intervals can be ## obtained only for the power spectral density of x; ## nothing else. ## @end table ## ## CONTROL-STRING ARGUMENTS -- each of these arguments is a character string. ## Control-string arguments must be after the other arguments but can be in ## any order. ## ## @table @asis ## @item range ## 'half', 'onesided' : frequency range of the spectrum is ## zero up to but not including Fs/2. Power from ## negative frequencies is added to the positive side of ## the spectrum, but not at zero or Nyquist (Fs/2) ## frequencies. This keeps power equal in time and ## spectral domains. See reference [2]. ## ## 'whole', 'twosided' : frequency range of the spectrum is ## -Fs/2 to Fs/2, with negative frequencies ## stored in "wrap around" order after the positive ## frequencies; e.g. frequencies for a 10-point 'twosided' ## spectrum are 0 0.1 0.2 0.3 0.4 0.5 -0.4 -0.3 -0.2 -0.1 ## ## 'shift', 'centered' : same as 'whole' but with the first half ## of the spectrum swapped with second half to put the ## zero-frequency value in the middle. See "help ## fftshift". ## ## If data (x and y) are real, the default range is 'half', ## otherwise default range is 'whole'. ## ## @item plot_type ## 'plot', 'semilogx', 'semilogy', 'loglog', 'squared' or 'db': ## specifies the type of plot. The default is 'plot', which ## means linear-linear axes. 'squared' is the same as 'plot'. ## 'dB' plots "10*log10(psd)". This argument is ignored and a ## spectrum is not plotted if the caller requires a returned ## value. ## ## @item detrend ## 'no-strip', 'none' -- do NOT remove mean value from the data ## ## 'short', 'mean' -- remove the mean value of each segment from ## each segment of the data. ## ## 'linear', -- remove linear trend from each segment of ## the data. ## ## 'long-mean' -- remove the mean value from the data before ## splitting it into segments. This is the default. ## ## @item sloppy ## 'sloppy': FFT length is rounded up to the nearest integer ## power of 2 by zero padding. FFT length is adjusted ## after addition of padding by explicit Nfft argument. ## The default is to use exactly the FFT and window/ ## segment lengths specified in argument list. ## ## @item results ## specifies what results to return (in the order specified ## and as many as desired). ## ## 'power' calculate power spectral density of "x" ## ## 'cross' calculate cross spectral density of "x" and "y" ## ## 'trans' calculate transfer function of a system with ## input "x" and output "y" ## ## 'coher' calculate coherence function of "x" and "y" ## ## 'ypower' calculate power spectral density of "y" ## ## The default is 'power', with argument "y" omitted. ## @end table ## ## RETURNED VALUES: ## ## If return values are not required by the caller, the results are ## plotted and nothing is returned. ## ## @table @asis ## @item spectra ## [real-or-complex matrix] columns of the matrix contain results ## in the same order as specified by "results" arguments. ## Each column contains one of the result vectors. ## ## @item Pxx_ci ## [real matrix] estimate of confidence interval for power ## spectral density of x. First column is the lower ## bound. Second column is the upper bound. ## ## @item freq ## [real column vector] frequency values ## @end table ## ## HINTS ## @enumerate ## @item EMPTY ARGS: ## if you don't want to use an optional argument you can leave it empty ## by writing its value as []. ## ## @item FOR BEGINNERS: ## The profusion of arguments may make pwelch difficult to use, and an ## unskilled user can easily produce a meaningless result or can easily ## mis-interpret the result. With real data "x" and sampling frequency ## "Fs", the easiest and best way for a beginner to use pwelch is ## probably "pwelch(x,[],[],[],Fs)". Use the "window" argument to ## control the length of the spectrum vector. For real data and integer ## scalar M, "pwelch(x,2*M,[],[],Fs)" gives an M+1 point spectrum. ## Run "demo pwelch" (octave only). ## ## @item WINDOWING FUNCTIONS: ## Without a window function, sharp spectral peaks can have strong ## sidelobes because the FFT of a data in a segment is in effect convolved ## with a rectangular window. A window function which tapers off ## (gradually) at the ends produces much weaker sidelobes in the FFT. ## Hann (hanning), hamming, bartlett, blackman, flattopwin etc are ## available as separate Matlab/sigproc or Octave functions. The sidelobes ## of the Hann window have a roll-off rate of 60dB/decade of frequency. ## The first sidelobe of the Hamming window is suppressed and is about 12dB ## lower than the first Hann sidelobe, but the roll-off rate is only ## 20dB/decade. You can inspect the FFT of a Hann window by plotting ## "abs(fft(postpad(hanning(256),4096,0)))". ## The default window is Hamming. ## ## @item ZERO PADDING: ## Zero-padding reduces the frequency step in the ## spectrum, and produces an artificially smoothed spectrum. For example, ## "Nfft=2*length(window)" gives twice as many frequency values, but ## adjacent PSD (power spectral density) values are not independent; ## adjacent PSD values are independent if "Nfft=length(window)", which is ## the default value of Nfft. ## ## @item REMOVING MEAN FROM SIGNAL: ## If the mean is not removed from the signal there is a large spectral ## peak at zero frequency and the sidelobes of this peak are likely to ## swamp the rest of the spectrum. For this reason, the default behavior ## is to remove the mean. However, the matlab pwelch does not do this. ## ## @item WARNING ON CONFIDENCE INTERVALS ## Confidence intervals are obtained by measuring the sample variance of ## the periodograms and assuming that the periodograms have a Gaussian ## probability distribution. This assumption is not accurate. If, for ## example, the data (x) is Gaussian, the periodogram has a Rayleigh ## distribution. However, the confidence intervals may still be useful. ## ## @item COMPATIBILITY WITH Matlab R11, R12, etc ## When used without the second data (y) argument, arguments are compatible ## with the pwelch of Matlab R12, R13, R14, 2006a and 2006b except that ## ## 1) overlap is expressed as a multiple of window length --- ## effect of overlap scales with window length ## ## 2) default values of length(window), Nfft and Fs are more sensible, and ## ## 3) Goertzel algorithm is not available so Nfft cannot be an array of ## frequencies as in Matlab 2006b. ## ## Pwelch has four persistent Matlab-compatibility levels. Calling pwelch ## with an empty first argument sets the order of arguments and defaults ## specified above in the USAGE and ARGUMENTS section of this documentation. ## @example ## prev_compat=pwelch([]); ## [Pxx,f]=pwelch(x,window,overlap,Nfft,Fs,conf,...); ## @end example ## Calling pwelch with a single string argument (as described below) gives ## compatibility with Matlab R11 or R12, or the R14 spectrum.welch ## defaults. The returned value is the PREVIOUS compatibility string. ## ## Matlab R11: For compatibility with the Matlab R11 pwelch: ## @example ## prev_compat=pwelch('R11-'); ## [Pxx,f]=pwelch(x,Nfft,Fs,window,overlap,conf,range,units); ## %% units of overlap are "number of samples" ## %% defaults: Nfft=min(length(x),256), Fs=2*pi, length(window)=Nfft, ## %% window=Hanning, do not detrend, ## %% N.B. "Sloppy" is not available. ## @end example ## ## Matlab R12: For compatibility with Matlab R12 to 2006a pwelch: ## @example ## prev_compat=pwelch('R12+'); ## [Pxx,f]=pwelch(x,window,overlap,nfft,Fs,...); ## %% units of overlap are "number of samples" ## %% defaults: length(window)==length(x)/8, window=Hamming, ## %% Nfft=max(256,NextPow2), Fs=2*pi, do not detrend ## %% NextPow2 is the next power of 2 greater than or equal to the ## %% window length. "Sloppy", "conf" are not available. Default ## %% window length gives very poor amplitude resolution. ## @end example ## ## To adopt defaults of the Matlab R14 "spectrum.welch" spectrum object ## associated "psd" method. ## @example ## prev_compat=pwelch('psd'); ## [Pxx,f] = pwelch(x,window,overlap,Nfft,Fs,conf,...); ## %% overlap is expressed as a percentage of window length, ## %% defaults: length(window)==64, Nfft=max(256,NextPow2), Fs=2*pi ## %% do not detrend ## %% NextPow2 is the next power of 2 greater than or equal to the ## %% window length. "Sloppy" is not available. ## %% Default window length gives coarse frequency resolution. ## @end example ## @end enumerate ## ## REFERENCES ## ## [1] Peter D. Welch (June 1967): ## "The use of fast Fourier transform for the estimation of power spectra: ## a method based on time averaging over short, modified periodograms." ## IEEE Transactions on Audio Electroacoustics, Vol AU-15(6), pp 70-73 ## ## [2] William H. Press and Saul A. Teukolsky and William T. Vetterling and ## Brian P. Flannery, ## "Numerical recipes in C, The art of scientific computing", 2nd edition, ## Cambridge University Press, 2002 --- Section 13.7. ## @end deftypefn function varargout = pwelch(x,varargin) ## ## COMPATIBILITY LEVEL ## Argument positions and defaults depend on compatibility level selected ## by calling pwelch without arguments or with a single string argument. ## native: compatib=1; prev_compat=pwelch(); prev_compat=pwelch([]); ## matlab R11: compatib=2; prev_compat=pwelch('R11-'); ## matlab R12: compatib=3; prev_compat=pwelch('R12+'); ## spectrum.welch defaults: compatib=4; prev_compat=pwelch('psd'); ## In each case, the returned value is the PREVIOUS compatibility string. ## compat_str = {[]; 'R11-'; 'R12+'; 'psd'}; persistent compatib; if ( isempty(compatib) || compatib<=0 || compatib>4 ) ## legal values are 1, 2, 3, 4 compatib = 1; endif if ( nargin <= 0 ) error( 'pwelch: Need at least 1 arg. Use "help pwelch".' ); elseif ( nargin==1 && (ischar(x) || isempty(x)) ) varargout{1} = compat_str{compatib}; if ( isempty(x) ) # native compatib = 1; elseif ( strcmp(x,'R11-') ) compatib = 2; elseif ( strcmp(x,'R12+') ) compatib = 3; elseif ( strcmp(x,'psd') ) compatib = 4; else error( 'pwelch: compatibility arg must be empty, R11-, R12+ or psd' ); endif ## return ## ## Check fixed argument elseif ( isempty(x) || ~isvector(x) ) error( 'pwelch: arg 1 (x) must be vector.' ); else ## force x to be COLUMN vector if ( size(x,1)==1 ) x=x(:); endif ## ## Look through all args to check if cross PSD, transfer function or ## coherence is required. If yes, the second arg is data vector "y". arg2_is_y = 0; x_len = length(x); nvarargin = length(varargin); for iarg=1:nvarargin arg = varargin{iarg}; if ( ~isempty(arg) && ischar(arg) && ... ( strcmp(arg,'cross') || strcmp(arg,'trans') || ... strcmp(arg,'coher') || strcmp(arg,'ypower') )) ## OK. Need "y". Grab it from 2nd arg. arg = varargin{1}; if ( nargin<2 || isempty(arg) || ~isvector(arg) || length(arg)~=x_len ) error( 'pwelch: arg 2 (y) must be vector, same length as x.' ); endif ## force COLUMN vector y = varargin{1}(:); arg2_is_y = 1; break; endif endfor ## ## COMPATIBILITY ## To select default argument values, "compatib" is used as an array index. ## Index values are 1=native, 2=R11, 3=R12, 4=spectrum.welch ## ## argument positions: ## arg_posn = varargin index of window, overlap, Nfft, Fs and conf ## args respectively, a value of zero ==>> arg does not exist arg_posn = [1 2 3 4 5; # native 3 4 1 2 5; # Matlab R11- pwelch 1 2 3 4 0; # Matlab R12+ pwelch 1 2 3 4 5]; # spectrum.welch defaults arg_posn = arg_posn(compatib,:) + arg2_is_y; ## ## SPECIFY SOME DEFAULT VALUES for (not all) optional arguments ## Use compatib as array index. ## Fs = sampling frequency Fs = [ 1.0 2*pi 2*pi 2*pi ]; Fs = Fs(compatib); ## plot_type: 1='plot'|'squared'; 5='db'|'dB' plot_type = [ 1 5 5 5 ]; plot_type = plot_type(compatib); ## rm_mean: 3='long-mean'; 0='no-strip'|'none' rm_mean = [ 3 0 0 0 ]; rm_mean = rm_mean(compatib); ## use max_overlap=x_len-1 because seg_len is not available yet ## units of overlap are different for each version: ## fraction, samples, or percent max_overlap = [ 0.95 x_len-1 x_len-1 95]; max_overlap = max_overlap(compatib); ## default confidence interval ## if there are more than 2 return values and if there is a "conf" arg conf = 0.95 * (nargout>2) * (arg_posn(5)>0); ## is_win = 0; # =0 means valid window arg is not provided yet Nfft = []; # default depends on segment length overlap = []; # WARNING: units can be #samples, fraction or percentage range = ~isreal(x) || ( arg2_is_y && ~isreal(y) ); is_sloppy = 0; n_results = 0; do_power = 0; do_cross = 0; do_trans = 0; do_coher = 0; do_ypower = 0; ## ## DECODE AND CHECK OPTIONAL ARGUMENTS end_numeric_args = 0; for iarg = 1+arg2_is_y:nvarargin arg = varargin{iarg}; if ( ischar(arg) ) ## first string arg ==> no more numeric args ## non-string args cannot follow a string arg end_numeric_args = 1; ## ## decode control-string arguments if ( strcmp(arg,'sloppy') ) is_sloppy = ~is_win || is_win==1; elseif ( strcmp(arg,'plot') || strcmp(arg,'squared') ) plot_type = 1; elseif ( strcmp(arg,'semilogx') ) plot_type = 2; elseif ( strcmp(arg,'semilogy') ) plot_type = 3; elseif ( strcmp(arg,'loglog') ) plot_type = 4; elseif ( strcmp(arg,'db') || strcmp(arg,'dB') ) plot_type = 5; elseif ( strcmp(arg,'half') || strcmp(arg,'onesided') ) range = 0; elseif ( strcmp(arg,'whole') || strcmp(arg,'twosided') ) range = 1; elseif ( strcmp(arg,'shift') || strcmp(arg,'centered') ) range = 2; elseif ( strcmp(arg,'long-mean') ) rm_mean = 3; elseif ( strcmp(arg,'linear') ) rm_mean = 2; elseif ( strcmp(arg,'short') || strcmp(arg,'mean') ) rm_mean = 1; elseif ( strcmp(arg,'no-strip') || strcmp(arg,'none') ) rm_mean = 0; elseif ( strcmp(arg, 'power' ) ) if ( ~do_power ) n_results = n_results+1; do_power = n_results; endif elseif ( strcmp(arg, 'cross' ) ) if ( ~do_cross ) n_results = n_results+1; do_cross = n_results; endif elseif ( strcmp(arg, 'trans' ) ) if ( ~do_trans ) n_results = n_results+1; do_trans = n_results; endif elseif ( strcmp(arg, 'coher' ) ) if ( ~do_coher ) n_results = n_results+1; do_coher = n_results; endif elseif ( strcmp(arg, 'ypower' ) ) if ( ~do_ypower ) n_results = n_results+1; do_ypower = n_results; endif else error( 'pwelch: string arg %d illegal value: %s', iarg+1, arg ); endif ## end of processing string args ## elseif ( end_numeric_args ) if ( ~isempty(arg) ) ## found non-string arg after a string arg ... oops error( 'pwelch: control arg must be string' ); endif ## ## first 4 optional arguments are numeric -- in fixed order ## ## deal with "Fs" and "conf" first because empty arg is a special default ## -- "Fs" arg -- sampling frequency elseif ( iarg == arg_posn(4) ) if ( isempty(arg) ) Fs = 1; elseif ( ~isscalar(arg) || ~isreal(arg) || arg<0 ) error( 'pwelch: arg %d (Fs) must be real scalar >0', iarg+1 ); else Fs = arg; endif ## ## -- "conf" arg -- confidence level ## guard against the "it cannot happen" iarg==0 elseif ( arg_posn(5) && iarg == arg_posn(5) ) if ( isempty(arg) ) conf = 0.95; elseif ( ~isscalar(arg) || ~isreal(arg) || arg < 0.0 || arg >= 1.0 ) error( 'pwelch: arg %d (conf) must be real scalar, >=0, <1',iarg+1 ); else conf = arg; endif ## ## skip all empty args from this point onward elseif ( isempty(arg) ) 1; ## ## -- "window" arg -- window function elseif ( iarg == arg_posn(1) ) if ( isscalar(arg) ) is_win = 1; elseif ( isvector(arg) ) is_win = length(arg); if ( size(arg,2)>1 ) # vector must be COLUMN vector arg = arg(:); endif else is_win = 0; endif if ( ~is_win ) error( 'pwelch: arg %d (window) must be scalar or vector', iarg+1 ); elseif ( is_win==1 && ( ~isreal(arg) || fix(arg)~=arg || arg<=3 ) ) error( 'pwelch: arg %d (window) must be integer >3', iarg+1 ); elseif ( is_win>1 && ( ~isreal(arg) ) ) error( 'pwelch: arg %d (window) vector must be real',iarg+1); endif window = arg; is_sloppy = 0; ## ## -- "overlap" arg -- segment overlap elseif ( iarg == arg_posn(2) ) if (~isscalar(arg) || ~isreal(arg) || arg<0 || arg>max_overlap ) error( 'pwelch: arg %d (overlap) must be real from 0 to %f', ... iarg+1, max_overlap ); endif overlap = arg; ## ## -- "Nfft" arg -- FFT length elseif ( iarg == arg_posn(3) ) if ( ~isscalar(arg) || ~isreal(arg) || fix(arg)~=arg || arg<0 ) error( 'pwelch: arg %d (Nfft) must be integer >=0', iarg+1 ); endif Nfft = arg; ## else error( 'pwelch: arg %d must be string', iarg+1 ); endif endfor if ( conf>0 && (n_results && ~do_power ) ) error('pwelch: can give confidence interval for x power spectrum only' ); endif ## ## end DECODE AND CHECK OPTIONAL ARGUMENTS. ## ## SETUP REMAINING PARAMETERS ## default action is to calculate power spectrum only if ( ~n_results ) n_results = 1; do_power = 1; endif need_Pxx = do_power || do_trans || do_coher; need_Pxy = do_cross || do_trans || do_coher; need_Pyy = do_coher || do_ypower; log_two = log(2); nearly_one = 0.99999999999; ## ## compatibility-options ## provides exact compatibility with Matlab R11 or R12 ## ## Matlab R11 compatibility if ( compatib==2 ) if ( isempty(Nfft) ) Nfft = min( 256, x_len ); endif if ( is_win > 1 ) seg_len = min( length(window), Nfft ); window = window(1:seg_len); else if ( is_win ) ## window arg is scalar seg_len = window; else seg_len = Nfft; endif ## make Hann window (don't depend on sigproc) xx = seg_len - 1; window = 0.5 - 0.5 * cos( (2*pi/xx)*[0:xx].' ); endif ## ## Matlab R12 compatibility elseif ( compatib==3 ) if ( is_win > 1 ) ## window arg provides window function seg_len = length(window); else ## window arg does not provide window function; use Hamming if ( is_win ) ## window arg is scalar seg_len = window; else ## window arg not available; use R12 default, 8 windows ## ignore overlap arg; use overlap=50% -- only choice that makes sense ## this is the magic formula for 8 segments with 50% overlap seg_len = fix( (x_len-3)*2/9 ); endif ## make Hamming window (don't depend on sigproc) xx = seg_len - 1; window = 0.54 - 0.46 * cos( (2*pi/xx)*[0:xx].' ); endif if ( isempty(Nfft) ) Nfft = max( 256, 2^ceil(log(seg_len)*nearly_one/log_two) ); endif ## ## Matlab R14 psd(spectrum.welch) defaults elseif ( compatib==4 ) if ( is_win > 1 ) ## window arg provides window function seg_len = length(window); else ## window arg does not provide window function; use Hamming if ( is_win ) ## window arg is scalar seg_len = window; else ## window arg not available; use default seg_len = 64 seg_len = 64; endif ## make Hamming window (don't depend on sigproc) xx = seg_len - 1; window = 0.54 - 0.46 * cos( (2*pi/xx)*[0:xx].' ); endif ## Now we know segment length, ## so we can set default overlap as number of samples if ( ~isempty(overlap) ) overlap = fix(seg_len * overlap / 100 ); endif if ( isempty(Nfft) ) Nfft = max( 256, 2^ceil(log(seg_len)*nearly_one/log_two) ); endif ## ## default compatibility level else # if ( compatib==1 ) ## calculate/adjust segment length, window function if ( is_win > 1 ) ## window arg provides window function seg_len = length(window); else ## window arg does not provide window function; use Hamming if ( is_win ) # window arg is scalar seg_len = window; else ## window arg not available; use default length: ## = sqrt(length(x)) rounded up to nearest integer power of 2 if ( isempty(overlap) ) overlap=0.5; endif seg_len = 2 ^ ceil( log(sqrt(x_len/(1-overlap)))*nearly_one/log_two ); endif ## make Hamming window (don't depend on sigproc) xx = seg_len - 1; window = 0.54 - 0.46 * cos( (2*pi/xx)*[0:xx].' ); endif ## Now we know segment length, ## so we can set default overlap as number of samples if ( ~isempty(overlap) ) overlap = fix(seg_len * overlap); endif ## ## calculate FFT length if ( isempty(Nfft) ) Nfft = seg_len; endif if ( is_sloppy ) Nfft = 2 ^ ceil( log(Nfft) * nearly_one / log_two ); endif endif ## end of compatibility options ## ## minimum FFT length is seg_len Nfft = max( Nfft, seg_len ); ## Mean square of window is required for normalizing PSD amplitude. win_meansq = (window.' * window) / seg_len; ## ## Set default or check overlap. if ( isempty(overlap) ) overlap = fix(seg_len /2); elseif ( overlap >= seg_len ) error( 'pwelch: arg (overlap=%d) too big. Must be 0 ) Vxx = xx; else Vxx = []; endif n_ffts = 0; for start_seg = [1:seg_len-overlap:x_len-seg_len+1] end_seg = start_seg+seg_len-1; ## Don't truncate/remove the zero padding in xx and yy if ( need_Pxx || need_Pxy ) if ( rm_mean==1 ) # remove mean from segment xx(1:seg_len) = window .* ( ... x(start_seg:end_seg) - sum(x(start_seg:end_seg)) / seg_len); elseif ( rm_mean == 2 ) # remove linear trend from segment xx(1:seg_len) = window .* detrend( x(start_seg:end_seg) ); else # rm_mean==0 or 3 xx(1:seg_len) = window .* x(start_seg:end_seg); endif fft_x = fft(xx); endif if ( need_Pxy || need_Pyy ) if ( rm_mean==1 ) # remove mean from segment yy(1:seg_len) = window .* ( ... y(start_seg:end_seg) - sum(y(start_seg:end_seg)) / seg_len); elseif ( rm_mean == 2 ) # remove linear trend from segment yy(1:seg_len) = window .* detrend( y(start_seg:end_seg) ); else # rm_mean==0 or 3 yy(1:seg_len) = window .* y(start_seg:end_seg); endif fft_y = fft(yy); endif if ( need_Pxx ) ## force Pxx to be real; pgram = periodogram pgram = real(fft_x .* conj(fft_x)); Pxx = Pxx + pgram; ## sum of squared periodograms is required for confidence interval if ( conf>0 ) Vxx = Vxx + pgram .^2; endif endif if ( need_Pxy ) ## Pxy (cross power spectrum) is complex. Do not force to be real. Pxy = Pxy + fft_y .* conj(fft_x); endif if ( need_Pyy ) ## force Pyy to be real Pyy = Pyy + real(fft_y .* conj(fft_y)); endif n_ffts = n_ffts +1; endfor ## ## Calculate confidence interval ## -- incorrectly assumes that the periodogram has Gaussian probability ## distribution (actually, it has a single-sided (e.g. exponential) ## distribution. ## Sample variance of periodograms is (Vxx-Pxx.^2/n_ffts)/(n_ffts-1). ## This method of calculating variance is more susceptible to round-off ## error, but is quicker, and for double-precision arithmetic and the ## inherently noisy periodogram (variance==mean^2), it should be OK. if ( conf>0 && need_Pxx ) if ( n_ffts<2 ) Vxx = zeros(Nfft,1); else ## Should use student distribution here (for unknown variance), but tinv ## is not a core Matlab function (is in statistics toolbox. Grrr) Vxx = (erfinv(conf)*sqrt(2*n_ffts/(n_ffts-1))) * sqrt(Vxx-Pxx.^2/n_ffts); endif endif ## ## Convert two-sided spectra to one-sided spectra (if range == 0). ## For one-sided spectra, contributions from negative frequencies are added ## to the positive side of the spectrum -- but not at zero or Nyquist ## (half sampling) frequencies. This keeps power equal in time and spectral ## domains, as required by Parseval theorem. ## if ( range == 0 ) if ( ~ rem(Nfft,2) ) # one-sided, Nfft is even psd_len = Nfft/2+1; if ( need_Pxx ) Pxx = Pxx(1:psd_len) + [0; Pxx(Nfft:-1:psd_len+1); 0]; if ( conf>0 ) Vxx = Vxx(1:psd_len) + [0; Vxx(Nfft:-1:psd_len+1); 0]; endif endif if ( need_Pxy ) Pxy = Pxy(1:psd_len) + conj([0; Pxy(Nfft:-1:psd_len+1); 0]); endif if ( need_Pyy ) Pyy = Pyy(1:psd_len) + [0; Pyy(Nfft:-1:psd_len+1); 0]; endif else # one-sided, Nfft is odd psd_len = (Nfft+1)/2; if ( need_Pxx ) Pxx = Pxx(1:psd_len) + [0; Pxx(Nfft:-1:psd_len+1)]; if ( conf>0 ) Vxx = Vxx(1:psd_len) + [0; Vxx(Nfft:-1:psd_len+1)]; endif endif if ( need_Pxy ) Pxy = Pxy(1:psd_len) + conj([0; Pxy(Nfft:-1:psd_len+1)]); endif if ( need_Pyy ) Pyy = Pyy(1:psd_len) + [0; Pyy(Nfft:-1:psd_len+1)]; endif endif else # two-sided (and shifted) psd_len = Nfft; endif ## end MAIN CALCULATIONS ## ## SCALING AND OUTPUT ## Put all results in matrix, one row per spectrum ## Pxx, Pxy, Pyy are sums of periodograms, so "n_ffts" ## in the scale factor converts them into averages spectra = zeros(psd_len,n_results); spect_type = zeros(n_results,1); scale = n_ffts * seg_len * Fs * win_meansq; if ( do_power ) spectra(:,do_power) = Pxx / scale; spect_type(do_power) = 1; if ( conf>0 ) Vxx = [Pxx-Vxx Pxx+Vxx]/scale; endif endif if ( do_cross ) spectra(:,do_cross) = Pxy / scale; spect_type(do_cross) = 2; endif if ( do_trans ) spectra(:,do_trans) = Pxy ./ Pxx; spect_type(do_trans) = 3; endif if ( do_coher ) ## force coherence to be real spectra(:,do_coher) = real(Pxy .* conj(Pxy)) ./ Pxx ./ Pyy; spect_type(do_coher) = 4; endif if ( do_ypower ) spectra(:,do_ypower) = Pyy / scale; spect_type(do_ypower) = 5; endif freq = [0:psd_len-1].' * ( Fs / Nfft ); ## ## range='shift': Shift zero-frequency to the middle if ( range == 2 ) len2 = fix((Nfft+1)/2); spectra = [ spectra(len2+1:Nfft,:); spectra(1:len2,:)]; freq = [ freq(len2+1:Nfft)-Fs; freq(1:len2)]; if ( conf>0 ) Vxx = [ Vxx(len2+1:Nfft,:); Vxx(1:len2,:)]; endif endif ## ## RETURN RESULTS or PLOT if ( nargout>=2 && conf>0 ) varargout{2} = Vxx; endif if ( nargout>=(2+(conf>0)) ) ## frequency is 2nd or 3rd return value, ## depends on if 2nd is confidence interval varargout{2+(conf>0)} = freq; endif if ( nargout>=1 ) varargout{1} = spectra; else ## ## Plot the spectra if there are no return variables. plot_title=['power spectrum x '; 'cross spectrum '; 'transfer function'; 'coherence '; 'power spectrum y ' ]; for ii = 1: n_results if ( conf>0 && spect_type(ii)==1 ) Vxxxx = Vxx; else Vxxxx = []; endif if ( n_results > 1 ) figure(); endif if ( plot_type == 1 ) plot(freq,[abs(spectra(:,ii)) Vxxxx]); elseif ( plot_type == 2 ) semilogx(freq,[abs(spectra(:,ii)) Vxxxx]); elseif ( plot_type == 3 ) semilogy(freq,[abs(spectra(:,ii)) Vxxxx]); elseif ( plot_type == 4 ) loglog(freq,[abs(spectra(:,ii)) Vxxxx]); elseif ( plot_type == 5 ) # db %ylabel( 'amplitude (dB)' ); plot(freq,[10*log10(abs(spectra(:,ii))) 10*log10(abs(Vxxxx))]); endif title( char(plot_title(spect_type(ii),:)) ); grid on; ylabel( 'power' ); xlabel( 'frequency' ); ## Plot phase of cross spectrum and transfer function if ( spect_type(ii)==2 || spect_type(ii)==3 ) figure(); if ( plot_type==2 || plot_type==4 ) semilogx(freq,180/pi*angle(spectra(:,ii))); else plot(freq,180/pi*angle(spectra(:,ii))); endif title( char(plot_title(spect_type(ii),:)) ); ylabel( 'phase' ); endif endfor endif endif endfunction %!demo %! a = [ 1.0 -1.6216505 1.1102795 -0.4621741 0.2075552 -0.018756746 ]; %! white = rand(1,16384); %! signal = detrend(filter(0.70181,a,white)); %! % frequency shift by modulating with exp(j.omega.t) %! skewed = signal.*exp(2*pi*i*2/25*[1:16384]); %! compat = pwelch ([]); %! hold on; %! pwelch(signal); %! pwelch(skewed); %! pwelch(signal,'shift','semilogy'); %! pwelch (compat); %! hold off; %!demo %! Fs = 25; %! a = [ 1.0 -1.6216505 1.1102795 -0.4621741 0.2075552 -0.018756746 ]; %! white = rand(1,16384); %! signal = detrend(filter(0.70181,a,white)); %! % frequency shift by modulating with exp(j.omega.t) %! skewed = signal.*exp(2*pi*i*2/25*[1:16384]); %! compat = pwelch ([]); %! pwelch(skewed,[],[],[],Fs,'shift','semilogy'); %! pwelch(skewed,[],[],[],Fs,0.95,'shift','semilogy'); %! pwelch('R12+'); %! pwelch(signal,'squared'); %! pwelch (compat); %!demo %! a = [ 1.0 -1.6216505 1.1102795 -0.4621741 0.2075552 -0.018756746 ]; %! white = rand(1,16384); %! signal = detrend(filter(0.70181,a,white)); %! compat = pwelch ([]); %! pwelch(signal,3640,[],4096,2*pi,[],'no-strip'); %! pwelch (compat); %!demo %! a = [ 1.0 -1.6216505 1.1102795 -0.4621741 0.2075552 -0.018756746 ]; %! white = rand(1,16384); %! signal = detrend(filter(0.70181,a,white)); %! compat = pwelch ([]); %! hold on; %! pwelch(signal,[],[],[],2*pi,0.95,'no-strip'); %! pwelch(signal,64,[],[],2*pi,'no-strip'); %! pwelch(signal,64,[],256,2*pi,'no-strip'); %! pwelch (compat); %! hold off; %!demo %! a = [ 1.0 -1.6216505 1.1102795 -0.4621741 0.2075552 -0.018756746 ]; %! white = rand(1,16384); %! signal = detrend(filter(0.70181,a,white)); %! compat = pwelch ('psd'); %! pwelch(signal,'squared'); %! pwelch({}); %! pwelch(white,signal,'trans','coher','short') %! pwelch (compat);