matlab开发-spectral

function varargout = spectral(x,dt,win,Nb,P,Wn,ftype,n) % Spectrum using periodogram method % % USAGE: % q = spectral(x,dt,win,Wn,ftype,n) % [psdf,f] = spectral(x,dt,win,Wn,ftype,n) % [psdf,psdfc,f] = spectral(x,dt,win,Wn,ftype,n) % % DESCRIPTION: % Calculates the spectrum for x % using the periodogram method % If a window other than boxcar is used % than the method is refered to as modified % periodogram method. % The confidence intervals are calculated % using the inverse of chi-square CDF. % Also includes a filtering option using the % butterworth filter to see the effect of the % filter on the spectrum % % INPUT VARIABLES: % x - Time series, [vector] % dt - Sampling Rate, [scalar] % win - Window, one of: % 'hanning', 'hamming', 'boxcar' % Nb - Band Averaging, number of bands to average % P - Probability for confidence intervals % Wn - Cut-Off frequencies, used for filtering % ftype - Type of filter, 'high', 'low' or 'stop' % ncb - Number of coefficients to use in % the Butterworth filter % % OUTPUT VARIABLES: % q - structure with the following fields: % xp - detrended x % f = Frequencies % T - Periods % m - Magnitude % a - Amplitude % s - Power spectrum, Sxx(win), [Power] % psdw - Power Spectral Density, Pxx(win), [Power/rad/sample] % psdf - Power Spectral Density, Pxx(f), [Power/sample-freq] % psdT - Power Spectral Density, Pxx(T), [Power*time-unit] % conf - Upper and Lower Confidence Interval multiplication % factors using chi-squared approach % %Copy-Left, Alejandro Sanchez if nargin<2 dt = 1; end if nargin<3 win = 'boxcar'; end if nargin<4 || isempty(Nb) Nb = 1; end if nargin<5 || isempty(P) P = 0.95; end if (nargin<7) && (nargin>5) ftype = 'low'; end if (nargin<8) && (nargin>5) n = 5; end if nargin==0 spectraldemo return end x = x(:); isf = isfinite(x); if nargout==0 figure t = dt*(0:(length(x)-1)); subplot(4,1,1) plot(t,x,'b') hold on t = t(isf); end x = x(isf); mx = mean(x); x = detrend(x); %also removes mean N = sum(isf); fn = (1/(2*dt)); %Nyquist frequency nfft = max(256, 2^(nextpow2(N))); f = (1:nfft/2)/(nfft/2)*fn; %frequencies, Note: starts from 1 f = f(:); %Frequencies %Makes a window try w = eval([win,'(',num2str(N),')']); catch error('Invalid Window') end %Butterworth Filter if nargin>5 [b1,b2] = butter(n,Wn*2*dt,ftype); %Wn/nyquist = Wn*2*dt x = filtfilt(b1,b2,x); xf = x; end if nargout==0 plot(t,x+mx,'r') axis tight set(gca,'xlim',[t(1),t(end)]) end q.xp = x; %save before applying window x = w.*x; %apply window % Evaluate the window normalization constant % A 1/nfft factor has been omitted since it will cancel below R = w'*w; % compute sum of squares of window %Compute Fourier Coefficients %Pad with zeros to length nfft, makes fft much faster %and increases frequency resolution, df=1/N/dt; c = fft(x,nfft); %The magnitude is multiplied by one of the following constants %to compensate for the lost energy when the window is applied %the values are taken from Emery & Thomson pag. 446-448 %In the case of the boxcar window the constant is 1 and is %therefore ommitted if strcmpi(win,'hanning') K = 1/2*sqrt(8/3); elseif strcmpi(win,'hamming') K = 1/2*sqrt(5/2); elseif strcmpi(win,'boxcar') K = 1; else error('Invalid Window') end %Magnitude m = K*abs(c(2:nfft/2+1)); %skip first coefficient which is the mean %Do band averaging now if Nb>1 f = bandaverage(f,Nb); m = bandaverage(m,Nb); end T = 1./f; %Periods %Amplitude a = m.*2./R; %Sxx(win), Angular-Frequency Power Spectrum, [Power] s = 2*(m.^2)./R; %Correction for one-sided spectrum s(2:end-1) = 2*s(2:end-1); %Pxx(win), Amgular-Frequency Spectrum %Power Spectrum Density (PSD), [Power/rad/sample] psdw = s/2/pi; %Pxx(f), Frequency Spectrum %Power Spectrum Density (PSD), [Power/sample-freq] psdf = s*dt; %Pxx(T), Period Spectrum %Power Spectrum Density (PSD), [Power*time-unit] psdT = 2*pi*psdw./T.^2; %Finally Calculate phase ph = unwrap(angle(c(1:fix(N/2)))); %Confidence Intervals using chi-squared approach v = 2*Nb; alpha2 = (1-P)/2; conf = v./chi2inv([alpha2,1-alpha2],v); %Collect variables switch nargout case 0 %do nothing case 1 varargout{1}.f = f; varargout{1}.T = T; varargout{1}.m = m; varargout{1}.a = a; varargout{1}.s = s; varargout{1}.psdw = psdw; varargout{1}.psdf = psdf; varargout{1}.psdT = psdT; varargout{1}.conf = conf; case 2 varargout{1} = psdf; varargout{2} = f; case 3 varargout{1} = psdf; varargout{2} = conf; varargout{3} = f; otherwise error('Invalid number of output variables') end %switch %Plotting if nargout==0 subplot(4,1,2) semilogx(T,a) set(gca,'xlim',[T(end),T(1)]) set(gca,'xticklabel',get(gca,'xtick')) xlabel('Period') ylabel('Amplitude') subplot(4,1,3) loglog(T,psdT) hold on loglog(T,psdT*conf,'--r') axis tight set(gca,'xticklabel',get(gca,'xtick')) xlabel('Period') ylabel('PSD(T)') %Filter or window subplot(4,1,4) if exist('Wn','var') H = freqz(b1,b2,f*pi*2*dt); semilogx(T,abs(H),'b'); axis tight set(gca,'xticklabel',get(gca,'xtick')) xlabel('Period') ylabel('Filter Response') else plot(w) ylabel('window') set(gca,'xticklabel',get(gca,'xtick')) axis tight % plot(ph) % xlabel('Period') end end return %---------------------------------------------------- function yy = bandaverage(y,Nb) %Does the band averaging Lb = floor(length(y)/Nb); yy = reshape(y(1:Nb*Lb),Nb,Lb); yy = mean(yy,1); y = y(:); return %------------------------------------------- function spectraldemo dt = 1/100; t=0:dt:200; sr = 2*randn(size(t)); %white noise s1 = 4*sin(2*pi*t/1); s2 = 2*sin(2*pi*t/4); s3 = 1*sin(t*2*pi/10); s4 = 3*sin(t*2*pi/20); s5 = -0.05*t; s0 = 5; x = s0 + s1 + s2 + s3 + s4 + s5 + sr; x(1000:6000) = NaN; Wn = 1; ftype = 'low'; win = 'hanning'; Nb = []; P = 0.95; n = 5; spectral(x,dt,win,Nb,P,Wn,ftype,n) return