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
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
% 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');
Matlab for Unwrapping Spectral Phase
Matlab for Finding Interpolated Spectral Peaks