function varargout=bayram(X,O,nfft,sint,R,L,X2) % BAYRAM(X,O,nfft,sint,R,L) % [tim,freq,SX,H,V]=BAYRAM(X,O,nfft,sint,R,L) % [tim,freq,SX,H,V,SX2,SX12,C2]=BAYRAM(X1,O,nfft,sint,R,L,X2) % % This program takes a time series and calculates the Hermite-multiple % window spectral estimate of it (including coherence for two signals) % % X analyzed signal (or X1 and X2) % O overlap of the windows in percent % nfft frequency sampling for each of the ffts of the windows % sint time sampling interval DX (default is 1) % R radius of conentration in the phase plane % L time resolution of the analysis (window length in samples) % X2 second signal in case coherence is required % % tim time domain of the signal % freq frequency domain % SX spectral estimate (or SX and SY and C2) % H Hermite matrix used % V matrix with eigenvalues % % See also COHEN % % EXAMPLE 1 % % nfft=1024; t=linspace(0,1,nfft); sint=indeks(diff(t),1); % X=chirp(t,0,1,200,'quadratic'); % X=sigmerge(X,rand(size(X)),20); % subplot(211); plot(t,X) % subplot(212); % bayram(X,99,256,sint,2); % % EXAMPLE 2 % % X=rand(2048,1); Y=rand(2048,1); sint=1; t=1:2048; R=3; L=128; % [tim,freq,SX,H,V,SX2,SX12,C2]=bayram(X,99,256,sint,R,L,Y); % imagesc(tim,freq,C2); axis xy; caxis([0 1]); % % EXAMPLE 3 % % subplot(211); [y,t,sint]=notch; subplot(212) % bayram(y,99,256,sint,2) % filename= 'lowry1'; % dapa= '/rosella/data14/fjsimons/GRAVITY/READYSETS'; % load(fullfile(dapa,filename)); % Bayram and Baraniuk, Proc. IEEE Sig. Proc. 173-176, 1996 % One or two signals? if exist('X2') & ~isempty(X2) flag=1; else; flag=0; end % Preliminaries - demean signals X=X(:)-mean(X); if flag; X2=X2(:)-mean(X2); end if flag ; if length(X)~=length(X2); error('Data vectors do not match'); end; end nx=length(X); if isempty(nfft); nfft=nx; end if isempty(sint); sint=1; end [H,V,K]=choice(R,L); disp([ 'Number of Hermite windows: ',num2str(K),'; Time window ',sprintf('%4.2f',L*sint)]) [fax,select]=fftaxis1D(X,nfft,sint*(nx-1)); O=floor(O/100*L); NWI=fix((nx-O)/(L-O)); % Number of windows disp([ 'Number of time windows: ',num2str(NWI)]) colindex = 1 + (0:(NWI-1))*(L-O); tim = (colindex-1)'*sint+(L*sint/2); % Middle time of window rowindex = (1:L)'; if length(X)<(L+colindex(NWI)-1) keyboard X(L+colindex(NWI)-1) = 0; % zero-pad last section of X if flag ; X2(L+colindex(NWI)-1) = 0; end end Ypre= zeros(L,NWI); if flag; Y2pre=Ypre; end % Put X into columns of Y with the proper offset Ypre= X(rowindex(:,ones(1,NWI))+colindex(ones(L,1),:)-1); if flag ; Y2pre= X2(rowindex(:,ones(1,NWI))+colindex(ones(L,1),:)-1); end % Apply the windows to the array of offset signal segments. SX=zeros(length(select),NWI); if flag; SX2=SX; SX12=SX;end for index=1:K WI=H(:,index); Y = WI(:,ones(1,NWI)).*Ypre; if flag; Y2 = WI(:,ones(1,NWI)).*Y2pre; end % Now FFT Y; this does the columns FY = fft(Y,nfft); if flag; FY2= fft(Y2,nfft); end if isreal(X) % Assume both signals are both REAL or COMPLEX FY = FY(select,:); if flag; FY2=FY2(select,:); end end % Make weighted sum SX=V(index)*(abs(FY).^2)+SX; if flag; SX2=V(index)*(abs(FY2).^2)+SX2; end if flag; SX12=V(index)*(FY.*conj(FY2))+SX12; end end SX=SX/sum(V); if flag SX2=SX2/sum(V); SX12=SX12/sum(V); C2=abs(SX12).^2./SX./SX2; end % Which property is plotted? prop=SX; %prop=decibel(SX); if ~nargout axes(gca) imagesc(tim,fax,prop); axis xy; colormap(jet) xlim([0 (nx-1)*sint]); xlabel('Time') ylabel('Frequency') title([ 'Number of Hermite windows: ',num2str(K),'; Time window ',sprintf('%4.2f',L*sint)]) else posout={ 'tim','fax','SX','H','V','SX2','SX12','C2'}; for index=1:nargout eval([ 'varargout{index}=',posout{index},';']) end end