% Usage: spectro(f, T, dT, overlap, 'windowtype') % Given a vector "f" of length N representing a signal of analog time % length "T" this computes a spectrogram for f using windows of length % dT with overlap percentage given by "overlap". function spectro(f, T, dT, overlap, window) % Convert f to a row vector. [n,N] = size(f); if n==1 f2 = f; else f2 = f'; end [n,N] = size(f2); %Set up grayscale color map. mappy = [(0:255)/256; (0:255)/256; (0:255)/256]'; colormap(mappy); %Bookkeeping m = floor(dT/T*N); %Window size r = floor((100-overlap)*m/100); %Index stepsize through signal H = 0.5*N/T; %Highest frequency dt = r/N*T; %Step in terms of time M = floor((N-m)/r); %Number of blocks we'll do G = zeros(m,M+1); %Array to hold transforms for j=0:M if strcmp(window,'rect') g = f2(r*j+1:r*j+m).*rectwindow(m,1,m); elseif strcmp(window,'tri') g = f2(r*j+1:r*j+m).*triwindow(m,1,m); elseif strcmp(window,'gauss') g = f2(r*j+1:r*j+m).*gausswindow(m,1,m,5.0); elseif strcmp(window,'hamming') g = f2(r*j+1:r*j+m).*hammingwindow(m,1,m); else %default to rectangular if j==0 display('Window not recognized---default to rectangular'); end; g = f2(r*j+1:r*j+m).*rectwindow(m,1,m); end; Gg = fft(g); G(:,j+1) = Gg'; end; ma = max(max(log(1+abs(G).^2))); X = ([0:M]+0.5)*dt; Y = [H:-1:0]; image(X,Y,255*log(1+abs(G(floor(m/2):-1:1,:).^2))/ma); axis xy; ylabel('Hertz'); xlabel('seconds');