m_spec5.m0000755000175000001440000003543210325242102011263 0ustar nghusersfunction m_spec5 close all hidden %% Use system background color for GUI components panelColor = get(0,'DefaultUicontrolBackgroundColor'); %% ------------ Callback Functions --------------- % Figure resize function function figResizeC(src,evt) fpos = get(fC,'Position'); for i=1:2 set(filePanelC(i),'Position',... [10/20+(i-1)*fpos(3)/2 fpos(4)*10/35 fpos(3)/2-.1 fpos(4)*25/35]) end set(botPanelC,'Position',... [1/20 1/20 fpos(3)-.1 fpos(4)*10/35]) end function figResizeP(src,evt) fpos = get(fP,'Position'); for i=1:2 set(botPanelP(i),'Position',... [1/20+(i-1)*fpos(3)/2 1/20 fpos(3)/2-.1 fpos(4)*10/35]) set(specPanel(i),'Position',... [1/20+(i-1)*fpos(3)/2 fpos(4)*10/35 fpos(3)/2-.1 fpos(4)*25/35]); end end function figResizeM(src,evt) fpos = get(fM,'Position'); for i=1:2 set(mapPanel(i),'Position',... [.1+(i-1)*fpos(3)/2 .1 fpos(3)/2 fpos(4)*.75-.1]); set(mapText(i),'Position',... [.2+(i-1)*fpos(3)/2 .2+fpos(4)*.75 fpos(3)/2-.4 fpos(4)*.25-.4]); end end % Bottom panel resize function function botPanelCResize(src, evt) bpos = get(botPanelC,'Position'); set(plotButton,'Position',... [bpos(3)*4/120 bpos(4)*6/8 bpos(3)*40/120 1.5]) set(EqualAxesToggle,'Position',... [bpos(3)*4/120 bpos(4)*5/8 bpos(3)*40/120 1.5]) set(RotaryToggle,'Position',... [bpos(3)*64/120 bpos(4)*5/8 bpos(3)*40/120 1.5]) set(popUpMap,'Position',... [bpos(3)*64/120 bpos(4)*6/8 bpos(3)*40/120 1.5]) set(popUpPlotType,'Position',... [bpos(3)*4/120 bpos(4)*3/8 bpos(3)*40/120 1.5]) set(popUpPlotTypeLabel,'Position',... [bpos(3)*4/120 bpos(4)*3.5/8 bpos(3)*40/120 2]) set(popUpPcLngth,'Position',... [bpos(3)*64/120 bpos(4)*3/8 bpos(3)*40/120 1.5]) set(popUpPcLngthLabel,'Position',... [bpos(3)*64/120 bpos(4)*3.5/8 bpos(3)*40/120 2]) end % Right panel resize function function filePanelCResize(src,evt,j) rpos = get(filePanelC(i),'Position'); set(setSpecLabel(i),'Position',... [rpos(3)*4/32 rpos(4)*25/27 rpos(3)*24/32 rpos(4)*1/27]); set(setSpec(i),'Position',... [rpos(3)*4/32 rpos(4)*24/27 rpos(3)*24/32 rpos(4)*1/27]); for i=1:2 set(listFile(i,j),'Position',... [rpos(3)*4/32 rpos(4)*(17-(i-1)*7)/27 rpos(3)*24/32 rpos(4)*5/27]); set(listFileLabel(i,j),'Position',... [rpos(3)*4/32 rpos(4)*(22-(i-1)*7)/27 rpos(3)*24/32 rpos(4)*1/27]); set(listDepth(i,j),'Position',... [rpos(3)*(3+(i-1)*13)/32 rpos(4)*5/27 rpos(3)*13/32 rpos(4)*3/27]); set(listDepthLabel(i,j),'Position',... [rpos(3)*(4+(i-1)*12)/32 rpos(4)*7/27 rpos(3)*12/32 rpos(4)*2/27]); set(listVar(i,j),'Position',... [rpos(3)*(4+(i-1)*12)/32 rpos(4)*1/27 rpos(3)*12/32 rpos(4)*2/27]); set(listVarLabel(i,j),'Position',... [rpos(3)*(4+(i-1)*12)/32 rpos(4)*3/27 rpos(3)*12/32 rpos(4)*1.5/27]); end end %% Callback for listing files after file spec change function setSpecCallback(src,evt,iP) for k=1:2 listFileCallback(listFile(k,iP),iP) listDepthCallback(listFile(k,iP),[],k,iP) end end % setSpecCallback %% Callback for list box function listFileCallback(src,evt) % Load workspace vars into list box fspc=get(setSpec(evt),'String'); fspc=strread(fspc,'%s'); vars=[]; for i=1:length(fspc) vars = [vars;dir([fspc{i},'.mat'])]; end for i=1:length(vars) varnames=vars(i).name; [dn,fn{i}]=fileparts(varnames); end set(src,'String',fn) end % listFileCallback %% Callback for depth choice buttons function listDepthCallback(src,evt,iD,iP) % Get workspace variables vars = get(src,'String'); var_index = get(src,'Value'); fi = load(vars{var_index}); depths={}; for i=1:length(fi.cmdat) depths{i}=int2str(fi.cmdat(i).depth); end val=min([length(depths),get(listDepth(iD,iP),'Value')]); set(listDepth(iD,iP),'String',depths,'Value',val) end % listDepthCallback %% Callback for plot button function plotButtonCallback(src,evt) % Get workspace variables plt_type=get(popUpPlotType,'String'); plt_type_index=get(popUpPlotType,'Value'); pc_lngth=get(popUpPcLngth,'String'); pc_lngth_index=get(popUpPcLngth,'Value'); axes_scale=get(EqualAxesToggle,'Value'); clr={'r','g';'m','c'};h=[];yl=[]; rot_ind=get(RotaryToggle,'Value'); for j=1:2 figure(fM) if ~isempty(get(mapPanel(j),'children')) delete(get(mapPanel(j),'children')) end cms = get(listFile(1,j),'String'); axes('parent',mapPanel(j)) cla mlat=[];mlon=[];ub=[];vb=[];mrg={}; for i=1:length(cms) cmstr=load(cms{i}); mlat(i)=cmstr.cmdat(1).lat; mlon(i)=cmstr.cmdat(1).lon; ub(i)=nanmean(cmstr.cmdat(1).u); vb(i)=nanmean(cmstr.cmdat(1).v); end lats=[min(mlat)-1,max(mlat)+1];lons=[min(mlon)-1,max(mlon)+1]; m_proj('mercator','lon',lons,'lat',lats); m_grid('linestyle','none','tickdir','out','linewidth',1) hold for i=1:length(cms) m_text(mlon(i),mlat(i),cms{i}(~isletter(cms{i}))); end m_quiver(mlon,mlat,ub,vb); m_plot(mlon,mlat,'rx'); cms_index = get(listFile(1,j),'Value'); slon=[];slat=[]; for i=1:2 cms = get(listFile(i,j),'String'); cms_index = get(listFile(i,j),'Value'); mrg{i}=cms{cms_index}; cmstr = load(cms{cms_index}); insts=get(listDepth(i,j),'String'); insts_index = get(listDepth(i,j),'Value'); inst{i}=insts{insts_index}; cm = cmstr.cmdat(insts_index); slon=[slon,cm.lon];slat=[slat,cm.lat]; m_plot(cm.lon,cm.lat,'g*'); tm{i}=cm.tm; var_index = get(listVar(i,j),'Value'); vars=get(listVar(i,j),'String'); tit{i}=[mrg{i},' ',vars{var_index},' ',inst{i}]; switch var_index case 1 if ~rot_ind var{i}=cm.u; else var{i}=complex(cm.u,-cm.v)/sqrt(2); end case 2 if ~rot_ind var{i}=cm.v; else var{i}=complex(cm.u,cm.v)/sqrt(2); end case 3 var{i}=cm.temp; end end set(mapText(j),'string',[num2str(sw_dist(slat,slon,'km'),'%.1f'),'km apart']) %% Add an axes to the plot panels figure(fP) axes('parent',botPanelP(j),'position',[.05,.1,.9,.85]); %cla reset h=plot(tm{1},var{1},clr{1,j},tm{2},var{2},clr{2,j}); %hold on legend(h,tit) datetick switch plt_type{plt_type_index} case 'Log-Log' if exist('a'),set(a,'visible','off'),end [ha(j),hp]=comp_spec(tm,var,str2num(pc_lngth{pc_lngth_index}),specPanel(j)); case 'Var Pres' end legend(ha(j),hp,tit,'location','southwest') end if ~axes_scale yl=get(ha,'ylim'); set(ha,'ylim',[min([yl{1}(1),yl{2}(1)]),max([yl{1}(2),yl{2}(2)])]) end end % plotButtonCallback %% Callback for map plot state toggle button function popUpMapCallback(src,evt) button_state = get(src,'Value'); % if button_state > 1 % for i=1:2 % axes(mapPanel(i)) % end end % popUpMapCallback %% Callback for axes scale state toggle button function EqualAxesToggleCallback(src,evt) button_state = get(src,'Value'); if button_state == get(src,'Max') % toggle button is depressed set(src,'String','Unequal Axes') elseif button_state == get(src,'Min') % toggle button is not depressed set(src,'String','Equal Axes') end end % EqualAxesToggleCallback %% Callback for cart/rotary toggle button function RotaryToggleCallback(src,evt) button_state = get(src,'Value'); if button_state == get(src,'Min') % toggle button is depressed set(src,'String','Cartesian') for i=1:2 for j=1:2 set(listVar(i,j),'String',{'U','V','T'}) end end elseif button_state == get(src,'Max') % toggle button is not depressed set(src,'String','Rotary') for i=1:2 for j=1:2 set(listVar(i,j),'String',{'CW','ACW','T'}) end end end end % RotaryToggleCallback %% ------------ GUI layout --------------- %% Set up the Plot figure and defaults fP = figure('Units','characters',... 'Position',[80 5 120 55],... 'PaperPositionMode','auto',... 'Color',panelColor,... 'HandleVisibility','callback',... 'IntegerHandle','off',... 'Renderer','painters',... 'Toolbar','figure',... 'NumberTitle','off',... 'Name','Time Series Analyzer',... 'ResizeFcn',@figResizeP); %% Create the bottom and center uipanels for i=1:2 botPanelP(i) = uipanel('BorderType','etchedin',... 'BackgroundColor',panelColor,... 'Units','characters',... 'Position',[1/20 1/20 119.9 8],... 'Parent',fP); specPanel(i) = uipanel('bordertype','etchedin',... 'Units','characters',... 'Position', [1/20+44*(i-1) 8 44 27],... 'Parent',fP); end %% Create mapping figure fM = figure('Units','characters',... 'Position',[180 5 50 15],... 'PaperPositionMode','auto',... 'Color',panelColor,... 'HandleVisibility','callback',... 'IntegerHandle','off',... 'Renderer','painters',... 'Toolbar','figure',... 'NumberTitle','off',... 'Name','Chart',... 'ResizeFcn',@figResizeM); %% Create the two mapping uipanels for i=1:2 mapPanel(i) = uipanel('BorderType','etchedin',... 'BackgroundColor','white',... 'Units','characters',... 'Position',[1/20 1/20 119.9 8],... 'Parent',fM); mapText(i) = uicontrol('Style','text','Units','characters',... 'Position',[1 1 20 4],'String',['Map parms'],... 'BackgroundColor',panelColor,'Parent',fM); end %% Set up the Control figure and defaults fC = figure('Units','characters',... 'Position',[3 5 50 55],'Color',panelColor,... 'HandleVisibility','callback',... 'IntegerHandle','off','Renderer','painters',... 'Toolbar','figure','NumberTitle','off',... 'Name','TS Control','ResizeFcn',@figResizeC); %% Create the right side panel for i=1:2 filePanelC(i) = uipanel('bordertype','etchedin',... 'BackgroundColor',panelColor,'Units','characters',... 'Position',[88 8 32 27],'Parent',fC,... 'ResizeFcn',{@filePanelCResize,i}); end %% Create the bottom uipanel botPanelC = uipanel('BorderType','etchedin',... 'BackgroundColor',panelColor,'Units','characters',... 'Position',[1/20 1/20 119.9 8],'Parent',fC,... 'ResizeFcn',@botPanelCResize); %% Add listbox and label for i=1:2 setSpecLabel(i) = uicontrol(fC,'Style','text','Units','characters',... 'Position',[4 24 24 1],'String',['Spec ',int2str(i)],... 'BackgroundColor',panelColor,'Parent',filePanelC(i)); setSpec(i) = uicontrol(fC,'Style','edit','Units','characters',... 'Position',[4 24 24 1],'String','*f','BackgroundColor','white',... 'Parent',filePanelC(i),'Callback',{@setSpecCallback,i}); for j=1:2 listFileLabel(i,j) = uicontrol(fC,'Style','text','Units','characters',... 'Position',[4 24 24 1],'String',['File ',int2str(i)],... 'BackgroundColor',panelColor,'Parent',filePanelC(j)); listFile(i,j) = uicontrol(fC,'Style','listbox','Units','characters',... 'Position',[4 12 24 10],'BackgroundColor','white',... 'Max',1,'Min',1,'Parent',filePanelC(j),... 'Callback',{@listDepthCallback,i,j}); listDepthLabel(i,j) = uicontrol(fC,'Style','text','Units','characters',... 'Position',[4 6 12 4],'String',['Depth ',int2str(i)],... 'BackgroundColor',panelColor,'Parent',filePanelC(j)); listDepth(i,j) = uicontrol(fC,'Style','listbox','Units','characters',... 'Position',[4 2 12 6],'BackgroundColor','white',... 'Max',1,'Min',1,'Parent',filePanelC(j)); listVar(i,j) = uicontrol(fC,'Style','listbox','Units','characters',... 'Position',[4 2 12 6],'BackgroundColor','white',... 'Max',1,'Min',1,'Parent',filePanelC(j),... 'String',{'U','V','T'}); listVarLabel(i,j) = uicontrol(fC,'Style','text','Units','characters',... 'Position',[4 6 12 4],'String',['Var ',int2str(i)],... 'BackgroundColor',panelColor,'Parent',filePanelC(j)); end end %% Add popups and labels popUpPlotTypeLabel = uicontrol(fC,'Style','text','Units','characters',... 'Position',[10 14 24 2],... 'String','Plot Type',... 'BackgroundColor',panelColor,... 'Parent',botPanelC); popUpPlotType = uicontrol(fC,'Style','popupmenu','Units','characters',... 'Position',[10 12 24 2],... 'BackgroundColor','white',... 'String',{'Log-Log','Var Pres'},... 'Parent',botPanelC); popUpPcLngthLabel = uicontrol(fC,'Style','text','Units','characters',... 'Position',[80 4 24 2],... 'String','Piece Length',... 'BackgroundColor',panelColor,... 'Parent',botPanelC); popUpPcLngth = uicontrol(fC,'Style','popupmenu','Units','characters',... 'Position',[80 2 24 2],... 'BackgroundColor','white',... 'String',{'16','32','64','128','256','512'},... 'Parent',botPanelC); %% Add mapPlot and plot buttons popUpMap = uicontrol(fC,'Style','popupmenu','Units','characters',... 'Position',[4 2 12 2],... 'String',{'No Maps','1 Map','2 Maps'}',... 'Parent',botPanelC,... 'Callback',{@popUpMapCallback,fM}); EqualAxesToggle = uicontrol(fC,'Style','toggle','Units','characters',... 'Position',[4 2 12 2],... 'String','Axes Scales',... 'Parent',botPanelC,... 'Callback',@EqualAxesToggleCallback); RotaryToggle = uicontrol(fC,'Style','toggle','Units','characters',... 'Position',[4 2 12 2],... 'String','Rotary/Cart.',... 'Parent',botPanelC,... 'Callback',@RotaryToggleCallback); plotButton = uicontrol(fC,'Style','pushbutton','Units','characters',... 'Position',[10 2 24 2],... 'String','Plot',... 'Parent',botPanelC,... 'Callback',@plotButtonCallback); %% Initialize list box and make sure for k=1:2 for iP=1:2 listFileCallback(listFile(k,iP),iP) listDepthCallback(listFile(k,iP),[],k,iP) end end end % uipanel1 comp_spec.m0000755000175000001440000000334110323472732011706 0ustar nghusersfunction [h1a,hl]=comp_spec(tm,v,len,cpan) conf_level=.95; if ~isempty(get(cpan,'child')),delete(get(cpan,'child')),end [tm,v,fs]=fix_timebase(tm,v); %hc=get(gcf,'children');%gets rid of present plots %delete(hc(strcmp(get(hc,'type'),'axes'))) len2=len/2; [P,f,ks]=my_spectrum(v{1},v{2},len,len2,[],fs); if any(~isreal(v{1}))%rotary spectrum ng=1+size(P,1)/2; P=P(1:ng,:);f=f(1:ng); end P(1,:)=[];f(1)=[];coh=sqrt(P(:,5));phs=angle(P(:,3))*180/pi; %[P,f,k(:,i)]=freqav(P(2:(len/2+1),:),f(2:(len/2+1)),ks,nfb); h1a=axes('parent',cpan,'position',[.05,.5,.9,.45]); mxy=max(P(:,1:2)); mny=min(P(:,1:2)); yl=[10^(floor(log10(min(mny)))),10^(ceil(log10(max(mxy))))]; hl=loglog(f,P(:,1),'r',f,P(:,2),'--g'); pos=get(h1a,'position');pos(1)=pos(1)+.05;pos(3)=pos(3)-.05; set(h1a,'xticklabel',[],'ylim',yl,'position',pos) h2a=axes('parent',cpan,'position',[.05,.3,.9,.2]); err=confphs(coh,ks); conf=confcoh(ks); nz=find(coh>conf); % keyboard semilogx(f,phs,'g'), if ~isempty(nz) errbarx(f(nz),phs(nz),err(nz)) end xl=get(h2a,'xlim'); hold on, semilogx(xl,[0,0],'r:',xl,[-90,-90],'r:',xl,[90,90],'r:'), hold off pos=get(h2a,'position');pos(1)=pos(1)+.05;pos(3)=pos(3)-.05; set(h2a,'ylim',[-180,180],'ytick',[-180,-90,0,90],'xticklabel',[],'position',pos); h3a=axes('parent',cpan,'position',[.05,.1,.9,.2]); semilogx(f,coh,'g',xl,[conf,conf],'--') pos=get(h3a,'position');pos(1)=pos(1)+.05;pos(3)=pos(3)-.05; set(h3a,'position',pos); %now do confidence limits np=fix(.9*length(f)); if iscell(yl) confsp=confspec(ks,conf_level)*yl{2}(2)/10; else confsp=confspec(ks,conf_level)*yl(2)/10; end axes(h1a); hold on plot(f(np),confsp(2),'+',[f(np),f(np)],[confsp(1),confsp(3)]) text(f(np),confsp(2),[' ',int2str(fix(100*conf_level)),'%']); hold off toolbox/mspec/my_spectrum.m0000755000175000001440000003166410321232073015101 0ustar nghusersfunction [Spec,f,Np] = my_spectrum(varargin) %SPECTRUM Power spectrum estimate of one or two data sequences. % P=SPECTRUM(X,NFFT,NOVERLAP,WIND) estimates the Power Spectral Density of % signal vector X using Welch's averaged periodogram method. The signal X % is divided into overlapping sections, each of which is detrended and % windowed by the WINDOW parameter, then zero padded to length NFFT. The % magnitude squared of the length NFFT DFTs of the sections are averaged % to form Pxx. P is a two column matrix P = [Pxx Pxxc]; the second column % Pxxc is the 95% confidence interval. The number of rows of P is NFFT/2+1 % for NFFT even, (NFFT+1)/2 for NFFT odd, or NFFT if the signal X is comp- % lex. If you specify a scalar for WINDOW, a Hanning window of that length % is used. % % [P,F] = SPECTRUM(X,NFFT,NOVERLAP,WINDOW,Fs) given a sampling frequency % Fs returns a vector of frequencies the same length as Pxx at which the % PSD is estimated. PLOT(F,P(:,1)) plots the power spectrum estimate % versus true frequency. % % [P, F] = SPECTRUM(X,NFFT,NOVERLAP,WINDOW,Fs,Pr) where Pr is a scalar % between 0 and 1, overrides the default 95% confidence interval and % returns the Pr*100% confidence interval for Pxx instead. % % SPECTRUM(X) with no output arguments plots the PSD in the current % figure window, with confidence intervals. % % The default values for the parameters are NFFT = 256 (or LENGTH(X), % whichever is smaller), NOVERLAP = 0, WINDOW = HANNING(NFFT), Fs = 2, % and Pr = .95. You can obtain a default parameter by leaving it out % or inserting an empty matrix [], e.g. SPECTRUM(X,[],128). % % P = SPECTRUM(X,Y) performs spectral analysis of the two sequences % X and Y using the Welch method. SPECTRUM returns the 8 column array % P = [Pxx Pyy Pxy Txy Cxy Pxxc Pyyc Pxyc] % where % Pxx = X-vector power spectral density % Pyy = Y-vector power spectral density % Pxy = Cross spectral density % Txy = Complex transfer function from X to Y = Pxy./Pxx % Cxy = Coherence function between X and Y = (abs(Pxy).^2)./(Pxx.*Pyy) % Pxxc,Pyyc,Pxyc = Confidence range. % All input and output options are otherwise exactly the same as for the % single input case. % % SPECTRUM(X,Y) with no output arguments will plot Pxx, Pyy, abs(Txy), % angle(Txy) and Cxy in sequence, pausing between plots. % % SPECTRUM(X,...,DFLAG), where DFLAG can be 'linear', 'mean' or 'none', % specifies a detrending mode for the prewindowed sections of X (and Y). % DFLAG can take the place of any parameter in the parameter list % (besides X) as long as it is last, e.g. SPECTRUM(X,'none'); % % See also PSD, CSD, TFE, COHERE, SPECGRAM, SPECPLOT, DETREND, PMTM, % PMUSIC. % ETFE, SPA, and ARX in the Identification Toolbox. % Author(s): J.N. Little, 7-9-86 % C. Denham, 4-25-88, revised % L. Shure, 12-20-88, revised % J.N. Little, 8-31-89, revised % L. Shure, 8-11-92, revised % T. Krauss, 4-15-93, revised % Copyright (c) 1988-1999 The MathWorks, Inc. All Rights Reserved. % $Revision: 1.3 $ $Date: 1999/01/11 19:05:13 $ % The units on the power spectra Pxx and Pyy are such that, using % Parseval's theorem: % % SUM(Pxx)/LENGTH(Pxx) = SUM(X.^2)/LENGTH(X) = COV(X) % % The RMS value of the signal is the square root of this. % If the input signal is in Volts as a function of time, then % the units on Pxx are Volts^2*seconds = Volt^2/Hz. % % Here are the covariance, RMS, and spectral amplitude values of % some common functions: % Function Cov=SUM(Pxx)/LENGTH(Pxx) RMS Pxx % a*sin(w*t) a^2/2 a/sqrt(2) a^2*LENGTH(Pxx)/4 %Normal: a*rand(t) a^2 a a^2 %Uniform: a*rand(t) a^2/12 a/sqrt(12) a^2/12 % % For example, a pure sine wave with amplitude A has an RMS value % of A/sqrt(2), so A = SQRT(2*SUM(Pxx)/LENGTH(Pxx)). % % See Page 556, A.V. Oppenheim and R.W. Schafer, Digital Signal % Processing, Prentice-Hall, 1975. error(nargchk(1,8,nargin)) [msg,x,y,nfft,noverlap,window,Fs,p,dflag]=specchk(varargin); error(msg) if isempty(p), p = .95; % default confidence interval even if not asked for end n = length(x); % Number of data points nwind = length(window); if n < nwind % zero-pad x (and y) if length less than the window length x(nwind)=0; n=nwind; if ~isempty(y), y(nwind)=0; end end x = x(:); % Make sure x and y are column vectors y = y(:); k = fix((n-noverlap)/(nwind-noverlap)); % Number of windows % (k = fix(n/nwind) for noverlap=0) index = 1:nwind; KMU = k*norm(window)^2; % Normalizing scale factor ==> asymptotically unbiased % KMU = k*sum(window)^2;% alt. Nrmlzng scale factor ==> peaks are about right if (isempty(y)) % Single sequence case. Pxx = zeros(nfft,1); Pxx2 = zeros(nfft,1); for i=1:k if strcmp(dflag,'linear') xw = window.*detrend(x(index)); elseif strcmp(dflag,'none') xw = window.*(x(index)); else xw = window.*detrend(x(index),0); end index = index + (nwind - noverlap); Xx = abs(fft(xw,nfft)).^2; Pxx = Pxx + Xx; Pxx2 = Pxx2 + abs(Xx).^2; end % Select first half if ~any(any(imag(x)~=0)), % if x and y are not complex if rem(nfft,2), % nfft odd select = [1:(nfft+1)/2]; else select = [1:nfft/2+1]; % include DC AND Nyquist end else select = 1:nfft; end Pxx = Pxx(select); Pxx2 = Pxx2(select); cPxx = zeros(size(Pxx)); if k > 1 c = (k.*Pxx2-abs(Pxx).^2)./(k-1); c = max(c,zeros(size(Pxx))); cPxx = sqrt(c); end ff = sqrt(2)*erfinv(p); % Equal-tails. Pxx = Pxx/KMU; Pxxc = ff.*cPxx/KMU; P = [Pxx Pxxc]; else Pxx = zeros(nfft,1); % Dual sequence case. Pyy = Pxx; Pxy = Pxx; Pxx2 = Pxx; Pyy2 = Pxx; Pxy2 = Pxx; for i=1:k if strcmp(dflag,'linear') xw = window.*detrend(x(index)); yw = window.*detrend(y(index)); elseif strcmp(dflag,'none') xw = window.*(x(index)); yw = window.*(y(index)); else xw = window.*detrend(x(index),0); yw = window.*detrend(y(index),0); end index = index + (nwind - noverlap); Xx = fft(xw,nfft); Yy = fft(yw,nfft); Yy2 = abs(Yy).^2; Xx2 = abs(Xx).^2; Xy = Yy .* conj(Xx); Pxx = Pxx + Xx2; Pyy = Pyy + Yy2; Pxy = Pxy + Xy; Pxx2 = Pxx2 + abs(Xx2).^2; Pyy2 = Pyy2 + abs(Yy2).^2; Pxy2 = Pxy2 + Xy .* conj(Xy); end % Select first half if ~any(any(imag([x y])~=0)), % if x and y are not complex if rem(nfft,2), % nfft odd select = [1:(nfft+1)/2]; else select = [1:nfft/2+1]; % include DC AND Nyquist end else select = 1:nfft; end Pxx = Pxx(select); Pyy = Pyy(select); Pxy = Pxy(select); Pxx2 = Pxx2(select); Pyy2 = Pyy2(select); Pxy2 = Pxy2(select); cPxx = zeros(size(Pxx)); cPyy = cPxx; cPxy = cPxx; if k > 1 c = max((k.*Pxx2-abs(Pxx).^2)./(k-1),zeros(size(Pxx))); cPxx = sqrt(c); c = max((k.*Pyy2-abs(Pyy).^2)./(k-1),zeros(size(Pxx))); cPyy = sqrt(c); c = max((k.*Pxy2-abs(Pxy).^2)./(k-1),zeros(size(Pxx))); cPxy = sqrt(c); end Txy = Pxy./Pxx; Cxy = (abs(Pxy).^2)./(Pxx.*Pyy); ff = sqrt(2)*erfinv(p); % Equal-tails. Pxx = Pxx/KMU; Pyy = Pyy/KMU; Pxy = Pxy/KMU; Pxxc = ff.*cPxx/KMU; Pxyc = ff.*cPxy/KMU; Pyyc = ff.*cPyy/KMU; P = [Pxx Pyy Pxy Txy Cxy Pxxc Pyyc Pxyc]; end P(:,1:3)=2*P(:,1:3)/Fs;%this corrects to get proper density - ngh may 2001 freq_vector = (select - 1)'*Fs/nfft; if nargout == 0, % do plots newplot; c = [max(Pxx-Pxxc,0) Pxx+Pxxc]; c = c.*(c>0); semilogy(freq_vector,Pxx,freq_vector,c(:,1),'--',... freq_vector,c(:,2),'--'); title('Pxx - X Power Spectral Density') xlabel('Frequency') if (isempty(y)), % single sequence case return end pause newplot; c = [max(Pyy-Pyyc,0) Pyy+Pyyc]; c = c.*(c>0); semilogy(freq_vector,Pyy,freq_vector,c(:,1),'--',... freq_vector,c(:,2),'--'); title('Pyy - Y Power Spectral Density') xlabel('Frequency') pause newplot; semilogy(freq_vector,abs(Txy)); title('Txy - Transfer function magnitude') xlabel('Frequency') pause newplot; plot(freq_vector,180/pi*angle(Txy)), ... title('Txy - Transfer function phase') xlabel('Frequency') pause newplot; plot(freq_vector,Cxy); title('Cxy - Coherence') xlabel('Frequency') elseif nargout ==1, Spec = P; elseif nargout ==2, Spec = P; f = freq_vector; elseif nargout ==3, Spec = P; f = freq_vector; Np=k; end function [msg,x,y,nfft,noverlap,window,Fs,p,dflag] = specchk(P) %SPECCHK Helper function for SPECTRUM % SPECCHK(P) takes the cell array P and uses each cell as % an input argument. Assumes P has between 1 and 7 elements. % Author(s): T. Krauss, 4-6-93 msg = []; if length(P{1})<=1 msg = 'Input data must be a vector, not a scalar.'; x = []; y = []; elseif (length(P)>1), if (all(size(P{1})==size(P{2})) & (length(P{1})>1) ) | ... length(P{2})>1, % 0ne signal or 2 present? % two signals, x and y, present x = P{1}; y = P{2}; % shift parameters one left P(1) = []; else % only one signal, x, present x = P{1}; y = []; end else % length(P) == 1 % only one signal, x, present x = P{1}; y = []; end % now x and y are defined; let's get the rest if length(P) == 1 nfft = min(length(x),256); window = hanning(nfft); noverlap = 0; Fs = 2; p = []; dflag = 'linear'; elseif length(P) == 2 if isempty(P{2}), dflag = 'linear'; nfft = min(length(x),256); elseif isstr(P{2}), dflag = P{2}; nfft = min(length(x),256); else dflag = 'linear'; nfft = P{2}; end window = hanning(nfft); noverlap = 0; Fs = 2; p = []; elseif length(P) == 3 if isempty(P{2}), nfft = min(length(x),256); else nfft=P{2}; end if isempty(P{3}), dflag = 'linear'; noverlap = 0; elseif isstr(P{3}), dflag = P{3}; noverlap = 0; else dflag = 'linear'; noverlap = P{3}; end window = hanning(nfft); Fs = 2; p = []; elseif length(P) == 4 if isempty(P{2}), nfft = min(length(x),256); else nfft=P{2}; end if isstr(P{4}) dflag = P{4}; window = hanning(nfft); else dflag = 'linear'; window = P{4}; window = window(:); % force window to be a column if length(window) == 1, window = hanning(window); end if isempty(window), window = hanning(nfft); end end if isempty(P{3}), noverlap = 0; else noverlap=P{3}; end Fs = 2; p = []; elseif length(P) == 5 if isempty(P{2}), nfft = min(length(x),256); else nfft=P{2}; end window = P{4}; window = window(:); % force window to be a column if length(window) == 1, window = hanning(window); end if isempty(window), window = hanning(nfft); end if isempty(P{3}), noverlap = 0; else noverlap=P{3}; end if isstr(P{5}) dflag = P{5}; Fs = 2; else dflag = 'linear'; if isempty(P{5}), Fs = 2; else Fs = P{5}; end end p = []; elseif length(P) == 6 if isempty(P{2}), nfft = min(length(x),256); else nfft=P{2}; end window = P{4}; window = window(:); % force window to be a column if length(window) == 1, window = hanning(window); end if isempty(window), window = hanning(nfft); end if isempty(P{3}), noverlap = 0; else noverlap=P{3}; end if isempty(P{5}), Fs = 2; else Fs = P{5}; end if isstr(P{6}) dflag = P{6}; p = []; else dflag = 'linear'; if isempty(P{6}), p = .95; else p = P{6}; end end elseif length(P) == 7 if isempty(P{2}), nfft = min(length(x),256); else nfft=P{2}; end window = P{4}; window = window(:); % force window to be a column if length(window) == 1, window = hanning(window); end if isempty(window), window = hanning(nfft); end if isempty(P{3}), noverlap = 0; else noverlap=P{3}; end if isempty(P{5}), Fs = 2; else Fs = P{5}; end if isempty(P{6}), p = .95; else p = P{6}; end if isstr(P{7}) dflag = P{7}; else msg = 'DFLAG parameter must be a string.'; return end end % NOW do error checking if (nfft= length(window)), msg = 'Requires NOVERLAP to be strictly less than the window length.'; end if (nfft ~= abs(round(nfft)))|(noverlap ~= abs(round(noverlap))), msg = 'Requires positive integer values for NFFT and NOVERLAP.'; end if ~isempty(p), if (prod(size(p))>1)|(p(1,1)>1)|(p(1,1)<0), msg = 'Requires confidence parameter to be a scalar between 0 and 1.'; end end if min(size(x))~=1, msg = 'Requires vector (either row or column) input.'; end if (min(size(y))~=1)&(~isempty(y)), msg = 'Requires vector (either row or column) input.'; end if (length(x)~=length(y))&(~isempty(y)), msg = 'Requires X and Y be the same length.'; end toolbox/mspec/confcoh.m0000755000175000001440000000033610321232073014141 0ustar nghusersfunction c=confcoh(k,percent) % function c=confcoh(#pieces,percent) % returns 95% confidence level for coherence amplitude if(nargin>1) alpha=(100-percent)/100; else alpha=.05; end c=sqrt(1-alpha.^(1 ./(k-1))); toolbox/mspec/confphs.m0000755000175000001440000000107110321232073014157 0ustar nghusersfunction c=confphs(coh,k,conf) %c=confphs(coherence,#pieces,conf) returns 95% confidence limit for phase % if conf(in %) given returns that confidence interval if nargin==3 conf=conf/100; else conf=.95; end gamsq=coh.^2; mu=2*(k-1); %from Abromowitz and Stegun a=[.27061 2.30753];b=[.04481 .99229 1]; p=(1-conf)/2; t=sqrt(log(1./p^2)); y=t-polyval(a,t)/polyval(b,t); y2=y^2;y3=y*y2;y5=y3*y2;y7=y5*y2; g1=.25*(y3+y); g2=(5*y5+16*y3+3*y)/96; g3=(3*y7+19*y5+17*y3-15*y)/384; t=y+(g1+(g2+g3./mu)./mu)./mu; ci=(180/pi)*asin(t.*sqrt((1-gamsq)./(mu.*gamsq))); c=real(ci); toolbox/mspec/confspec.m0000755000175000001440000000067010321232073014323 0ustar nghusersfunction conf=confspec(npc,prob) % function conf=confspec(npc,prob) % gives confidence limits for #pieces at prob(e.g. 0.95) % where conf returns (upper, mid and lower) positions of confidence limits suitable for use % on a log-log plot %int=(1-prob)/2; for i=1:length(npc) % confchi=chi2inv([int,.5,1-int],2*npc(i));%routine from statistics toolbox confchi=chisq(2*npc(i),prob);%my routine, faster conf(i,:)=2*npc(i)./confchi; end