diff options
Diffstat (limited to 'macros/pwelch.sci')
| -rw-r--r-- | macros/pwelch.sci | 882 |
1 files changed, 817 insertions, 65 deletions
diff --git a/macros/pwelch.sci b/macros/pwelch.sci index 3e4aa8d..e1c0e17 100644 --- a/macros/pwelch.sci +++ b/macros/pwelch.sci @@ -1,67 +1,819 @@ -function [spectra, Pxx_ci, freqn] = pwelch(x, varargin) -//Estimate power spectral density of data "x" by the Welch (1967) periodogram/FFT method. -//Calling Sequence -//[spectra,freq] = pwelch(x, window, overlap, Nfft, Fs, range, plot_type, detrend, sloppy) -//[spectra,freq] = pwelch(x, y, window, overlap, Nfft, Fs, range, plot_type, detrend, sloppy, results) -//[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) -//Parameters -//x: [non-empty vector] system-input time-series data -// -//y: [non-empty vector] system-output time-series data -// -//window: [real vector] of window-function values between 0 and 1; 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. -// -//overlap: [real scalar] segment overlap expressed as a multiple of window or segment length. 0 <= overlap < 1, The default is overlap=0.5 . -// -//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. -// -//Fs:[real scalar] sampling frequency (Hertz); default=1.0 -// -//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. -// -//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', 'centerdc' : 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'. -// -//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. -// -//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. -// -//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. -//Description -//Estimate power spectral density of data "x" by the Welch (1967) periodogram/FFT method. 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. -funcprot(0); -rhs=argn(2); -lhs=argn(1) -if(rhs<9 | rhs>9) then - error("Wrong number of input arguments."); -end -if(rhs<11 | rhs>11) then - error("Wrong number of input arguments."); -end -if(lhs<2 | lhs>3) then - error("Wrong number of output arguments."); -end -select(rhs) -case 9 then - select(lhs) - case 2 then - [spectra, freqn] = callOctave("pwelch", x, varargin(1), varargin(2), varargin(3), varargin(4), varargin(5), varargin(6), varargin(7), varargin(8)); - case 3 then - [spectra, Pxx_ci, freqn] = callOctave("pwelch", x, varargin(1), varargin(2), varargin(3), varargin(4), varargin(5), varargin(6), varargin(7), varargin(8)); +/* +Dependencies : fft1 +Description: + Estimate power spectral density of data "x" by the Welch (1967) periodogram/FFT method. + All arguments except "x" are optional. +Calling Sequence: + 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". + x : [non-empty vector] system-input time-series data + y : [non-empty vector] system-output time-series data + window : [real vector] of window-function values between 0 and 1; 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. + overlap: [real scalar] segment overlap expressed as a multiple of window or segment length. 0 <= overlap < 1, The default is overlap=0.5 . + 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. + Fs : [real scalar] sampling frequency (Hertz); default=1.0 + 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. + 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. + 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’, ’centerdc’ : 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’. + 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. + 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. + 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. + 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. + RETURNED VALUES: + If return values are not required by the caller, the results are plotted and nothing is returned. + 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. + 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. + freq : [real column vector] frequency values + HINTS + EMPTY ARGS: if you don’t want to use an optional argument you can leave it empty by writing its value as []. + 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 "help pwelch". + 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. + 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. + 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. + 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. + 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. + prev_compat=pwelch([]); + [Pxx,f]=pwelch(x,window,overlap,Nfft,Fs,conf,...); + 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: + 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. + Matlab R12: For compatibility with Matlab R12 to 2006a pwelch: + 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. + To adopt defaults of the Matlab R14 "spectrum.welch" spectrum object associated "psd" method. + 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. + */ +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. + // + nargin = argn(2) + compat_str = {[]; 'R11-'; 'R12+'; 'psd'}; + global compatib + if ( isempty(compatib) || compatib<=0 || compatib>4 ) + // legal values are 1, 2, 3, 4 + compatib = 1; end -case 11 then - select(lhs) - case 2 then - [spectra, freqn] = callOctave("pwelch", x, varargin(1), varargin(2), varargin(3), varargin(4), varargin(5), varargin(6), varargin(7), varargin(8), varargin(9), varargin(10)); - case 3 then - [spectra, Pxx_ci, freqn] = callOctave("pwelch", x, varargin(1), varargin(2), varargin(3), varargin(4), varargin(5), varargin(6), varargin(7), varargin(8), varargin(9), varargin(10)); + if ( nargin <= 0 ) + error( 'pwelch: Need at least 1 arg. Use help pwelch.' ); + elseif ( nargin==1 && (type(x) == 10 ) || 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' ); + end + // 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(:); + end + // + // 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 = max(size(x)); + nvarargin = max(size(varargin)); + for iarg=1:nvarargin + arg = varargin(iarg); + if ( ~isempty(arg) && (type(arg) == 10 ) && ... + ( ~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) || max(size(arg))~=x_len ) + error( 'pwelch: arg 2 (y) must be vector, same length as x.' ); + end + // force COLUMN vector + y = varargin(1)(:); + arg2_is_y = 1; + break; + end + end + // + // 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 ( ( type (arg) == 10 ) ) + // 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,'centerdc') ) + 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; + end + elseif ( ~strcmp(arg, 'cross' ) ) + if ( ~do_cross ) + n_results = n_results+1; + do_cross = n_results; + end + elseif ( ~strcmp(arg, 'trans' ) ) + if ( ~do_trans ) + n_results = n_results+1; + do_trans = n_results; + end + elseif ( ~strcmp(arg, 'coher' ) ) + if ( ~do_coher ) + n_results = n_results+1; + do_coher = n_results; + end + elseif ( ~strcmp(arg, 'ypower' ) ) + if ( ~do_ypower ) + n_results = n_results+1; + do_ypower = n_results; + end + else + error( 'pwelch: string arg %d illegal value: %s', iarg+1, arg ); + end + // 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' ); + end + // + // 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; + end + // + // -- "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; + end + // + // 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 = max(size(arg)); + if ( size(arg,2)>1 ) // vector must be COLUMN vector + arg = arg(:); + end + else + is_win = 0; + end + 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 and >=0',iarg+1); + end + 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 ); + end + 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 ); + end + Nfft = arg; + // + else + error( 'pwelch: arg %d must be string', iarg+1 ); + end + end + if ( conf>0 && (n_results && ~do_power ) ) + error('pwelch: can give confidence interval for x power spectrum only' ); + end + // + // 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; + end + 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 ); + end + if ( is_win > 1 ) + seg_len = min( max(size(window)), Nfft ); + window = window(1:seg_len); + else + if ( is_win ) + // window arg is scalar + seg_len = window; + else + seg_len = Nfft; + end + // make Hann window (don't depend on sigproc) + xx = seg_len - 1; + window = 0.5 - 0.5 * cos( (2*%pi/xx)*[0:xx].' ); + end + // + // Matlab R12 compatibility + elseif ( compatib==3 ) + if ( is_win > 1 ) + // window arg provides window function + seg_len = max(size(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 ); + end + // make Hamming window (don't depend on sigproc) + xx = seg_len - 1; + window = 0.54 - 0.46 * cos( (2*%pi/xx)*[0:xx].' ); + end + if ( isempty(Nfft) ) + Nfft = max( 256, 2^ceil(log(seg_len)*nearly_one/log_two) ); + end + // Matlab R14 psd(spectrum.welch) defaults + elseif ( compatib==4 ) + if ( is_win > 1 ) + // window arg provides window function + seg_len = max(size(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; + end + // make Hamming window (don't depend on sigproc) + xx = seg_len - 1; + window = 0.54 - 0.46 * cos( (2*%pi/xx)*[0:xx].' ); + end + // Now we know segment length, + // so we can set default overlap as number of samples + if ( ~isempty(overlap) ) + overlap = fix(seg_len * overlap / 100 ); + end + if ( isempty(Nfft) ) + Nfft = max( 256, 2^ceil(log(seg_len)*nearly_one/log_two) ); + end + // + // default compatibility level + else // if ( compatib==1 ) + // calculate/adjust segment lenght, window function + if ( is_win > 1 ) + // window arg provides window function + seg_len = max(size(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(max(size(x))) rounded up to nearest integer power of 2 + if ( isempty(overlap) ) + overlap=0.5; + end + seg_len = 2 ^ ceil( log(sqrt(x_len/(1-overlap)))*nearly_one/log_two ); + end + // make Hamming window (don't depend on sigproc) + xx = seg_len - 1; + window = 0.54 - 0.46 * cos( (2*%pi/xx)*[0:xx].' ); + end + // Now we know segment length, + // so we can set default overlap as number of samples + if ( ~isempty(overlap) ) + overlap = fix(seg_len * overlap); + end + // + // calculate FFT length + if ( isempty(Nfft) ) + Nfft = seg_len; + end + if ( is_sloppy ) + Nfft = 2 ^ ceil( log(Nfft) * nearly_one / log_two ); + end + end + // 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 <max(size(window)=%d',... + overlap, seg_len ); + end + // + // Pad data with zeros if shorter than segment. This should not happen. + if ( x_len < seg_len ) + x = [x; zeros(seg_len-x_len,1)]; + if ( arg2_is_y ) + y = [y; zeros(seg_len-x_len,1)]; + end + x_len = seg_len; + end + // end SETUP REMAINING PARAMETERS + // + // + // MAIN CALCULATIONS + // Remove mean from the data + if ( rm_mean == 3 ) + n_ffts = max( 0, fix( (x_len-seg_len)/(seg_len-overlap) ) ) + 1; + x_len = min( x_len, (seg_len-overlap)*(n_ffts-1)+seg_len ); + if ( need_Pxx || need_Pxy ) + x = x - sum( x(1:x_len) ) / x_len; + end + if ( arg2_is_y || need_Pxy) + y = y - sum( y(1:x_len) ) / x_len; + end + end + // + // Calculate and accumulate periodograms + // xx and yy are padded data segments + // Pxx, Pyy, Pyy are periodogram sums, Vxx is for confidence interval + xx = zeros(Nfft,1); + yy = xx; + Pxx = xx; + Pxy = xx; + Pyy = xx; + if ( conf>0 ) + Vxx = xx; + else + Vxx = []; + end + 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); + end + fft_x = fft1(xx); + end + 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); + end + fft_y = fft1(yy); + end + 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; + end + end + if ( need_Pxy ) + // Pxy (cross power spectrum) is complex. Do not force to be real. + Pxy = Pxy + fft_y .* conj(fft_x); + end + if ( need_Pyy ) + // force Pyy to be real + Pyy = Pyy + real(fft_y .* conj(fft_y)); + end + n_ffts = n_ffts +1; + end + // + // 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); + end + end + // + // 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 ) + if (modulo(Nfft,2) == 0 ) // 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]; + end + end + if ( need_Pxy ) + Pxy = Pxy(1:psd_len) + conj([0; Pxy(Nfft:-1:psd_len+1); 0]); + end + if ( need_Pyy ) + Pyy = Pyy(1:psd_len) + [0; Pyy(Nfft:-1:psd_len+1); 0]; + end + 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)]; + end + end + if ( need_Pxy ) + Pxy = Pxy(1:psd_len) + conj([0; Pxy(Nfft:-1:psd_len+1)]); + end + if ( need_Pyy ) + Pyy = Pyy(1:psd_len) + [0; Pyy(Nfft:-1:psd_len+1)]; + end + end + else // two-sided (and shifted) + psd_len = Nfft; + end + // 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; + end + end + if ( do_cross ) + spectra(:,do_cross) = Pxy / scale; + spect_type(do_cross) = 2; + end + if ( do_trans ) + spectra(:,do_trans) = Pxy ./ Pxx; + spect_type(do_trans) = 3; + end + if ( do_coher ) + // force coherence to be real + spectra(:,do_coher) = real(Pxy .* conj(Pxy)) ./ Pxx ./ Pyy; + spect_type(do_coher) = 4; + end + if ( do_ypower ) + spectra(:,do_ypower) = Pyy / scale; + spect_type(do_ypower) = 5; + end + 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,:)]; + end + end + // + // RETURN RESULTS or PLOT + // droping non postive values for plots + function res= postives(x) + j = 1 ;res = []; + for i=1:size(x,2) + if or( x(:,i) < 0 ) + j = j - 1 + warning("warning: axis: omitting non-positive data in log plot") + else + res(:,j) = x (:,i) + end + j = j + 1 + end + endfunction + // -------------------- + if ( nargout>=2 && conf>0 ) + varargout(2) = Vxx; + end + if ( nargout>=(2+(conf>0)) ) + // frequency is 2nd or 3rd return value, + // depends on if 2nd is confidence interval + varargout(2+(conf>0)) = freq; + end + 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 = []; + end + if ( n_results > 1 ) + figure(); + end + 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)) postives(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))]); + end + title( char(plot_title(spect_type(ii),:)) ); + ylabel( 'amplitude' ); + // 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))); + end + title( char(plot_title(spect_type(ii),:)) ); + ylabel( 'phase' ); + end + end + end end -end - -endfunction + endfunction + /* + +// demo 1 // use "hold on" and "hold off" for octave... + pi = %pi ; i = %i; + 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)); + skewed = signal.*exp(2*pi*i*2/25*[1:16384]); + compat = pwelch ([]); + hold on; + pwelch(skewed,[],[],[],Fs,'shift','semilogy'); + pwelch(skewed,[],[],[],Fs,0.95,'shift','semilogy'); + pwelch('R12+'); + pwelch(signal,'squared'); + pwelch (compat); + hold off; +// demo 2 // use "hold on" and "hold off" for octave... + a = [ 1.0 -1.6216505 1.1102795 -0.4621741 0.2075552 -0.018756746 ]; + white = rand(1,16384); + signal = detrend(filter(0.70181,a,white)); + 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 3 + 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 4 // use "hold on" and "hold off" for octave... + 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 5 + 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); + */ |