DSPRelated.com
   Code Snippets

2D Quantizer with Error Output

Ron ● April 1, 2011 ● Coded in Matlab 
clf;
truncate = 5;
x = imread('input.jpg');
subplot(2,2,1);
imshow(x);
subplot(2,2,2);
x_t = (x/2^truncate) * 2^truncate;
%if above line doesn't work your MATLAB is old! use next line:
% x_t = uint8(double(uint8((double(x)/2^truncate))) * 2^truncate);
imshow(x_t);
subplot(2,2,3);
imshow(uint8(abs(double(x_t) - double(x))));

WAV Quantizer with Error Output

Ron ● April 1, 2011 ● Coded in Matlab 
truncate = 6;
[f, fs, nbits] = wavread('in.wav');
f_int = int16(f*2^(nbits-1-truncate));
f_back_to_float = (double(f_int))/2^(nbits-1-truncate);
f_diff = f - f_back_to_float;
wavwrite(f_back_to_float, fs, nbits, 'out.wav');
wavwrite(f_diff, fs, nbits, 'error.wav');

Quantization Demo

Ron ● April 1, 2011 ● Coded in Matlab 
a = 10;
x = 0:2*pi/1000:2*pi;
y = a*sin(x);
y_int = int8(y);

subplot(3,1,1);
plot(x,y);
grid;
axis([0 2*pi -a*1.2 a*1.2]);
title('Continuous');

subplot(3,1,2);
plot(x,y_int);
grid;
axis([0 2*pi -a*1.2 a*1.2]);
title('Quantized');

subplot(3,1,3);
plot(x,abs(double(y_int)-y));
grid;
axis([0 2*pi 0 1.2/2]);
title('Error');

WAV Subsample

Ron ● April 1, 2011 ● Coded in Matlab 
factor = 8;
[X, fs, nbits] = wavread('in.wav');
Y = downsample(X,factor);
wavwrite (Y, fs/factor, nbits, 'out.wav');

2D Image Subsample

Ron ● April 1, 2011 ● Coded in Matlab 
A = imread('in.jpg');

M =   [ 1 4  6  4  1;
        4 16 24 16 4;
        6 24 36 24 6;
        4 16 24 16 4;
        1 4  6  4  1 ] / 256;

B = uint8(conv2(double(A),M));
imwrite(B,'out.jpg','JPG');

2D Image FIR Subsample

Ron ● April 1, 2011 ● Coded in Matlab 
A = imread('in.jpg');

M =   [ 1 4  6  4  1;
        4 16 24 16 4;
        6 24 36 24 6;
        4 16 24 16 4;
        1 4  6  4  1 ] / 256;

Y = filter2(M,A);
B = uint8(conv2(double(Y),M));

imwrite(B,'out.jpg','JPG');

WAV Butterworth Subsample

Ron ● April 1, 2011 ● Coded in Matlab 
[X, fs, nbits] = wavread('file.wav');
[N, Wn] = buttord(0.12, 0.125, 3, 60);
[B,A] = butter(N,Wn,'low');
Y = filter(B,A,X);
Z = downsample(Y,8);
wavwrite (Z, fs/8, nbits, 'file-sampled.wav');

Aliasing Demo

Ron ● April 1, 2011 ● Coded in Matlab 
clf;
% First plot:
subplot(2,2,1)
t = 0:0.005:10;
xa = 2*t.*exp(-t);
plot(t,xa);grid
xlabel('Time, sec');ylabel('Amplitude');
title('Continuous-time signal x_{a}(t)');

% Second plot:
subplot(2,2,2)
wa = 0:10/511:10;
ha = freqs(2,[1 2 1],wa);
plot(wa/(2*pi),abs(ha));grid;
xlabel('Frequency, Hz');ylabel('Magnitude');
title('|X_{a}(j\Omega)|');
axis([0 5/pi 0 2]);

% Third plot:
subplot(2,2,3)
T = .5;
n = 0:T:10;
xs = 2*n.*exp(-n);
k = 0:length(n)-1;
stem(k,xs);grid;
xlabel('Time index n');ylabel('Amplitude');
title('Discrete-time signal x[n]');

% Fourth plot:
subplot(2,2,4)
wd = 0:pi/255:4*pi;
hd = freqz(xs,1,wd);
plot(wd/(2*T*pi), T*abs(hd));grid;
xlabel('Frequency, Hz');ylabel('Magnitude');
title('|X_{\delta}(j\Omega)|');
axis([0 4/(2*T) 0 2])

Sampling Demo

Ron ● April 1, 2011 ● Coded in Matlab 
clf;
% First plot
t = 0:0.0005:1;
f = 5;
x_t = cos(2*pi*f*t);
subplot (4,1,1);
plot (t,x_t);
grid;
xlabel ('Time [sec]');
ylabel ('Amplitude');
title ('Continuous-time signal x(t)');
axis([0 1 -1.2 1.2])

% Second plot
n = 0:0.125:1;
x_n = cos(2*pi*f*n);
subplot(4,1,2);
y = zeros(1,length(t));
for i = 1:length(n)
    y = y + x_n(i)*sinc(t/0.125 - i + 1);
end
plot(t,y);
grid;
xlabel ('Time, sec');
ylabel ('Amplitude');
title ('Reconstructed continuous-time signal y(t) [fs < bound] (T=0.125)');
axis ([0 1 -1.2 1.2]);

% Third plot
n = 0:0.1:1;
x_n = cos(2*pi*f*n);
subplot(4,1,3);
y = zeros(1,length(t));
for i = 1:length(n)
    y = y + x_n(i)*sinc(t/0.1 - i + 1);
end
plot(t,y);
grid;
xlabel ('Time, sec');
ylabel ('Amplitude');
title ('Reconstructed continuous-time signal y(t) [fs = bound] (T=0.1)');
axis ([0 1 -1.2 1.2]);

% Fourth plot
n = 0:0.075:1;
x_n = cos(2*pi*f*n);
subplot(4,1,4);
y = zeros(1,length(t));
for i = 1:length(n)
    y = y + x_n(i)*sinc(t/0.075 - i + 1);
end
plot(t,y);
grid;
xlabel ('Time, sec');
ylabel ('Amplitude');
title ('Reconstructed continuous-time signal y(t) [fs > bound] (T=0.075)');
axis ([0 1 -1.2 1.2]);

Interpolation Demo

Ron ● April 1, 2011 ● Coded in Matlab 
clf;
for steps = 3:20
    x = 0:2*pi/steps:2*pi;
    y = sin(x);
    plot(x,y,'o-');
    grid;
    axis([0 2*pi -1.2 1.2]);
    pause(0.5)
end