Code covered by the BSD License  

Highlights from
NICAM

NICAM

by

 

Near Instantaneously Companded Audio Multiplex

NICAM(x,fs)
function [y,x14] = NICAM(x,fs)
% Near Instantaneously Companded Audio Multiplex
% 
%   y = NICAM(x,fs)
%   [y,x14] = NICAM(...)
% 
%   y = NICAM(x,fs) applies the NICAM audio bit-rate
%   reduction algorithm to the input column vector or matrix
%   x sampled at fs Hz. If x is a matrix, then the NICAM
%   algorithm is applied down the columns of x. The input x
%   should be in the range [-1,1]. The NICAM algorithm has
%   the following steps:
% 
%   1. Convert samples to 14-bit accuracy.
%   2. Divide the signal into 1 ms frames.
%   3. Define the "coding range" for each frame from the
%      largest (absolute) sample in the frame.
%   4. Remove a number of LSBs for each sample in a frame
%      according to the coding range for the frame. The
%      number of LSBs is:
%       - 4 if coding range 1 (max abs amplitude is
%         1XXXXXXXXXXXXX);
%       - 3 if coding range 2 (max abs amplitude is
%         01XXXXXXXXXXXX);
%       - 2 if coding range 3 (max abs amplitude is
%         001XXXXXXXXXXX);
%       - 1 if coding range 4 (max abs amplitude is
%         0001XXXXXXXXXX);
%       - 0 if coding range 5 (max abs amplitude is
%         00001XXXXXXXXX).
%       (Note that this notation uses binary integer rather
%       than twos complement)
%   5. Removed bits are replaces with intermediate values
%   and samples are converted back to 14-bit.
% 
%   [y,x14] = NICAM(...) returns the 14-bit version of the
%   input. The algorithm re-dithers the audio, using TPDF
%   dither, before it is companded to 14-bit resolution.

%  (c) Chris Hummersone & University of Surry, 2013
%  c.hummersone@surrey.ac.uk

% check input
assert(nargin==2,'Incorrect number of inputs')
assert(iscolumn(x) | ismatrix(x),'x must be a column vector or matrix')
assert(max(abs(x(:)))<=1,'x must be in the range [-1,1]')

y = zeros(size(x));
x14 = zeros(size(x));

for n = 1:size(x,2)
    [y(:,n),x14(:,n)] = do_nicam(x(:,n),fs);
end

function [y,x14] = do_nicam(x,fs)

x = double(x); % ensure input is double
N = 14; % NICAM resolution

% convert to 14 bit integer (but keep as double); re-dither
dither = rand(size(x))-rand(size(x)); % TPDF
x14 = round((x.*(2^(N-1)))+dither);
% clip if dither exceeds permitted values
x14max = (2^(N-1))-1;
IX = abs(x14)>x14max;
x14(IX) = sign(x14(IX)).*x14max;

% store polarity of x (i.e. sign bit carried separately)
polarity = x14<0;

% convert to binary of N-1 bits (input must be nonnegative)
x14 = convert_to_N_bit(abs(x14),N-1);

% convert back to int
x14 = bin2dec(x14);

% restore polarity
x14(polarity) = -x14(polarity);

% frame length (1 ms) in samples
frame_length = round(0.001*fs);

% calculate number of frames and zero pad input if not whole
% number of frames
num_frames = ceil(length(x14)/frame_length);
x14 = [x14; zeros((num_frames*frame_length)-length(x14),1)];

% reshape x to make each row a frame
z = reshape(x14,[frame_length length(x14)/frame_length])';

% convert peak value in each frame (row) to binary
peak = convert_to_N_bit(max(abs(z),[],2),N-1);
peak_ones = strfind(cellstr(peak),'1');
coding_range = 5.*ones(num_frames,1); % default is 5 (no compansion)
IX = ~cellfun(@isempty,peak_ones);

% Determine coding range for each frame
coding_range(IX) = cellfun(@(x) x(1),peak_ones(IX));
coding_range(coding_range>5) = 5;

%% Do the companding

% convert to binary
z = convert_to_N_bit(abs(x14),N-1);

% apply compansion to each frame
for j = 1:num_frames
    samples = ((j-1)*frame_length)+1:j*frame_length;
    switch coding_range(j)
        case 1
            z(samples,end-3:end) = repmat('0111',frame_length,1);
        case 2
            z(samples,end-2:end) = repmat('011',frame_length,1);
        case 3
            z(samples,end-1:end) = repmat('01',frame_length,1);
        case 4
            z(samples,end) = repmat('0',frame_length,1);
    end
end

% convert back to int, crop, and correct polarity
y = bin2dec(z);
x14 = x14(1:length(x));
y = y(1:length(x));
y(polarity) = -y(polarity);

function y = convert_to_N_bit(x,N)
% convert a signal to N-bit binary

assert(iscolumn(x),'x must be a column vector')

y = dec2bin(x,N);

assert(size(y,2)==N,'Input requires more bits than were specified')

Contact us