annotate Code/Descriptors/Matlab/MPEG7/FromWeb/h_specgram2.m @ 4:92ca03a8fa99 tip

Update to ICASSP 2013 benchmark

function [yo,fo,to] = h_specgram2(varargin)
%SPECGRAM Calculate spectrogram from signal.
%   B = SPECGRAM(A,NFFT,Fs,WINDOW,SHIFT) calculates the spectrogram for 
%   the signal in vector A.  SPECGRAM splits the signal into overlapping 
%   segments, windows each with the WINDOW vector and forms the columns of
%   B with their zero-padded, length NFFT discrete Fourier transforms.  Thus
%   each column of B contains an estimate of the short-term, time-localized
%   frequency content of the signal A.  Time increases linearly across the 
%   columns of B, from left to right.  Frequency increases linearly down 
%   the rows, starting at 0.  If A is a length NX complex signal, B is a 
%   complex matrix with NFFT rows and 
%        k = fix((NX-NOVERLAP)/(length(WINDOW)-NOVERLAP)) 
%   columns, where NOVERLAP = length(WINDOW)-mean(SHIFT)
%   If A is real, B still has k columns but the higher frequency
%   components are truncated (because they are redundant); in that case,
%   SPECGRAM returns B with NFFT/2+1 rows for NFFT even and (NFFT+1)/2 rows 
%   for NFFT odd.  If you specify a scalar for WINDOW, SPECGRAM uses a 
%   Hanning window of that length.  WINDOW must have length smaller than
%   or equal to NFFT and greater than NOVERLAP.  NOVERLAP is the number of
%   samples the sections of A overlap.  Fs is the sampling frequency
%   which does not effect the spectrogram but is used for scaling plots.
%
%   [B,F,T] = SPECGRAM(A,NFFT,Fs,WINDOW,NOVERLAP) returns a column of 
%   frequencies F and one of times T at which the spectrogram is computed.
%   F has length equal to the number of rows of B, T has length k. If you 
%   leave Fs unspecified, SPECGRAM assumes a default of 2 Hz.
%
%   B = SPECGRAM(A) produces the spectrogram of the signal A using default
%   settings; the defaults are NFFT = minimum of 256 and the length of A, a
%   Hanning window of length NFFT, and NOVERLAP = length(WINDOW)/2.  You 
%   can tell SPECGRAM to use the default for any parameter by leaving it 
%   off or using [] for that parameter, e.g. SPECGRAM(A,[],1000)
%
%   See also PWELCH, CSD, COHERE and TFE.

%   Author(s): L. Shure, 1-1-91
%              T. Krauss, 4-2-93, updated
%   Copyright 1988-2000 The MathWorks, Inc.
%   $Revision: 1.6 $  $Date: 2000/06/09 22:07:35 $

error(nargchk(1,5,nargin))
[msg,x,nfft,Fs,window,shift]=specgramchk(varargin);
error(msg)
    
nx = length(x);
nwind = length(window);
if nx < nwind    % zero-pad x if it has length less than the window length
    x(nwind)=0;  nx=nwind;
end
x = x(:); % make a column vector for ease later
window = window(:); % be consistent with data set
noverlap = nwind - mean(shift);
ncol = fix((nx-noverlap)/(nwind-noverlap));

totalshift = sum(shift);
numshift = length(shift);
progshift = (cumsum(shift)-shift(1))';
overshift = 1:totalshift:nx;
shifts = progshift(:,ones(1,length(overshift)))+overshift(ones(length(progshift),1),:);
colindex = reshape(shifts,1,prod(size(shifts)));
colindex = colindex(1:ncol);
%colindex = 1 + (0:(ncol-1))*(nwind-noverlap);
rowindex = (1:nwind)';
if length(x)<(nwind+colindex(ncol)-1)
    x(nwind+colindex(ncol)-1) = 0;   % zero-pad x
end

y = zeros(nwind,ncol);

% put x into columns of y with the proper offset
% should be able to do this with fancy indexing!
y(:) = x(rowindex(:,ones(1,ncol))+colindex(ones(nwind,1),:)-1);

% Apply the window to the array of offset signal segments.
y = window(:,ones(1,ncol)).*y;

% now fft y which does the columns
y = fft(y,nfft);
if ~any(any(imag(x)))    % x purely real
    if rem(nfft,2),    % nfft odd
        select = [1:(nfft+1)/2];
    else
        select = [1:nfft/2+1];
    end
    y = y(select,:);
else
    select = 1:nfft;
end
f = (select - 1)'*Fs/nfft;

t = (colindex-1)'/Fs;

% take abs, and use image to display results
if nargout == 0
%     newplot;
%     if length(t)==1
%         imagesc([0 1/f(2)],f,20*log10(abs(y)+eps));axis xy; colormap(jet)
%     else
%         imagesc(t,f,20*log10(abs(y)+eps));axis xy; colormap(jet)
%     end
%     xlabel('Time')
%     ylabel('Frequency')
elseif nargout == 1,
    yo = y;
elseif nargout == 2,
    yo = y;
    fo = f;
elseif nargout == 3,
    yo = y;
    fo = f;
    to = t;
end

function [msg,x,nfft,Fs,window,shift] = specgramchk(P)
%SPECGRAMCHK Helper function for SPECGRAM.
%   SPECGRAMCHK(P) takes the cell array P and uses each cell as 
%   an input argument.  Assumes P has between 1 and 5 elements.

msg = [];

x = P{1}; 
if (length(P) > 1) & ~isempty(P{2})
    nfft = P{2};
else
    nfft = min(length(x),256);
end
if (length(P) > 2) & ~isempty(P{3})
    Fs = P{3};
else
    Fs = 2;
end
if length(P) > 3 & ~isempty(P{4})
    window = P{4}; 
else
    if length(nfft) == 1
        window = hanning(nfft);
    else
        msg = 'You must specify a window function.';
    end
end
if length(window) == 1, window = hanning(window); end
if (length(P) > 4) & ~isempty(P{5})
    shift = P{5};
else
    shift = ceil(length(window)/2);
end

% NOW do error checking
if (length(nfft)==1) & (nfft<length(window)), 
    msg = 'Requires window''s length to be no greater than the FFT length.';
end
if (min(shift) <= 0),
    msg = 'Requires SHIFT to be strictly positive.';
end
if (length(nfft)==1) & (nfft ~= abs(round(nfft)))
    msg = 'Requires positive integer values for NFFT.';
end
if sum(shift ~= abs(round(shift))),
    msg = 'Requires positive integer value for NOVERLAP.';
end
if min(size(x))~=1,
    msg = 'Requires vector (either row or column) input.';
end