STFT and ISTFT
Wednesday, September 17th, 2008 (6:38 pm), under MATLAB, cocktail party, dsp, utilities.Problem: I need to invert a spectrogram to listen to the signal. And I want to do it the best way possible.
Solution: It isn’t all that sophisticated. Here’s a limited implementation of Yang’s ICASSP 2008 solution; which is — surprise surprise — quite nearly the same as the ubiquitous heuristic overlap and add istft.
It’s limited in the following ways:
- The least-squares solution is not implemented. This means that LS solution is available only for regularly spaced frequency bands since in that case LS solution reduces to p=2 solution.
- MATLAB fft/ifft routines are used; and they take care of the P matrix (excess zeros to round to powers of 2).
- The matrix algebra is not exploited — this is a quick and dirty test run.
stft.m
function [sgram, freqs, tt]=stft(x, fs, pp)
% get stft of signal
if(size(x, 1)~=1 && size(x, 2)~=1), error('1D vector expected'); end;
if(size(x, 1)<size(x, 2)), x=x'; end;
hoplen=ceil(fs*pp.hoplen);
winlen=ceil(fs*pp.winlen);
if(winlen<hoplen)
winlen=hoplen;
error('winlen less than hoplen; parts of signal are naked!');
end;
if(rem(winlen, 2)==0), winlen=winlen+1; end;
if(length(x)<N) x(end+1:N)=0.0; end;
win=hamming(winlen);
nhops=floor(length(x)/hoplen);
N=(nhops-1)*hoplen+winlen;
if(length(x)<N) x(end+1:N)=0.0; end;
wby2=ceil(winlen/2);
sgram=zeros(wby2, nhops);
for ii=1:nhops,
%[ii length(win) (ii-1)*nhops (ii-1)*hoplen+winlen N]
tmpft=fft(win.*x((ii-1)*hoplen+[1:winlen]));
sgram(:, ii)=tmpft(1:wby2);
end;
freqs=([1:wby2]/wby2)*fs*0.5;
tt=([1:nhops]-1)*pp.hoplen;
istft.m
function [sig, winseq]=istft(sgram, freqs, fs, pp)
% implement yang 2008 ICASSP
% prepare parameters
hoplen=ceil(fs*pp.hoplen);
winlen=ceil(fs*pp.winlen);
if(winlen<hoplen)
winlen=hoplen;
error('winlen less than hoplen; parts of signal are naked!');
end;
%winlen=2*wby2;
if(rem(winlen, 2)==0), winlen=winlen+1; end;
winlen=2*size(sgram, 1)-1;
%wby2=ceil(winlen/2);
win=hamming(winlen);
%nhops=floor(length(x)/hoplen);
nhops=size(sgram, 2);
N=(nhops-1)*hoplen+winlen;
winpm1=win.^(pp.p-1);
winp=win.^(pp.p);
winseq=zeros(N, 1);
sig=zeros(N, 1);
for ii=1:nhops,
winseq((ii-1)*hoplen+[1:winlen])=winp+...
winseq((ii-1)*hoplen+[1:winlen]);
end;
for ii=1:nhops,
tmpsig=(ifft([sgram(1:end, ii); conj(sgram(end-1:-1:1, ii))], 'symmetric'));
tmpsig=tmpsig.*winpm1;
sig((ii-1)*hoplen+[1:winlen])=sig((ii-1)*hoplen+[1:winlen])+...
tmpsig;
end;
sig=sig./winseq;
test_istft.m
fs=1000;
x=chirp(t, 0, 1, 150);
gg=specgram(x);
pp.hoplen=128/fs; pp.winlen=2*128/fs;
[sgram, ff, tt]=stft(x, fs, pp);
figure_docked;
subplot(2, 1, 1); surfo(10*log10(abs(gg)));
subplot(2, 1, 2); surfo(10*log10(abs(sgram)));
pp.p=2;
[xn, ws]=istft(sgram, ff, fs, pp);
figure_docked;
col={'r'; 'g'; 'b'; 'c'; 'y'; 'm'; 'k'};
subplot(2, 1, 1); plot(x, col{1}); hold on;
for ppp=3:-1:0,
pp.p=ppp;
[xn, ws]=istft(sgram, ff, fs, pp);
subplot(2, 1, 1); hold on; plot(xn, col{ppp+2});
subplot(2, 1, 2); hold on; plot(ws, col{ppp+2});
end;
sdr.m
function d=sdr(xx, xn)
len=min(length(xx), length(xn));
if(size(xx, 1)<size(xx, 2)) xx=xx'; end;
if(size(xn, 1)<size(xn, 2)) xn=xn'; end;
[size(xx) size(xn)]
d=10*log10(sum((xx(1:len).^2))./sum((xx(1:len)-xn(1:len)).^2));