
%==============================================================
% SMS-Matlab like emulation
%%==============================================================
clear all
close all
%==== USER DATA =====
DAFx_in = wavread('love.wav'); % wave file
%DAFx_in = wavread('ohwhere_22050.wav'); % wave file
SR = 22050; % sampling rate
w1Length = 1028; % analysis window size
n1 = 128; % analysis window hop size
nPeaks = 100; % number of peaks detected
nSines = 50; % number of sinuosoids to track (and synthetise)
minSpacePeaks = 2; % minimum space (bins) between two picked peaks
zp = 2; % zero-padding coefficient
rgain = 1.; % gain for the residual component
MaxFreq = 11000; % maximum frequency, in Hertz, for plottings
MinMag = -100; % minimum magnitude, in dB, for plottings
%--- figure data
fig5 = 'yes'; % for plotting the sine tracks
fig6 = 'yes'; % if uncommented, will plot the original signal,
% its sine and residual part, and the transformed signal
%=== Definition of the Windows ===
%--- definition of the analysis window
fConst=2*pi/(w1Length+1-1);
w1=[1:w1Length]';
w1=.35875 -.48829*cos(fConst*w1)+.14128*cos(fConst*2*w1) ...
-.01168*cos(fConst*3*w1);
w1=w1/sum(w1)*2;
N=w1Length*zp; % new size of the window
%--- synthesis window
w2=w1;
n2=n1;
%--- triangular window
wt2=triang(n2*2+1); % triangular window
%--- main lobe table of bh92
[bh92SINE2SINE,bh92SINE2SINEsize]=bh92SINE2SINEgeneration;
%--- data for the loops
frametime = n1/SR;
pin = 0;
pout = 0;
TuneLength=length(DAFx_in);
pend=TuneLength-w1Length;
%=== Definition of the data arrays ===
DAFx_in = [zeros(w1Length/2-n1-1,1); DAFx_in];
DAFx_outsine = zeros(TuneLength,1);
DAFx_outres = zeros(TuneLength,1);
%--- arrays for the partial tracking
iloc = zeros(nSines,1);
ival = zeros(nSines,1);
iphase = zeros(nSines,1);
previousiloc = zeros(nSines,1);
previousival = zeros(nSines,1);
maxSines = 400; % maximum voices for harmonizer
syniloc = zeros(maxSines,1);
synival = zeros(maxSines,1);
previoussyniloc = zeros(maxSines,1);
previousiphase = zeros(maxSines,1);
currentiphase = zeros(maxSines,1);
%--- arrays for the sinus' frequencies and amplitudes
SineFreq = zeros(nSines,ceil(TuneLength/n2));
SineAmp = zeros(nSines,ceil(TuneLength/n2));
pitch = zeros(1,1+ceil(pend/n1));
pitcherr = zeros(1,1+ceil(pend/n1));
%--- creating figures ---


if(exist('fig5'))
h = figure(5); %set(h,'position', [220, 45, 550, 200]);
axisFig5 = [1 1+ceil(pend/n1) 0 MaxFreq]; 
end

