# HG changeset patch
# User Katerina <katkost@gmail.com>
# Date 1366457751 -3600
# Node ID 881c3acf1164262c4f8de82dd2449e46f5a206f2
# Parent  59c6861155d9aa48b4f956c60f8ed981b7656467
matlab code Kat

diff -r 59c6861155d9 -r 881c3acf1164 code Kat/hpsmodel.m
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/code Kat/hpsmodel.m	Sat Apr 20 12:35:51 2013 +0100
@@ -0,0 +1,101 @@
+function [y,yh,ys,fr0] = hpsmodel(x,fs,w,N,t,nH,minf0,maxf0,f0et,maxhd,stocf) 
+%=> analysis/synthesis of a sound using the sinusoidal harmonic model 
+% x: input sound, fs: sampling rate, w: analysis window (odd size),  
+% N: FFT size (minimum 512), t: threshold in negative dB,  
+% nH: maximum number of harmonics, minf0: minimum f0 frequency in Hz,  
+% maxf0: maximim f0 frequency in Hz,  
+% f0et: error threshold in the f0 detection (ex: 5), 
+% maxhd: max. relative deviation in harmonic detection (ex: .2) 
+% stocf: decimation factor of mag spectrum for stochastic analysis 
+% y: output sound, yh: harmonic component, ys: stochastic component 
+M = length(w);   % analysis window size 
+Ns = 1024;                               % FFT size for synthesis 
+H = 256;                                 % hop size for analysis and synthesis 
+N2 = N/2+1;                              % half-size of spectrum 
+soundlength = length(x);                 % length of input sound array 
+hNs = Ns/2;                              % half synthesis window size 
+hM = (M-1)/2;                            % half analysis window size 
+pin = max(hNs+1,1+hM);   % initialize sound pointer to middle of analysis window 
+pend = soundlength-max(hM,hNs);            % last sample to start a frame 
+fftbuffer = zeros(N,1);                  % initialize buffer for FFT 
+yh = zeros(soundlength+Ns/2,1);          % output sine component 
+ys = zeros(soundlength+Ns/2,1);          % output residual component 
+w = w/sum(w);                            % normalize analysis window 
+sw = zeros(Ns,1); 
+ow = triang(2*H-1);                      % overlapping window 
+ovidx = Ns/2+1-H+1:Ns/2+H;               % overlap indexes 
+sw(ovidx) = ow(1:2*H-1); 
+bh = blackmanharris(Ns);                 % synthesis window 
+bh = bh ./ sum(bh);                      % normalize synthesis window 
+wr = bh;                                 % window for residual  
+sw(ovidx) = sw(ovidx) ./ bh(ovidx); 
+sws = H*hanning(Ns);                     % synthesis window for stochastic 
+
+i = 0;
+while pin<pend 
+  i = i+1;
+  %-----analysis-----% 
+  xw = x(pin-hM:pin+hM).*w(1:M);         % window the input sound 
+  fftbuffer(1:(M+1)/2) = xw((M+1)/2:M);  % zero-phase window in fftbuffer 
+  fftbuffer(N-(M-1)/2+1:N) = xw(1:(M-1)/2); 
+  X = fft(fftbuffer);                    % compute the FFT 
+  mX = 20*log10(abs(X(1:N2)));           % magnitude spectrum  
+  pX = unwrap(angle(X(1:N/2+1)));        % unwrapped phase spectrum  
+  ploc = 1 + find((mX(2:N2-1)>t) .* (mX(2:N2-1)>mX(3:N2)) ... 
+                  .* (mX(2:N2-1)>mX(1:N2-2)));          % find peaks 
+  [ploc,pmag,pphase] = peakinterp(mX,pX,ploc);          % refine peak values 
+  f0 = f0detection(mX,fs,ploc,pmag,f0et,minf0,maxf0);   % find f0 
+  fr0(i)=f0;
+  hloc = zeros(nH,1);                    % initialize harmonic locations 
+  hmag = zeros(nH,1)-100;                % initialize harmonic magnitudes 
+  hphase = zeros(nH,1);                  % initialize harmonic phases 
+  hf = (f0>0).*(f0.*(1:nH));             % initialize harmonic frequencies 
+  hi = 1;                                % initialize harmonic index 
+  npeaks = length(ploc);                 % number of peaks found 
+  while (f0>0 && hi<=nH && hf(hi)<fs/2)  % find harmonic peaks 
+    [dev,pei] = min(abs((ploc(1:npeaks)-1)/N*fs-hf(hi)));    % closest peak 
+    if ((hi==1 || ~any(hloc(1:hi-1)==ploc(pei))) && dev<maxhd*hf(hi)) 
+      hloc(hi) = ploc(pei);              % harmonic locations 
+      hmag(hi) = pmag(pei);              % harmonic magnitudes 
+      hphase(hi) = pphase(pei);          % harmonic phases 
+    end 
+    hi = hi+1;                           % increase harmonic index 
+  end 
+  hloc(1:hi-1) = (hloc(1:hi-1)~=0).*((hloc(1:hi-1)-1)*Ns/N+1);  % synth. locs 
+  ri= pin-hNs;                     % input sound pointer for residual analysis 
+  xr = x(ri:ri+Ns-1).*wr(1:Ns);          % window the input sound 
+  Xr = fft(fftshift(xr));                % compute FFT for residual analysis 
+  Yh = genspecsines(hloc(1:hi-1),hmag,hphase,Ns);             % generate sines 
+  Yr = Xr-Yh;                            % get the residual complex spectrum 
+  mYr = abs(Yr(1:Ns/2+1));               % magnitude spectrum of residual 
+  %mYs = stochenvelope(mYr,stocf);
+  %-----transformations-----% 
+  mYsenv = decimate(mYr,stocf,1); % decimate the magnitude spectrum
+  %-----synthesis-----%
+  mYs = interp(mYsenv,stocf,1); % interpolate to original size
+  
+% n=1:N/2+1; 
+% plot(n/N*Ns,mX); %plotting the original spectrum
+% hold on;
+% plot(20*log10(abs(mYs)), 'r'); %plotting the approximation done by the decimate function
+% 
+% hold on;
+% plot(hloc, hmag, 'g*');
+% hold off;
+% pause
+  
+  roffset = ceil(stocf/2)-1; % interpolated array offset
+  mYs = [ mYs(1)*ones(roffset,1); mYs(1:Ns/2+1-roffset) ];
+  pYs = 2*pi*rand(Ns/2+1,1);                 % generate phase random values 
+  mYs1 = [mYs(1:Ns/2+1); mYs(Ns/2:-1:2)];    % create magnitude spectrum 
+  pYs1 = [pYs(1:Ns/2+1); -1*pYs(Ns/2:-1:2)]; % create phase spectrum 
+  Ys = mYs1.*cos(pYs1)+1i*mYs1.*sin(pYs1);   % compute complex spectrum 
+  yhw = fftshift(real(ifft(Yh)));            % sines in time domain using IFFT 
+  ysw = fftshift(real(ifft(Ys)));            % stoc. in time domain using IFFT 
+  yh(ri:ri+Ns-1) = yh(ri:ri+Ns-1)+yhw(1:Ns).*sw;   % overlap-add for sines 
+  ys(ri:ri+Ns-1) = ys(ri:ri+Ns-1)+ysw(1:Ns).*sws;  % overlap-add for stoch. 
+  pin = pin+H;                                     % advance the sound pointer 
+end 
+
+%ys=tanh(10*ys);
+y= yh+ys; % sum sines and stochastic
\ No newline at end of file
diff -r 59c6861155d9 -r 881c3acf1164 code Kat/spsmodel.m
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/code Kat/spsmodel.m	Sat Apr 20 12:35:51 2013 +0100
@@ -0,0 +1,125 @@
+function [y,yh,ys] = spsmodel(x,fs,w,N,t,maxnS,stocf)
+%=> analysis/synthesis of a sound using the sinusoidal harmonic model
+% x: input sound, fs: sampling rate, w: analysis window (odd size),
+% N: FFT size (minimum 512), t: threshold in negative dB,
+% maxnS: maximum number of sinusoids,
+% stocf: decimation factor of mag spectrum for stochastic analysis
+% y: output sound, yh: harmonic component, ys: stochastic component
+M = length(w); % analysis window size
+Ns = 1024; % FFT size for synthesis
+H = 256; % hop size for analysis and synthesis
+N2 = N/2+1; % half-size of spectrum
+soundlength = length(x); % length of input sound array
+hNs = Ns/2; % half synthesis window size
+hM = (M-1)/2; % half analysis window size
+pin = max(N2,1+hM) % initialize sound pointer to middle of analysis window
+pend = soundlength-max(N2,hM); % last sample to start a frame
+fftbuffer = zeros(N,1); % initialize buffer for FFT
+yh = zeros(soundlength+Ns/2,1); % output sine component
+ys = zeros(soundlength+Ns/2,1); % output residual component
+w = w/sum(w); % normalize analysis window
+sw = zeros(Ns,1);
+ow = triang(2*H-1); % overlapping window
+ovidx = Ns/2+1-H+1:Ns/2+H; % overlap indexes
+sw(ovidx) = ow(1:2*H-1);
+bh = blackmanharris(Ns); % synthesis window
+bh = bh ./ sum(bh); % normalize synthesis window
+wr = bh; % window for residual
+sw(ovidx) = sw(ovidx) ./ bh(ovidx);
+sws = H*hanning(Ns)/2; % synthesis window for stochastic
+lastysloc = zeros(maxnS,1); % initialize synthesis harmonic locations
+ysphase = 2*pi*rand(maxnS,1); % initialize synthesis harmonic phases
+fridx = 0;
+
+
+
+
+
+while pin<pend
+    %-----analysis-----%
+    xw = x(pin-hM:pin+hM).*w(1:M); % window the input sound
+    fftbuffer = fftbuffer*0; % reset buffer;
+    fftbuffer(1:(M+1)/2) = xw((M+1)/2:M); % zero-phase window in fftbuffer
+    fftbuffer(N-(M-1)/2+1:N) = xw(1:(M-1)/2);
+    X = fft(fftbuffer); % compute the FFT
+    mX = 20*log10(abs(X(1:N2))); % magnitude spectrum
+    pX = unwrap(angle(X(1:N/2+1))); % unwrapped phase spectrum
+    ploc = 1 + find((mX(2:N2-1)>t) .* (mX(2:N2-1)>mX(3:N2)) ...
+        .* (mX(2:N2-1)>mX(1:N2-2))); % find peaks
+    [ploc,pmag,pphase] = peakinterp(mX,pX,ploc); % refine peak values
+    % sort by magnitude
+    [smag,I] = sort(pmag(:),1,'descend');
+    nS = min(maxnS,length(find(smag>t)));
+    sloc = ploc(I(1:nS));
+    sphase = pphase(I(1:nS));
+    if (fridx==0)
+        % update last frame data for first frame
+        lastnS = nS;
+        lastsloc = sloc;
+        lastsmag = smag;
+        lastsphase = sphase;
+    end
+    % connect sinusoids to last frame lnS (lastsloc,lastsphase,lastsmag)
+    sloc(1:nS) = (sloc(1:nS)~=0).*((sloc(1:nS)-1)*Ns/N); % synth. locs
+    lastidx = zeros(1,nS);
+    for i=1:nS
+        [dev,idx] = min(abs(sloc(i) - lastsloc(1:lastnS)));
+        lastidx(i) = idx;
+    end
+    ri= pin-hNs; % input sound pointer for residual analysis
+    xr = x(ri:ri+Ns-1).*wr(1:Ns); % window the input sound
+    Xr = fft(fftshift(xr)); % compute FFT for residual analysis
+    Xh = genspecsines(sloc,smag,sphase,Ns); % generate sines
+    Xr = Xr-Xh; % get the residual complex spectrum
+    mXr = 20*log10(abs(Xr(1:Ns/2+1))); % magnitude spectrum of residual
+    mXsenv = decimate(max(-200,mXr),stocf); % decimate the magnitude spectrum
+    % and avoid -Inf
+    %-----synthesis data-----%
+    ysloc = sloc; % synthesis locations
+    ysmag = smag(1:nS); % synthesis magnitudes
+    mYsenv = mXsenv; % synthesis residual envelope
+    
+    
+%-----transformations-----%
+
+%-----time mapping-----%
+outsoundlength = soundlength*2;
+TM=[ 0 soundlength/fs; % input time (sec)
+0 outsoundlength/fs ]; % output time (sec)
+
+%%-----frequency stretch-----%
+%fstretch = 1.1;
+%ysloc = ysloc .* (fstretch.^[0:length(ysloc)-1]');
+   
+%-----synthesis-----%
+    if (fridx==0)
+        lastysphase = ysphase;
+    end
+    if (nS>lastnS)
+        lastysphase = [ lastysphase ; zeros(nS-lastnS,1) ];
+        lastysloc = [ lastysloc ; zeros(nS-lastnS,1) ];
+    end
+    ysphase = lastysphase(lastidx(1:nS)) + 2*pi*(lastysloc(lastidx(1:nS))+ysloc)/2/Ns*H; % propagate phases
+    lastysloc = ysloc;
+    lastysphase = ysphase;
+    lastnS = nS; % update last frame data
+    lastsloc = sloc; % update last frame data
+    lastsmag = smag; % update last frame data
+    lastsphase = sphase; % update last frame data
+    Yh = genspecsines(ysloc,ysmag,ysphase,Ns); % generate sines
+    mYs = interp(mYsenv,stocf); % interpolate to original size
+    roffset = ceil(stocf/2)-1; % interpolated array offset
+    mYs = [ mYs(1)*ones(roffset,1); mYs(1:Ns/2+1-roffset) ];
+    mYs = 10.^(mYs/20); % dB to linear magnitude
+    pYs = 2*pi*rand(Ns/2+1,1); % generate phase spectrum with random values
+    mYs1 = [mYs(1:Ns/2+1); mYs(Ns/2:-1:2)]; % create complete magnitude spectrum
+    pYs1 = [pYs(1:Ns/2+1); -1*pYs(Ns/2:-1:2)]; % create complete phase spectrum
+    Ys = mYs1.*cos(pYs1)+1i*mYs1.*sin(pYs1); % compute complex spectrum
+    yhw = fftshift(real(ifft(Yh))); % sines in time domain using inverse FFT
+    ysw = fftshift(real(ifft(Ys))); % stochastic in time domain using IFFT
+    yh(ri:ri+Ns-1) = yh(ri:ri+Ns-1)+yhw(1:Ns).*sw; % overlap-add for sines
+    ys(ri:ri+Ns-1) = ys(ri:ri+Ns-1)+ysw(1:Ns).*sws; % overlap-add for stochastic
+    pin = pin+H; % advance the sound pointer
+    fridx = fridx+1;
+end
+y= yh+ys; % sum sines and stochastic
\ No newline at end of file