katkost@1: function [y,yh,ys] = spsmodel(x,fs,w,N,t,maxnS,stocf) katkost@1: %=> analysis/synthesis of a sound using the sinusoidal harmonic model katkost@1: % x: input sound, fs: sampling rate, w: analysis window (odd size), katkost@1: % N: FFT size (minimum 512), t: threshold in negative dB, katkost@1: % maxnS: maximum number of sinusoids, katkost@1: % stocf: decimation factor of mag spectrum for stochastic analysis katkost@1: % y: output sound, yh: harmonic component, ys: stochastic component katkost@1: M = length(w); % analysis window size katkost@1: Ns = 1024; % FFT size for synthesis katkost@1: H = 256; % hop size for analysis and synthesis katkost@1: N2 = N/2+1; % half-size of spectrum katkost@1: soundlength = length(x); % length of input sound array katkost@1: hNs = Ns/2; % half synthesis window size katkost@1: hM = (M-1)/2; % half analysis window size katkost@1: pin = max(N2,1+hM) % initialize sound pointer to middle of analysis window katkost@1: pend = soundlength-max(N2,hM); % last sample to start a frame katkost@1: fftbuffer = zeros(N,1); % initialize buffer for FFT katkost@1: yh = zeros(soundlength+Ns/2,1); % output sine component katkost@1: ys = zeros(soundlength+Ns/2,1); % output residual component katkost@1: w = w/sum(w); % normalize analysis window katkost@1: sw = zeros(Ns,1); katkost@1: ow = triang(2*H-1); % overlapping window katkost@1: ovidx = Ns/2+1-H+1:Ns/2+H; % overlap indexes katkost@1: sw(ovidx) = ow(1:2*H-1); katkost@1: bh = blackmanharris(Ns); % synthesis window katkost@1: bh = bh ./ sum(bh); % normalize synthesis window katkost@1: wr = bh; % window for residual katkost@1: sw(ovidx) = sw(ovidx) ./ bh(ovidx); katkost@1: sws = H*hanning(Ns)/2; % synthesis window for stochastic katkost@1: lastysloc = zeros(maxnS,1); % initialize synthesis harmonic locations katkost@1: ysphase = 2*pi*rand(maxnS,1); % initialize synthesis harmonic phases katkost@1: fridx = 0; katkost@1: katkost@1: katkost@1: katkost@1: katkost@1: katkost@1: while pint) .* (mX(2:N2-1)>mX(3:N2)) ... katkost@1: .* (mX(2:N2-1)>mX(1:N2-2))); % find peaks katkost@1: [ploc,pmag,pphase] = peakinterp(mX,pX,ploc); % refine peak values katkost@1: % sort by magnitude katkost@1: [smag,I] = sort(pmag(:),1,'descend'); katkost@1: nS = min(maxnS,length(find(smag>t))); katkost@1: sloc = ploc(I(1:nS)); katkost@1: sphase = pphase(I(1:nS)); katkost@1: if (fridx==0) katkost@1: % update last frame data for first frame katkost@1: lastnS = nS; katkost@1: lastsloc = sloc; katkost@1: lastsmag = smag; katkost@1: lastsphase = sphase; katkost@1: end katkost@1: % connect sinusoids to last frame lnS (lastsloc,lastsphase,lastsmag) katkost@1: sloc(1:nS) = (sloc(1:nS)~=0).*((sloc(1:nS)-1)*Ns/N); % synth. locs katkost@1: lastidx = zeros(1,nS); katkost@1: for i=1:nS katkost@1: [dev,idx] = min(abs(sloc(i) - lastsloc(1:lastnS))); katkost@1: lastidx(i) = idx; katkost@1: end katkost@1: ri= pin-hNs; % input sound pointer for residual analysis katkost@1: xr = x(ri:ri+Ns-1).*wr(1:Ns); % window the input sound katkost@1: Xr = fft(fftshift(xr)); % compute FFT for residual analysis katkost@1: Xh = genspecsines(sloc,smag,sphase,Ns); % generate sines katkost@1: Xr = Xr-Xh; % get the residual complex spectrum katkost@1: mXr = 20*log10(abs(Xr(1:Ns/2+1))); % magnitude spectrum of residual katkost@1: mXsenv = decimate(max(-200,mXr),stocf); % decimate the magnitude spectrum katkost@1: % and avoid -Inf katkost@1: %-----synthesis data-----% katkost@1: ysloc = sloc; % synthesis locations katkost@1: ysmag = smag(1:nS); % synthesis magnitudes katkost@1: mYsenv = mXsenv; % synthesis residual envelope katkost@1: katkost@1: katkost@1: %-----transformations-----% katkost@1: katkost@1: %-----time mapping-----% katkost@1: outsoundlength = soundlength*2; katkost@1: TM=[ 0 soundlength/fs; % input time (sec) katkost@1: 0 outsoundlength/fs ]; % output time (sec) katkost@1: katkost@1: %%-----frequency stretch-----% katkost@1: %fstretch = 1.1; katkost@1: %ysloc = ysloc .* (fstretch.^[0:length(ysloc)-1]'); katkost@1: katkost@1: %-----synthesis-----% katkost@1: if (fridx==0) katkost@1: lastysphase = ysphase; katkost@1: end katkost@1: if (nS>lastnS) katkost@1: lastysphase = [ lastysphase ; zeros(nS-lastnS,1) ]; katkost@1: lastysloc = [ lastysloc ; zeros(nS-lastnS,1) ]; katkost@1: end katkost@1: ysphase = lastysphase(lastidx(1:nS)) + 2*pi*(lastysloc(lastidx(1:nS))+ysloc)/2/Ns*H; % propagate phases katkost@1: lastysloc = ysloc; katkost@1: lastysphase = ysphase; katkost@1: lastnS = nS; % update last frame data katkost@1: lastsloc = sloc; % update last frame data katkost@1: lastsmag = smag; % update last frame data katkost@1: lastsphase = sphase; % update last frame data katkost@1: Yh = genspecsines(ysloc,ysmag,ysphase,Ns); % generate sines katkost@1: mYs = interp(mYsenv,stocf); % interpolate to original size katkost@1: roffset = ceil(stocf/2)-1; % interpolated array offset katkost@1: mYs = [ mYs(1)*ones(roffset,1); mYs(1:Ns/2+1-roffset) ]; katkost@1: mYs = 10.^(mYs/20); % dB to linear magnitude katkost@1: pYs = 2*pi*rand(Ns/2+1,1); % generate phase spectrum with random values katkost@1: mYs1 = [mYs(1:Ns/2+1); mYs(Ns/2:-1:2)]; % create complete magnitude spectrum katkost@1: pYs1 = [pYs(1:Ns/2+1); -1*pYs(Ns/2:-1:2)]; % create complete phase spectrum katkost@1: Ys = mYs1.*cos(pYs1)+1i*mYs1.*sin(pYs1); % compute complex spectrum katkost@1: yhw = fftshift(real(ifft(Yh))); % sines in time domain using inverse FFT katkost@1: ysw = fftshift(real(ifft(Ys))); % stochastic in time domain using IFFT katkost@1: yh(ri:ri+Ns-1) = yh(ri:ri+Ns-1)+yhw(1:Ns).*sw; % overlap-add for sines katkost@1: ys(ri:ri+Ns-1) = ys(ri:ri+Ns-1)+ysw(1:Ns).*sws; % overlap-add for stochastic katkost@1: pin = pin+H; % advance the sound pointer katkost@1: fridx = fridx+1; katkost@1: end katkost@1: y= yh+ys; % sum sines and stochastic