E71 DSP - Lab 1: MATLAB Review

Allison Barlow and Alex Benn

Due 30 January 2008


Part 1: Delayed Sample Function

Here we present the code and resulting graph.

% E71 Lab 1
% Allison Barlow, Alex Benn
% Due Wed, 01/30/08

%% Part 1 - Delayed Sample Function

N = 8; % length of sample sequence
M = 4; % delay

% d[n], general form of the unit sample of length N:
d = [1 zeros(1,N-1)];

% dm[n], general form of delayed sample sequence (d[n] delayed by M samples):
dm = [zeros(1,M) 1 zeros(1,N-M-1)]; % TODO

% clear plot
clf;

% plot d[n]
subplot(2,1,1), stem(d);
axis([1,10,0,1.5]);
title('Unit Sample Function, N=8');
xlabel('n');
ylabel('Unit Sample');

% plot dm[n]
subplot(2,1,2), stem(dm);
axis([1,10,0,1.5]);
title('Delayed Unit Sample Function, N=8, M=4');
xlabel('n');
ylabel('Delayed Unit Sample');

Figure 1: Unit sample function with and without delay

Part 2 - More Simple Functions

% E71 Lab 1
% Allison Barlow, Alex Benn
% Due Wed, 01/30/08

%% Part 2 - More Simple Functions
N = 10; % length of sample sequence

% d[n], general form of the unit sample of length N:
d = [1 zeros(1,N-1)];

% ds[n], unit step:
ds = [ones(1,N)];

% dr[n], unit ramp:
dr = [zeros(1,N)];
for i=1:N
    dr(i)=i-1;
end

% clear plot
clf;

% plot d[n]
hold on;
stem(dr, 'g');
stem(ds, 'b');
stem(d, 'r');
axis([1,10,0,10]);
xlabel('n');

Figure 2: Three Simple Functions

Part 3 - Sawtooth and Square Functions

% E71 Lab 1
% Allison Barlow, Alex Benn
% Due Wed, 01/30/08

%% Part 3
L = input('Input L:');
A = input('Input A:');
N = input('Input N:');

S = 2*pi/N;
T = 0:S:(L-1)*S;
X = 0:L-1; % graph against this X axis, not the default

% define sawtooth and square wave arrays
dsaw = sawtooth(T)*A;
dsquare = square(T)*A;

% clear plot
clf;

% plot d[n]
hold on;
subplot(2,1,1), stem(X,dsaw);
xlabel('Time');
ylabel('Amplitude');

subplot(2,1,2), stem(X,dsquare);
xlabel('Time');
ylabel('Amplitude');

Figure 3: Sawtooth and Square Waves

Part 4 - Aliasing

Below is the code and resulting graph. Notice that the sampled signals of the original 3Hz and 13Hz sinusoids are identical, and the 7Hz sampled signal is simply the additive inverse of the 3Hz and 13Hz samples. This is a result of aliasing: the 7Hz and 13Hz signals are above the Nyquist frequency of 10Hz/2 = 5Hz, and therefore cannot be uniquely represented by a 10Hz-sampled signal.

% E71 Lab 1
% Allison Barlow, Alex Benn
% Due Wed, 01/30/08

%% Part 4
T = 0:.001:1;
Ts = 0:.1:1;
X = 0:L-1; % graph against this X axis, not the default

% define sine waves
S1 = sin(2*pi*3*T);
S2 = sin(2*pi*7*T);
S3 = sin(2*pi*13*T);

% sampled sine waves
S1_samp = sin(2*pi*3*Ts);
S2_samp = sin(2*pi*7*Ts);
S3_samp = sin(2*pi*13*Ts);

% clear plot
clf;

% plot d[n]
hold on;
subplot(4,1,1);
hold on;
plot(T,S1,'r');
plot(T,S2,'g');
plot(T,S3,'b');
xlabel('Time');
ylabel('Amplitude');
title('Three sinusoids with frequencies of 3 Hz, 7 Hz and 13 Hz');

subplot(4,1,2);
hold on;
stem(Ts,S1_samp,'r');
plot(T,S1,'--','Color',[.8,.8,.8]);
xlabel('Time');
ylabel('Amplitude');
title('3 Hz Sinusoid Sampled at 10 Hz');

