% Usage: [ycomp,P] = audiocompress(y, m, r).  Given a row or column vector y,  
% blocksize m and scaling parameter r this command applied JPEG-like
% compression to the signal in blocks of size m, using the DCT and
% frequency domain quantization Y_k -> r*round(Y_k/r). Simulated compressed
% audio signal returned in vector "ycomp" (same dimensions as y). Fraction
% of non-zero block DCT coefficients returned in "P".

% If you don't have the signal processing toolbox substitute kdct() and
% kidct() for dct() and idct().

function [ycomp,P] = audiocompress(y, m, r)

% Get dimensions
[M0,N0] = size(y);
N = length(y);
blocks = ceil(N/m);

% Zero pad y to even blocksize.
y((N+1):(m*(blocks))) = 0.0;

% Reset N, set up new audio array for "compressed" signal
N = m*(blocks);
if M0 == 1
    ycomp = zeros(1,N);
else
    ycomp = zeros(N,1);
end

% DCT each block, then quantize as Y_k -> r*round(Y_k/r)
s = 0; %Counter for non-zero DCT coefficients
for k = 0:(blocks-1)
    y3 = y((k*m+1):((k+1)*m));
    Y = dct(y3);
    Y = round(Y/(r))*r; %Quantization step.
    s = s + sum(abs(Y)>0); %Keep count of non-zero coefficients
    ycomp((k*m+1):((k+1)*m)) = idct(Y);
end

% Restrict back to original size
ycomp = ycomp(1:max(M0,N0));

%Crude measure of compression ratio
P = s/N;