function vibrato(input_file, output_file, frequency)

% Reads a WAVE file specified by the parameter input_file, returning the effected sound in a file specified by the parameter output_file.
% Amplitude values are in the range [-1,+1].
% Fs is the sample rate in Hertz. NBITS is the number of bits per sample
% used to encode the data
[in_wave, Fs, NBITS ] = wavread('sample.wav');

ya_alt=0;
Delay = 0.0008; % basic delay of input sample in sec
Width = Delay;
%Modfreq = 2; %10 Khz
Modfreq=frequency;

% Basic delay in # samples
DELAY=round(Delay*Fs); 
% Modulation width in # samples
WIDTH=round(Width*Fs); 
if WIDTH>DELAY 
  error('delay greater than basic delay !!!');
  return;
end
% Modulation frequency in # samples
MODFREQ=Modfreq/Fs; 
% Number of samples in WAV-file
LEN=length(in_wave);        
% Length of the entire delay
L=2+DELAY+WIDTH*2;
% Memory allocation for delay
Delayline=zeros(L,1); 
% Memory allocation for output vector
out_wave = zeros(length(in_wave),1);     

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=[in_wave(n);Delayline(1:L-1)]; 
   % Linear Interpolation
   out_wave(n,1)=Delayline(i+1)*frac+Delayline(i)*(1-frac); 

end 

% Write vibrato output to vibrato_out file
wavwrite(out_wave, Fs, NBITS, output_file);

% Plot the before and after waves
figure(3)
hold on
plot(in_wave,'r');
plot(out_wave,'b');
% Plot first 500 samples of vibrato output waveform
%plot(x(1:500),'r');
%plot(y(1:500),'b');
title('Sample wave form and new vibrato wave form');