Sign in

username:

password:



Not a member?

Search Online Books



Search tips

Free Online Books



Chapters

See Also

Embedded SystemsFPGAElectronics
Chapter Contents:

Search Spectral Audio Signal Processing

  

Book Index | Global Index


Would you like to be notified by email when Julius Orion Smith III publishes a new entry into his blog?

  

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


Previous: Matlab listing: invmyspectrogram.m
Next: Matlab for Unwrapping Spectral Phase

Order a Hardcopy of Spectral Audio Signal Processing


About the Author: Julius Orion Smith III
Julius Smith's background is in electrical engineering (BS Rice 1975, PhD Stanford 1983). He is presently Professor of Music and Associate Professor (by courtesy) of Electrical Engineering at Stanford's Center for Computer Research in Music and Acoustics (CCRMA), teaching courses and pursuing research related to signal processing applied to music and audio systems. See http://ccrma.stanford.edu/~jos/ for details.


Comments


No comments yet for this page


Add a Comment
You need to login before you can post a comment (best way to prevent spam). ( Not a member? )