% This script simulates the distortion from an R-bit quantizer and 
% computes SNR for the signal to quantization noise.
% The signal to be quantized is either a sine wave or signal from a
% wave file (typically quantized with 16-bit words).
% See comments in the script for how to set a flag to chose between
% the 2 signal options.  Since the wave file is scaled between -1 and 1 as
%  it is read in (i.e. divided by 2^15 if 16 bit quantization used), this script
% scales amplitudes back to the integer values so it will fit within the dyanmic
% range of the simulated quantizer.  Quantization (or requantization) is
% applied using R quantization  bits, and the signal is scaled up by 2^(R-1).
% Note that since one bit is used for the signed bit (in 2-complement
% format), the largest magnitude corresponds to 2^(R-1) (not 2^R). The
% script also plots the average spectra for both the original and
% quantized signals. The resulting SNR is printed in the plot title. 
%
%   Written by Kevin D. Donohue, 2-5-2008 ( donohue@engr.uky.edu )
%   updated 8-29-2010
clear     % Remove all variables from previous programs
R = 9;    % Quanitzation bits for simulated quantization

L = 2^(R-1);  % Largest quantization level magnitude (used for normalization)

flag = 1;  % Set flag equal to 1 use a sine wave for analysis
           % If equal to any other value it will read in a wave
           % file specified by the path/filenamein the second condition
           % of the following if statement
if flag ==1
    %  Create sinudoid
    dur = 2;  % Sine wave duration in seconds
    fs = 16e3;  %  Set sampling rate in Hz
    t = [0:round(2*fs)-1]'/fs;  % Time axis
    fo = 413.2015; %  Freqency of sinusoid
    g = sin(2*pi*fo*t);  % Generate sine wave
    s = g/max(abs(g));  % Normalized to ensure max absolute value is 1.
else
%  Read in wave file to be quanitzed (normalized amplitudes is part of the
%  wavread funtion)
    filname = 'mozart-1.wav';  % String denoting path and file name to be read in
                               % if no path is specifed, it will look for file
                               % in the current directory
    %  Read in signal (s), sampling frequency (fs) and bits per sample (B)
    [s, fs, B] = wavread(filname);
end

% Perform the quantization through scaling to integer values corresponding
% to the simulated number of quantization bit and rounding off 
sq = (ceil(s*(L-.5))-.5)/L;
%  Compute error relative to original signal (note, vector nq is the quantization noise signal) 
nq = sq-s;
%  compute signal to noise ratio
snr = 10*log10(mean(s.^2)/mean(nq.^2));
%  Concatenate sounds into one vector to hear one after the 
% other with a little pause/silence between them
stest = [s; zeros(round(fs*.1),1); sq];
soundsc(stest,fs)


% Plot the average spectra, called power spectral density (PSD), of the original
% and simulated quantized sound
wlen = 512;  %  Window length in samples for extracting data for FFT magnitude computation
olap = fix(wlen/2); % Number of samples (overlap) between sucessive windows
nfft = 2*wlen;  %  Number of FFT points (achieved through zero padding if longer than WLEN)
[ps,f] = pwelch(s,hamming(wlen),olap,nfft,fs);  % Compute PSD of original 
[psq,fq] = pwelch(sq,hamming(wlen),olap,nfft,fs); % Compute PSD of quantized
%  Plot both spectra for comparison
plot(f,10*log10(ps*fs/wlen),'--r',fq,10*log10(psq*fs/wlen),':k')
xlabel('Hertz')
ylabel('dB')
title(['Original (red) and quantized (black) signals: bits = ' int2str(R) ' SNR = ' num2str(snr)])