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