Why does 24-bit wave file introduce errors?

I have found that writing and reading 24-bit data to a .wav file introduces errors even though the format is uncompressed. This does not occur with 32-bit data. The following code

for bitsPerSample = [32 24]
    % Generate random integers that fit within the sample width
    B = 2^(bitsPerSample - 1);
    xlims = [-B, B-1];
    x = randi(xlims, [100 1], 'int32');

    % Write data to and read data from file
    audiowrite('tmp.wav', x, 8000, 'BitsPerSample', bitsPerSample);
    info = audioinfo('tmp.wav');
    y = audioread('tmp.wav', 'native');
    
    %  Compare written and read data
    fprintf('Bits per sample is %d\n', bitsPerSample)
    fprintf('CompressionMethod is ''%s''\n', info.CompressionMethod)
    fprintf('Std. deviation between written and read data is %.1f\n\n', ...
        std(double(x) - double(y)))
end

generates the following output in R2015b

Bits per sample is 32
CompressionMethod is 'Uncompressed'
Std. deviation between written and read data is 0.0

Bits per sample is 24
CompressionMethod is 'Uncompressed'
Std. deviation between written and read data is 77.6

What causes the precision loss?

Answer 1

[Thomas]

I deduced my own answer. When audiowrite is given int32 data to write as 24-bit samples, it considers the data to have 32 bits and drops the least significant 8 bits. Thus, 24-bit data should be multiplied by 2^8 before sending to audiowrite. Inserting a scaling statement into the code

for bitsPerSample = [32 24]
    % Generate random integers that fit within the sample width
    B = 2^(bitsPerSample - 1);
    xlims = [-B, B-1];
    x = randi(xlims, [100 1], 'int32');
    
    if bitsPerSample == 24, x = x * 2^8; end  % Scale 24-bit data to 32-bit limits
    
    % Write data to and read data from file
    audiowrite('tmp.wav', x, 8000, 'BitsPerSample', bitsPerSample);
    info = audioinfo('tmp.wav');
    y = audioread('tmp.wav', 'native');
    
    %  Compare written and read data
    fprintf('Bits per sample is %d\n', bitsPerSample)
    fprintf('CompressionMethod is ''%s''\n', info.CompressionMethod)
    fprintf('Std. deviation between written and read data is %.1f\n\n', ...
        std(double(x) - double(y)))
 end

eliminates the precision loss:

 Bits per sample is 32
CompressionMethod is 'Uncompressed'
Std. deviation between written and read data is 0.0

Bits per sample is 24
CompressionMethod is 'Uncompressed'
Std. deviation between written and read data is 0.0