Copy quantize.m from aperitif_matlab git.

applications/lofar2/model/quantize.m

+%-----------------------------------------------------------------------------
+%
+% Copyright (C) 2016
+% ASTRON (Netherlands Institute for Radio Astronomy) <http://www.astron.nl/>
+% P.O.Box 2, 7990 AA Dwingeloo, The Netherlands
+%
+% This program is free software: you can redistribute it and/or modify
+% it under the terms of the GNU General Public License as published by
+% the Free Software Foundation, either version 3 of the License, or
+% (at your option) any later version.
+%
+% This program is distributed in the hope that it will be useful,
+% but WITHOUT ANY WARRANTY; without even the implied warranty of
+% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+% GNU General Public License for more details.
+%
+% You should have received a copy of the GNU General Public License
+% along with this program.  If not, see <http://www.gnu.org/licenses/>.
+%
+%-----------------------------------------------------------------------------
+% Author: E. Kooistra, 2016
+%
+% quantize - quantize two's-complement input data
+%     The output has nof_bits. Definitions:
+%     . w = nof_bits   (output data width)
+%     . q_full_scale = 2**(w-1)
+%     . q_max =  q_full_scale - 1
+%     . q_min = -q_full_scale
+%     Magnitude and amplitude are synonyms for full scale.
+%
+%     The quantization process involves the following steps:
+%
+%     1 scale the input data magnitude
+%     2 round the data
+%     3 handle overflow
+%     4 scale back to the original input data magnitude
+%
+%     The in_full_scale internally maps to in_q_max. Choose in_q_max =
+%     . q_full_scale : default, implies that in_data = in_full_scale will
+%                      just clip to q_max (or just wrap dependent on the
+%                      overflow option) and in_data = -in_full_scale will
+%                      map to q_min.
+%     . q_max        : in_data = in_full_scale will map to q_max, so no
+%                      overflow.
+%     Typically use in_q_max = q_full_scale for data and use in_q_max =
+%     q_max for coefficients. Use in_q_max < q_max if more backoff is
+%     needed for some other reason.
+%
+%     The 'rounding' determines how -0.5 is rounded:
+%     . 'half_away': Round half away from zero, so +0.5 --> 1, -0.5 --> -1.
+%     . 'half_up'  : Round half up to +infinity, so +0.5 --> 1, -0.5 --> 0.
+%
+%     Quantized data that is out of range -q_min : q_max gets treated 
+%     dependent on 'overflow':
+%     . 'clip'     : Clip to -q_min or q_max
+%     . 'clip_sym' : Clip symmetrical to -q_max or q_max
+%     . 'wrap'     : Wrap within nof_bits signed integer range
+%     . 'no_limit' : No range limit, so allow quantized data to get out
+%                    of range -q_min : q_max
+%
+%     The output data range is scaled back to the original in_full_scale.
+
+function out_data = quantize(in_data, in_full_scale, nof_bits, rounding, overflow, in_q_max)
+
+q_bit = 1;
+q_full_scale = 2^(nof_bits-1);  % maximum amplitude, magnitude
+q_max =  q_full_scale-1;
+q_min = -q_full_scale;
+q_period = 2^nof_bits;
+
+% Default options
+if ~exist('rounding', 'var'); rounding = 'half_away'; end;
+if ~exist('overflow', 'var'); overflow = 'clip'; end;
+if ~exist('in_q_max', 'var'); in_q_max = q_full_scale; end;
+
+% Scale
+q_data = in_q_max * in_data/in_full_scale;
+
+% Round
+if strcmp(rounding, 'half_away')
+    q_data = round(q_data);
+else
+    q_data = floor(q_data + q_bit/2);
+end
+
+% Overflow
+if isreal(q_data)
+    q_re = q_data;
+else
+    q_re = real(q_data);
+    q_im = imag(q_data);
+end
+if strcmp(overflow, 'clip')
+    q_re(q_re>q_max) = q_max;
+    q_re(q_re<q_min) = q_min;
+    if ~isreal(q_data)
+        q_im(q_im>q_max) = q_max;
+        q_im(q_im<q_min) = q_min;
+    end
+elseif strcmp(overflow, 'clip_sym')
+    q_re(q_re> q_max) =  q_max;
+    q_re(q_re<-q_max) = -q_max;
+    if ~isreal(q_data)
+        q_im(q_im> q_max) =  q_max;
+        q_im(q_im<-q_max) = -q_max;
+    end
+elseif strcmp(overflow, 'wrap')
+    q_re = mod(q_re - q_min, q_period) + q_min;
+    if ~isreal(q_data)
+        q_im = mod(q_im - q_min, q_period) + q_min;
+    end
+end
+if isreal(q_data)
+    q_data = q_re;
+else
+    q_data = complex(q_re, q_im);
+end
+% Back to original full scale
+out_data = in_full_scale * q_data/q_full_scale;