Matlab spectrogram code

function X = spectrogram(x,nfft,fs,window,noverlap,doplot,dbdown);

%SPECTROGRAM Calculate spectrogram from signal.

% B = SPECTROGRAM(A,NFFT,Fs,WINDOW,NOVERLAP) calculates the

% spectrogram for the signal in vector A.

%

% NFFT is the FFT size used for each frame of A. It should be a

% power of 2 for fastest computation of the spectrogram.

%

% Fs is the sampling frequency. Since all processing parameters are

% in units of samples, Fs does not effect the spectrogram itself,

% but it is used for axis scaling in the plot produced when

% SPECTROGRAM is called with no output argument (see below).

%

% WINDOW is the length M window function applied, IN ZERO-PHASE

% FORM, to each frame of A. M cannot exceed NFFT. For M<NFFT,

% NFFT-M zeros are inserted in the FFT buffer (for interpolated

% zero-phase processing). The window should be supplied in CAUSAL

% FORM.

%

% NOVERLAP is the number of samples the sections of A overlap, if

% nonnegative. If negative, -NOVERLAP is the "hop size", \ie, the

% number of samples to advance successive windows. (The overlap is

% the window length minus the hop size.) The hop size is called

% NHOP below. NOVERLAP must be less than M.

%

% Thus, SPECTROGRAM splits the signal into overlapping segments of

% length M, windows each segment with the length M WINDOW vector, in

% zero-phase form, and forms the columns of B with their

% zero-padded, length NFFT discrete Fourier transforms.

%

% With no output argument B, SPECTROGRAM plots the dB magnitude of

% the spectrogram in the current figure, using

% IMAGESC(T,F,20*log10(ABS(B))), AXIS XY, COLORMAP(JET) so the low

% frequency content of the first portion of the signal is displayed

% in the lower left corner of the axes.

%

% 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 returned as a complex

% matrix with NFFT rows and

% k = floor((NX-NOVERLAP)/(length(WINDOW)-NOVERLAP))

% = floor((NX-NOVERLAP)/NHOP)

% columns. When A is real, only the NFFT/2+1 rows are needed when

% NFFT even, and the first (NFFT+1)/2 rows are sufficient for

% inversion when NFFT is odd.

%

% If DOPLOT is nonzero, the spectrogram will be plotted in any case.

% DBDOWN (default=100) sets the clip level below the maximum in dB.

%

% See also: Matlab's SPECGRAM function.

% 02/04/02/jos: Created

% 02/12/04/jos: Added dbdown

if nargin<7, dbdown=100; end

if nargin<6, doplot=0; end

if nargin<5, noverlap=256; end

if nargin<4, window=Hamming(512); end

if nargin<3, fs=1; end

if nargin<2, nfft=2048; end

x=x(:);

nfft

M = length(window)

if (M<2) error(...

'spectrogram: Expect complete window, not just its length');

end;

Modd = mod(M,2); % 0 if M even, 1 if odd

Mo2 = (M-Modd)/2;

w = window(:); % Make sure its a column

zp = zeros(nfft-M,1);

wzp = [w(Mo2+1:M);zp;w(1:Mo2)];

if noverlap<0

nhop = - noverlap;

noverlap = M-nhop;

else

nhop = M-noverlap;

end

nx = length(x);

nframes = 1+floor((nx-noverlap)/nhop);

X = zeros(nfft,nframes);

x = x(:); % make sure it's a column

xoff = 0;

for m=1:nframes-1

xframe = x(xoff+1:xoff+M); % extract frame of input data

xoff = xoff + nhop; % advance in-pointer by hop size

xzp = [xframe(Mo2+1:M);zp;xframe(1:Mo2)];

xw = wzp .* xzp;

X(:,m) = fft(xw);

end

if (nargout==0) | doplot

t = (0:nframes-1)*nhop/fs;

f = 0.001*(0:nfft-1)*fs/nfft;

Xdb = 20*log10(abs(X));

Xmax = max(max(Xdb));

% Clip lower limit to -dbdown dB so nulls don't dominate:

clipvals = [Xmax-dbdown,Xmax];

imagesc(t,f,Xdb,clipvals);

% grid;

axis xy;

colormap(jet);

xlabel('Time (sec)');

ylabel('Freq (kHz)');

end