DSPRelated.com
Free Books

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 listing: unwrap.m
Previous Section:
Matlab listing: invmyspectrogram.m