tic
%UUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUU
disp('analazing frame ...');
while pin<pend
%--- windowing
grain = DAFx_in(pin+1:pin+w1Length).*w1(1:w1Length);
%--- zero padding
padgrain = zeros(N,1);
padgrain(1:w1Length/2) = grain(w1Length/2+1:w1Length);
padgrain(N-w1Length/2+1:N) = grain(1:w1Length/2);
%--- fft computation
f = fft(padgrain);
r = abs(f);
phi = angle(f);
ft = r.*exp(j*phi);
%===== Analysis =====
%--- peak detection (and their plottings)
[ftloc, ftval]=PickPeaks(r(1:N/2),nPeaks,minSpacePeaks);
%--- calculate interpolated values (peak position,phase,amplitude)
[iftloc, iftphase, iftval] = ...
interpolatedValues (r,phi,N,zp,ftloc,ftval);
%--- pitch detection
[pitchvalue,pitcherror,isHarm] = ...
pitchDetection (r,N,SR,nPeaks,iftloc,iftval);
pitch(1+pin/n1) = pitchvalue*isHarm;
pitcherr(1+pin/n1) = pitcherror;
%--- peaks tracking
if (pin==0) %--- for the first frame
nNewPeaks = nSines;
else %--- creating new born tracks
for i=1:nSines
if (previousiloc(i)==0)
[previousiloc(i), previousival(i)] = CreateNewTrack ...
(iftloc, iftval, previousiloc, previousival, nSines, MinMag);
nNewPeaks = nNewPeaks - 1;
end
end
%--- simple Peak tracker
[iloc,ival,iphase,previousiloc,previousival,distminindex] = ...
peakTrackSimple(nSines,nPeaks,N,SR,pitchvalue,iftloc, ...
iftval,iftphase,isHarm,previousiloc,previousival);
end
%--- savings
previousival = ival;
previousiloc = iloc;
SineFreq(:,1+pin/n1)=max((iloc-1)/N*SR,0.);
% frequency of the partials
SineAmp(:,1+pin/n1)=max(ival, MinMag);
% amplitudes of the partials
syniloc(1:nSines) = max(1,iloc);
synival(1:nSines) = ival;
if(exist('fig3')) % plot: the trackings of partials
figure(3); clf; hold on
PlotTracking(SineFreq(:,1:1+pin/n1), pitch(1:1+pin/n1));
xlabel('Frame number');ylabel('Frequency (Hz)');
axis(axisFig3);title('Peak tracking'); drawnow
end
%--- residual computation
resfft = ft;
if(isHarm==1)
%alex
%resfft=resfft-sinefillSpectrum(iloc,ival,iphase,nSines,...
resfft=resfft-sinefillspectrum(iloc,ival,iphase,nSines,...
w1Length, zp, bh92SINE2SINE, bh92SINE2SINEsize);
end
nSynSines = nSines;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%===== Transformations =====
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%===== Synthesis =====
%--- phase computation
if (pin > 0)
for i=1:nSynSines
if (syniloc(i)~=0)
ifreq = (previoussyniloc(distminindex(i))+ syniloc(i))/2;
% average bin
freq = (ifreq-1)/N*SR; % freq in Hz (if loc=1 --> freq=0)
currentiphase(i)=unwrap2pi(previousiphase(distminindex(i))+...
2*pi*freq*frametime);
end
end
end
previoussynival = synival;
previoussyniloc = syniloc;
previousiphase = currentiphase;
%--- compute sine spectrum
padsynthft=sinefillspectrum(syniloc, synival, currentiphase,nSynSines,w1Length, zp, bh92SINE2SINE, bh92SINE2SINEsize);
if (isHarm==0)
padsynthft = zeros(size(padsynthft));
end
%--- residual computation
respadgrain=real(ifft(resfft));
resgrain=[respadgrain(N-w1Length/2+1:N); ...
respadgrain(1:w1Length/2)]./w2(1:w1Length);
ressynthgrain=wt2(1:n2*2).*resgrain(w1Length/2-n2:w1Length/2+n2-1); 
DAFx_outres(pout+1:pout+n2*2)=DAFx_outres(pout+1:pout+n2*2)+ ...
   ressynthgrain;
%--- sinusoidal computation
sinpadgrain=real(ifft(padsynthft));
singrain=[sinpadgrain(N-w1Length/2+1:N); ...
sinpadgrain(1:w1Length/2)]./w2(1:w1Length);
sinsynthgrain=wt2(1:n2*2).*singrain(w1Length/2-n2:w1Length/2+n2-1); 
DAFx_outsine(pout+1:pout+n2*2)=DAFx_outsine(pout+1:pout+n2*2)+ ...
sinsynthgrain;
%--- figure with original signal and transformed signal FFT
synthr = abs(fft(respadgrain + sinpadgrain));
if(exist('fig4'))
figure(4); clf; hold on
plot((1:N/2)/N*SR, 20*log10(r(1:N/2)),'b:'); axis(axisFig4);
plot((1:N/2)/N*SR, 20*log10(synthr(1:N/2)),'r');
figure(4);
xlabel('Frequency (Hz)');ylabel('Magnitude (dB)');axis(axisFig4);
title('FFT of the original (blue) and the')
title('transformed (red) signals');
drawnow
end
%--- increment loop indexes
pin = pin + n1;
pout = pout + n2;
disp(pin/n1);
end
%UUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUU
toc
%===== write output sounds =====
DAFx_in = DAFx_in(w1Length/2-n1:length(DAFx_in));
% remove the zeros added for the process
DAFx_outresynth = DAFx_outsine(1:TuneLength)+ ...
rgain*DAFx_outres(1:TuneLength);
mm = max(abs(DAFx_outresynth));
mi = max(abs(DAFx_in));
wavwrite(DAFx_outresynth/mm*mi, SR, 'DAFx_out.wav');
wavwrite(DAFx_outsine/mm*mi, SR, 'DAFx_outsine.wav');
wavwrite(DAFx_outres/mm*mi, SR ,'DAFx_outres.wav');

if(exist('fig6')) % plot the input signal, its sinus
% and its residual part, and the transformed signal
figure(6)
t = ((1:TuneLength)-1)/SR;
subplot(411)
plot(t, DAFx_in,'k');
ylabel('amplitud')
%gset title "Original Signal"
title('Original')
%figure(2)
subplot(412)
plot(t, DAFx_outresynth,'k');
%gset title "Resynthesized Signal"
ylabel('amplitud')
title('Resintesis')
%figure(3)
subplot(413)
plot(t, DAFx_outsine,'k');
%gset title "Sinusodal Synthesis"
ylabel('amplitud')
title('Sintesis sinusoidal')
%figure(4)
subplot(414)
plot(t, DAFx_outres,'k');
ylabel('amplitud')
xlabel('tiempo (s)')
%gset title "Residuum Synthesis"
title('Sintesis del residuo')
end

print -depsc 'sms_sintesis.eps'

if(exist('fig5'))% plot: trackings of partials
% only at the end of the process
figure(5); 
PlotTracking(SineFreq(:,1:round(1+pend/n1)), pitch(1:round(1+pend/n1)));
xlabel('Numero de trama'); ylabel('Frequencia (Hz)'); axis(axisFig5);
title('Seguimiento de picos');
axis tight
end

print -depsc 'sms_tracking.eps'
