function Vrms = replaygain(filename,a1,b1,a2,b2)
% Determine the perceived loudness of a file
% METHOD:
% 1) Calculate Vrms every 50ms
% 2) Sort in ascending order of loudness
% 3) Pick the 95% interval (i.e. go 95% up the list, and choose the value at this point)
% 4) Convert this value into dB
% 5) return this value.
% Back in the main program...
% 6) Subtract it from that calculated for -20dB FS RMS pink noise
% Result = required correction to replay gain (relative to 83dB reference)

% David Robinson, 10th July 2001. http://www.David.Robinson.org/


% Get information about file
lngth=wavread(filename,'size');
samples=lngth(1);
channels=lngth(2);
% Read sampling rate and No. of bits
[dummy,fs,bs]=wavread(filename,[1 2]);

% The the file isn't CD sample rate, try to generate appropriate equal loudness filter
if fs~=44100 | nargin<2,
   [a1,b1,a2,b2]=equalloudfilt(fs);
end

% Set the Vrms window to 50ms
rms_window_length=round(50*(fs/1000));

% Set the interval to 95%
percentage=95;	% Which rms value to take as typical of whole file

% Set amount of data (in seconds) which Matlab on my PC happily copes with at once
block_length=2;	% chunk data in from wave file in 2 second blocks - file less than this length will cause an error

% Determine how many rms value to calculate per block of data
rms_per_block=fix((fs*block_length)/rms_window_length);

% Check that the file is long enough to process in block_length blocks
if lngth<(fs*block_length),
   warning(['skipping ' filename ' because it is too short']);
   Vrms=0;
   Vrms_all=0;
   return
end

% Display a Waitbar to show user how far into file we are
wbh=waitbar(0,'Processing...');

% Loop through all the file in blocks a defined above
for audio_block=0:fix(samples/(fs*block_length))-1,
   % Update the waitbar display to reflect progress
   waitbar(audio_block/(fix(samples/(fs*block_length))-1));
   % Grab a section of audio
   inaudio=wavread(filename,[(fs*block_length*audio_block)+1 fs*block_length*(audio_block+1)]);
   % Filter it using the equal loudness curve filter:
   inaudio=filter(b1,a1,inaudio);
   inaudio=filter(b2,a2,inaudio);
   % Calculate Vrms:
   for rms_block=0:rms_per_block-1,
      % Mono signal: just do the one channel
      if channels==1,
         Vrms_all((audio_block*rms_per_block)+rms_block+1)=mean(inaudio((rms_block*rms_window_length)+1:(rms_block+1)*rms_window_length).^2);
      % Stereo signal: take average Vrms of both channels   
      elseif channels==2,
         Vrms_left=mean(inaudio((rms_block*rms_window_length)+1:(rms_block+1)*rms_window_length,1).^2);
         Vrms_right=mean(inaudio((rms_block*rms_window_length)+1:(rms_block+1)*rms_window_length,2).^2);
         Vrms_all((audio_block*rms_per_block)+rms_block+1)=(Vrms_left+Vrms_right)/2;
      end
   end
end

% Close the waitbar
close(wbh);

% Convert to dB
Vrms_all=10*log10(Vrms_all+10^-10);
% Sort the Vrms values into numerical order
Vrms_all=sort(Vrms_all);
% Pick the 95% value
Vrms=Vrms_all(round(length(Vrms_all)*percentage/100));

return