subplot(4,1,3);
hold on;
stem(Ts,S2_samp,'g');
plot(T,S2,'--','Color',[.8,.8,.8]);
xlabel('Time');
ylabel('Amplitude');
title('7 Hz Sinusoid Sampled at 10 Hz');

subplot(4,1,4);
hold on;
stem(Ts,S3_samp,'b');
plot(T,S3,'--','Color',[.8,.8,.8]);
xlabel('Time');
ylabel('Amplitude');
title('13 Hz Sinusoid Sampled at 10 Hz');

Figure 4: Three sinusoid waves sampled at 10Hz

Part 5 - Quantization Error

Below are the calculated SQNR values, code, and resulting graph for Question 1.16. Notice that error is larger and SQNR values are lower for truncation than for round-to-nearest; compare Figures 5 and 6 and the SQNR values listed in Table 1.

Table 1. SQNR Values for Different Rounding Methods
bits Truncate Round Theoretical
6 31.8223 36.0690 37.88
7 40.4275 42.1022 43.90
8 46.8486 49.4231 49.92
SQNR Equation for Truncation:
SQNR = -12.89 + 7.51*b
This equation was calculated using a linear fit regression on the calculated truncation values.

% E71 Lab 1
% Allison Barlow, Alex Benn
% Due Wed, 01/30/08

%% Part 5
%Quantize a signal to n bits.  This code assumes the signal is between -1
%and +1.

b=6; %Number of bits;
m=200; %Number of samples;
f0=1/50; %sin frequency
n=1:m; %n vector!

x=sin(2*pi*f0*n);   %signal between -1 and 1.
                                %Trying "sin()" instead of "sawtooth"
                                %results in more interesting error(to the
                                %extent that error is interesting).
                                
%When we quantize our signal, we must choose whether to truncate (round 
%toward zero) or round to the nearest quantized value. Swap the commenting
%of the following two lines to switch between truncation and nearest.

q = quantizer('fixed', 'nearest', 'saturate', [b, b-1]); %round to nearest
%q = quantizer('fixed', 'fix', 'saturate', [b, b-1]);    %truncate

xq = quantize(q,x-1/(2^b))+1/(2^b);

xe=x-xq;                   %Error

Pq = 1/m * sum(xe.^2);
Px = 1/m * sum(x.^2);

SQNR = 10*log10(Px/Pq)
clf;
subplot(2,1,1);
stem(x,'b');
hold on;
stem(xq,'r');
legend('exact','quantized','Location','Southeast')
subplot(2,1,2);
hold on;
stem(xe,'g');
title(sprintf('Signal, Quantized signal and Error for %g bits, %g quantization levels',b,2^b));
hold off

Figure 5: Quantization using truncation

Figure 6: Quantization using round-to-nearest

Part 6 (Extra) - Elegance in Quantization

We modified Prof. Cheever's quantization code to make it cleaner and more readable. Our modified version is presented below, along with the resulting graph to demonstrate correctness.

% E71 Lab 1
% Allison Barlow, Alex Benn
% Due Wed, 01/30/08

%% Part 6 - Extra
%Quantize a signal to n bits.  This code assumes the signal is between -1
%and +1.

n=3;                       %Number of bits;
m=120;                     %Number of samples;

x=sawtooth(2*pi*(0:(m-1))/m);   %signal between -1 and 1.
                                %Trying "sin()" instead of "sawtooth"
                                %results in more interesting error(to the
                                %extent that error is interesting).
x(find(x>=1))=(1-eps);     %Make  signal from -1 to just less than 1.

q = quantizer('fixed', 'nearest', 'saturate', [n, n-1]); %define quantizer
xq = quantize(q,x-1/(2^n))+1/(2^n);


xe=x-xq;                   %Error

stem(x,'b');
hold on;
stem(xq,'r');
hold on;
stem(xe,'g');
legend('exact','quantized','error','Location','Southeast')
title(sprintf('Signal, Quantized signal and Error for %g bits, %g quantization levels',n,2^n));
hold off

Figure 7: Quantization on Sawtooth Signal

Figure 8: Quantization on Sine Signal