## 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');

**Next Section:**

Matlab for Unwrapping Spectral Phase

**Previous Section:**

Matlab for Finding Interpolated Spectral Peaks