 function y = vusoft(x, fs)
 % Copyright 2006, Bryce Lobdell
 % Parameters for the system:
 % Td = 1/fs/d => oversample by 8x to prevent aliasing.
 wn = 13.5119; eta = .8127; D = 8;Td = 1/fs/D; 
 % Parameters for the filter:
 B = Td^2 * wn^2 * [ 1 2 1 ];
 A = [ ( 4 + 4*eta*wn*Td + wn^2*Td^2) (-8 + 2*wn^2*Td^2)...
      (4 - 4*eta*wn*Td + wn^2*Td^2) ];  
 % Scale:
 scaling = pi / 2 / sqrt(600*0.001*2);
 % Upsample the input signal by 8x.
 x_u = resample( x, D, 1, 50 );
 % Apply the absolute value, and the filter.
 yl_u = scaling * filter(B,A,abs(x_u));
 % Downsample back to the original rate.
 yl = yl_u(1:D:end); y = 20*log10( yl );
 %end
