function y=vibrato(flute,fs,Modfreq,Width)
ya_alt=0;
[flute, fs, nbits] = wavread('flute_nonvib_F#5.wav');
Delay = 0.0003; % basic delay of input sample in sec (3 milliseconds recommended)
DELAY = round(Delay*fs); % basic delay in # samples
Width = 0.0003; %(recommended width)
WIDTH = round(Width*fs); % modulation width in # samples
if WIDTH>DELAY
error( 'delay greater than basic delay !!!');
return;
end;
Modfreq = 5; % (approx 5 Hz recommended)
MODFREQ=Modfreq/fs; % modulation frequency in # samples
LEN=length(flute); % # of samples in WAV-file
L=2+DELAY+WIDTH*2; % length of the entire delay
Delayline=zeros(L,1); % memory allocation for delay
y=zeros(size(flute)); % memory allocation for output vector320
for n=1:(LEN-1);
M=MODFREQ;
MOD=sin(M*2*pi*n);
ZEIGER=1+DELAY+WIDTH*MOD;
i=floor(ZEIGER);
frac=ZEIGER-i;
Delayline=[flute(n);Delayline(1:L-1)];
%---Linear Interpolation-----------------------------
y(n,1)=Delayline(i+1)*frac+Delayline(i)*(1-frac);
%---Allpass Interpolation------------------------------
%y(n,1)=(Delayline(i+1)+(1-frac)*Delayline(i)-(1-frac)*ya_alt);
%ya_alt=ya(n,1);

end

wavwrite(y, fs, 'flute_vibrato.wav');
