gusucode.com > 信号处理工具箱 - signal源码程序 > signal\signal\signal\pmtm.m
function varargout=pmtm(x,varargin); %PMTM Power Spectrum estimate via the Thomson multitaper method. % Pxx = PMTM(X) is the Power Spectral Density (PSD) estimate, Pxx(w), % of time series X using MTM (multitaper method). % % The single-sided PSD is returned for real signals and the % double-sided PSD is returned for complex signals. % % Pxx = PMTM(X,NW,NFFT) uses NW as the "time-bandwidth product" for % the discrete prolate spheroidal sequences (or Slepian sequences) % used as data windows. Typical choices for NW are 2, 5/2, 3, 7/2, or % 4 (the default is 4). The number of sequences used to form Pxx is % 2*NW-1. NFFT defines the frequency grid NFFT. When X is real, Pxx % is length (NFFT/2+1) for NFFT even or (NFFT+1)/2 for NFFT odd; when % X is complex, Pxx is length NFFT. NFFT is optional; it defaults to % the greater of 256 and the next power of 2 greater than the length % of X. % % [Pxx,W] = PMTM(X,NW,NFFT) returns the frequency vector, W, in % rads/sample at which the PSD is estimated. The PSD estimate is % computed over the interval [0, Pi] for a real signal X and over the % interval [0, 2*Pi] for a complex X. % % [Pxx,F] = PMTM(X,NW,NFFT,Fs) return the PSD estimate and the vector % of frequencies, F, in Hz, at which the PSD is estimated. Fs is the % sampling frequency. In this case, the PSD estimate is computed over % the interval [0, Fs/2] for a real signal X and over the interval % [0, Fs] for a complex X. If left empty, Fs defaults to 1 Hz. % % [Pxx,F] = PMTM(X,NW,NFFT,Fs,method) uses the algorithm in method for % combining the individual spectral estimates: % 'adapt' - Thomson's adaptive non-linear combination (default). % 'unity' - linear combination with unity weights. % 'eigen' - linear combination with eigenvalue weights. % % [Pxx,Pxxc,F] = PMTM(X,NW,NFFT,Fs,method) returns the 95% confidence % interval Pxxc for Pxx. [Pxx,Pxxc,F] = PMTM(X,NW,NFFT,Fs,method,P) % where P is a scalar between 0 and 1, returns the P*100% confidence % interval for Pxx. Confidence intervals are computed using a % chi-squared approach. Pxxc(:,1) is the lower bound of the confidence % interval, Pxxc(:,2) is the upper bound. % % PMTM with no output arguments plots the PSD in the current or next % available figure, with confidence intervals. % % [Pxx,Pxxc,F] = PMTM(X,E,V,NFFT,Fs,method,P) is the PSD estimate, % confidence interval, and frequency vector from the data tapers in E % and their concentrations V. % % [Pxx,Pxxc,F] = PMTM(X,DPSS_PARAMS,NFFT,Fs,method,P) is the PSD % estimate, confidence interval, and frequency vector from the data % tapers computed with the function DPSS with parameters in the cell % array DPSS_PARAMS. The elements of DPSS_PARAMS contain the % parameters to DPSS starting with the second one. For example, % PMTM(X,{3.5,'trace'},512,Fs) calculates the Slepian sequences and % returns the method DPSS uses. See HELP DPSS for options. % % You can obtain default parameters by inserting an empty matrix [], % e.g. PMTM(X,[],[],10000) uses defaults for NW and NFFT. % % EXAMPLE % Fs = 1000; t = 0:1/Fs:.3; x = cos(2*pi*t*200)+randn(size(t)); % [Pxx,Pxxc,f] = pmtm(x,3.5,512,Fs,[],.99); % plot(f,10*log10([Pxx Pxxc])) % % See also DPSS, PWELCH, PMUSIC, PBURG, PYULEAR, PCOV, PEIG and PMCOV. % Author: Eric Breitenberger, version date 10/1/95. % Copyright (c) 1988-98 by The MathWorks, Inc. % $Revision: 1.6.1.2 $ $Date: 1999/01/22 03:42:33 $ % defaults: N=length(x); NW = 4; nfft = max(256,2^nextpow2(N)); Fs = 1; method = 'adapt'; conf = .95; range = 'half'; samprateflag = 0; magunits = 'db'; if length(varargin)>0 if length(varargin{1})>0 NW = varargin{1}; end if iscell(NW) [E,V]=dpss(N,NW{:}); NW = NW{1}; elseif length(NW)>1 E = NW; if length(varargin)<2 error('Must provide V with E matrix.') else V = varargin{2}; end varargin(2) = []; if size(E,2)~=length(V) error('Number of columns of E and length of V do not match.') end NW = size(E,2)/2; % only used for computation of # of tapers k else % Get the dpss, one way or another: [E,V] = dpss(N,NW); end else % Get the dpss, one way or another: [E,V] = dpss(N,NW); end if length(varargin)>1 & ~isempty(varargin{2}), nfft = varargin{2}; end if length(varargin)>2 samprateflag = 1; if ~isempty(varargin{3}), Fs = varargin{3}; end end if length(varargin)>3 & ~isempty(varargin{4}), method = lower(varargin{4}); end if length(varargin)>4 & ~isempty(varargin{5}), conf = varargin{5}; end % error checking: switch method case {'adapt','unity','eigen'} % do nothing otherwise error('Method must be ''adapt'',''unity'', or ''eigen''.') end if ~all(size(conf)==[1 1]) error('Confidence interval parameter P must be a scalar.') end x=x(:); k=min(round(2*NW),N); % By convention, the first 2*NW % eigenvalues/vectors are stored k=max(k-1,1); V=V(1:k); % Compute the windowed dfts and the % corresponding spectral estimates: if N<=nfft Sk=abs(fft(E(:,1:k).*x(:,ones(1,k)),nfft)).^2; else % use CZT to compute DFT on nfft evenly spaced samples around the % unit circle: Sk=abs(czt(E(:,1:k).*x(:,ones(1,k)),nfft)).^2; end % Select the proper points from fft: if isreal(x) if rem(nfft,2)==0, M=nfft/2+1; else M=(nfft+1)/2; end else range = 'whole'; M = nfft; end Sk=Sk(1:M,:); switch method case 'adapt' % Set up the iteration to determine the adaptive weights: sig2=x'*x/N; % Power S=(Sk(:,1)+Sk(:,2))/2; % Initial spectrum estimate Stemp=zeros(M,1); S1=zeros(M,1); % The algorithm converges so fast that results are % usually 'indistinguishable' after about three iterations. % This version uses the equations from P&W pp 368-370 % Set tolerance for acceptance of spectral estimate: tol=.0005*sig2/M; i=0; a=sig2*(1-V); % Do the iteration: while sum(abs(S-S1)/M)>tol i=i+1; % calculate weights b=(S*ones(1,k))./(S*V'+ones(M,1)*a'); % calculate new spectral estimate wk=(b.^2).*(ones(M,1)*V'); S1=sum(wk'.*Sk')./ sum(wk'); S1=S1'; Stemp=S1; S1=S; S=Stemp; % swap S and S1 end case {'unity','eigen'} % Compute the averaged estimate: simple arithmetic % averaging is used. The Sk can also be weighted % by the eigenvalues, as in Park et al. Eqn. 9.; % note that that eqn. apparently has a typo. as % the weights should be V and not 1/V. if strcmp(method,'eigen') wt = V(:); % Park estimate else wt = ones(k,1); end S = Sk*wt/k; end % Scale by 1/Fs or 1/2pi in order to get Power per unit of frequency if samprateflag, % Linear freq (Hz) specified scaleFactor = Fs; else % Using default normalized, angular freq scaleFactor = 2*pi; end S = S ./ scaleFactor; % For real signals return the single-sided spectrum with full power. if isreal(x), S = [S(1); 2*S(2:end-1); S(end)]; end % Calculate confidence limits if nargout==0 | nargout==3 % don't calculate these unless needed, since it can take a while. c=S*chi2conf(conf,k); end % Default frequency vector is normalized angular frequency; either [0, pi) % or [0,2*pi). If user specifies Fs, linear frequency in Hz is returned. [ff,xlab,xtickFlag,xlim] = calcfreqvector(M,Fs,samprateflag,range); if nargout == 0 % do plots newplot; pxxplot('MTM',S,ff,range,xlab,xtickFlag,xlim,magunits); end if nargout > 0, varargout{1} = S; end if nargout > 1, varargout{2} = ff; end if nargout > 2, varargout{2} = c; varargout{3} = ff; end % [EOF] pmtm.m