% 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.
%
clear % Remove all variables from previous programs
R = 6; % Quanitzation bits for simulated quantization
L = 2^(R-1); % Largest signed 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 = 600; %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] = audioread(filname);
s = s/max(abs(s)); % Set max amplitude to 1
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];
sound(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)])