Matlab for Computing Spectrograms

The Matlab Signal Processing Toolbox provides the command spectrogram for computing and displaying a spectrogram (and Octave has the command stft). As a side effect, it returns the complex STFT data in a matrix.

The myspectrogram function below illustrates computation of a spectrogram in matlab for purposes of basic spectrum analysis. It is compatible with the basic functionality of Matlab's spectrogram function, but does not implement all of its features.

We also include invmyspectrogram, which inverts a spectrogram to reproduce the corresponding time-domain signal. Finally, testmyspectrogram is included to illustrate their performance on test signal.

Matlab listing: myspectrogram.m

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

%MYSPECTROGRAM Calculate spectrogram from signal.
% B = MYSPECTROGRAM(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
% MYSPECTROGRAM 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", i.e., 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.
%
% If doplot is nonzero, or if there is no output argument, the
% spectrogram is displayed.
%
% When the spectrogram is displayed, it is ``clipped'' dbdown dB
% below its maximum magnitude.  The default clipping level is 100
% dB down.
%
% Thus, MYSPECTROGRAM 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, MYSPECTROGRAM 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.
%
% See also: Matlab's SPECTROGRAM and Octave's STFT function.
% 02/04/02/jos: Created
% 02/12/04/jos: Added dbdown
% 07/23/08/jos: Changed name from SPECTROGRAM to MYSPECTROGRAM

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(:); % make sure it's a column

M = length(window);
if (M<2) error(...
  'myspectrogram: Expect complete window, not just its length');
end;
if (M<2) error(...
  'myspectrogram: Expect complete window, not just its length');
end;
if length(x)<M % zero-pad to fill a window:
  x = [x;zeros(M-length(x),1)];
end;
Modd = mod(M,2); % 0 if M even, 1 if odd
Mo2 = (M-Modd)/2;
w = window(:); % Make sure it's a column

if noverlap<0
  nhop = - noverlap;
  noverlap = M-nhop;
else
  nhop = M-noverlap;
end

nx = length(x);
%nframes = 1+floor((nx-noverlap)/nhop);
nframes = 1+ceil(nx/nhop);

X = zeros(nfft,nframes); % allocate output spectrogram

zp = zeros(nfft-M,1); % zero-padding for each FFT
xframe = zeros(M,1);
xoff = 0 - Mo2; % input time offset = half a frame
for m=1:nframes
%  M,Mo2,xoff,nhop
  if xoff<0
    xframe(1:xoff+M) = x(1:xoff+M); % partial input data frame
  else
    if xoff+M > nx
      xframe = [x(xoff+1:nx);zeros(xoff+M-nx,1)];
    else
      xframe = x(xoff+1:xoff+M); % input data frame
    end
  end
  xw = w .* xframe; % Apply window
  xwzp = [xw(Mo2+1:M);zp;xw(1:Mo2)];
  X(:,m) = fft(xwzp);
  xoff = xoff + nhop; % advance input offset by hop size
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

Matlab listing: invmyspectrogram.m

function a = invmyspectrogram(b,hop)
%INVMYSPECTROGRAM Resynthesize a signal from its spectrogram.
%   A = INVMYSPECTROGRAM(B,NHOP)
%   B = complex array of STFT values as generated by MYSPECTROGRAM.
%   The number of rows of B is taken to be the FFT size, NFFT.
%   INVMYSPECTROGRAM resynthesizes A by inverting each frame of the
%   FFT in B, and overlap-adding them to the output array A.
%   NHOP is the overlap-add offset between successive IFFT frames.
%
%   See also: MYSPECTROGRAM

[nfft,nframes] = size(b);

No2 = nfft/2; % nfft assumed even
a = zeros(1, nfft+(nframes-1)*hop);
xoff = 0 - No2; % output time offset = half of FFT size
for col = 1:nframes
  fftframe = b(:,col);
  xzp = ifft(fftframe);
  % xzp = real(xzp); % if signal known to be real
  x = [xzp(nfft-No2+1:nfft); xzp(1:No2)];
  if xoff<0 % FFT's "negative-time indices" are out of range
    ix = 1:xoff+nfft;
    a(ix) = a(ix) + x(1-xoff:nfft)'; % partial frames out
  else
    ix = xoff+1:xoff+nfft;
    a(ix) = a(ix) + x';  % overlap-add reconstruction
  end
  xoff = xoff + hop;
end

Matlab listing: testmyspectrogram.m

% testmyspectrogram.m (tested July 2008 in Octave 3.0)

fs=22050;	% sampling frequency
D=1;		% test signal duration in seconds
L = ceil(fs*D)+1; % signal duration in samples
n = 0:L-1;      % discrete-time axis (samples)
t = n/fs;       % discrete-time axis (sec)
%c = ones(1,L); % dc test (good first COLA check)
c = chirp(t,0,D,fs/2); % sine sweep from 0 Hz to fs/2 Hz
windur = 0.01;  % window length in sec
M = 2*round((windur*fs-1)/2); % win length in samples (even)
%M = 2*round((windur*fs-1)/2)+1; % win length in samples (odd)
Modd = mod(M,2);
hop = (M-Modd)/2;

if Modd
  w = 0.54 - 0.46 * cos(2*pi*(0:M-1)'/(M-1)); % causal Hamming window
  w(1)=w(1)/2;  w(M)=w(M)/2;  % modify for constant-overlap-add
else
  w = 0.54 - 0.46 * cos(2*pi*(0:M-1)'/M); % causal Hamming window
end
w = w/(2*0.54); % scale for unity overlap-add
nfft = 2^(3+nextpow2(M));

% For zero-phase processing, the signal should have at least
% half a window of leading and trailing zeros
zp = zeros(1,(M+Modd)/2);
x = [zp,c,zp];

figure(1); clf;

X = myspectrogram(x,nfft,fs,w,-hop,1);
%or myspectrogram(x,nfft,fs,w,M-hop,1);

title('Spectrogram of test signal');
% print -depsc 'testmyspectrogram.eps'; % write plot to disk

%xh = invmyspectrogram(X,hop,8);
xh = invmyspectrogram(X,hop);
xh = real(xh); % imaginary part is round-off noise since x was real
xmxh = x - xh(1:length(x)); % Have extra zeros at end of xh
err=norm(xmxh)/norm(x);
disp(sprintf('L2 norm of relative reconstruction error  = %g',err));
% I observe 8E-16 for both odd and even win length cases - jos

figure(2); clf;
%n1 = round(L/8); n2 = n1+100;
n1 = 1; n2 = n1+3000;
pn = n1:n2; % plot indices
%plot(n(pn),[x(pn)',xh(pn)']);
%legend('original','resynthesized');
%title('Original and resynthesized test signal');
plot(n(pn),x(pn)-xh(pn));
title('Original minus resynthesized test signal');

figure(3); clf;
Xh=myspectrogram(xh,nfft,fs,w,-hop,1);
title('Spectrogram of resynthesized test signal');

Sign in