## Frequency Response Creator

April 11, 2011 Coded in Matlab
[header, s] = hdrload('data_file.txt');

num = 2;
x = zeros(ceil(44100/num),1);                 %channel 1
y = zeros(ceil(44100/num),1);                 %channel 2

Fs = 44100;                     % Sampling frequency
T = 1/Fs;                       % Sample time
L = Fs/num;                  	  % Length of window
g = floor(num*length(s)/Fs);    % number of windows
H = 0;
count = 50;
max_index = 1;                  % pointer for max value
max = 0;                        % max value
max_index2 = 1;                 % pointer for max value
max2 = 0;                       % max value

z = zeros(g,2);
z2 = zeros(g,2);

for i = 1:1:g-3

% creating windows for channel one and two
for j = 1:1:L
k = j + L*i;
x(j,1) = s(k+i,1);
y(j,1) = s(k+i,2);
end

NFFT = 2^nextpow2(L);            % Next power of 2 from length of y
YY = fft(y,NFFT)/L;
Y = fft(x,NFFT)/L;
f = Fs/2*linspace(0,1,NFFT/2+1); % creating frequency vector
Y2 = 2*abs(Y(1:NFFT/2+1));
YY2 = 2*abs(YY(1:NFFT/2+1));

%%channel 1
max = Y(max_index,1);

%%channel 2
if(max_index2 < length(YY2))
max2 = YY(max_index2,1);

b = length(YY2) - max_index2 - 2*count;
if b > 0
a = max_index2 + 2*count;
else
a = length(YY2);
end

for j = max_index2:1:(a-1)
if max2 < YY2(j+1,1)
max2 = YY2(j+1,1);
max = Y2(j+1,1);
max_index2 = j+1;
end
end
end
z2(i+1,1) = abs(max2);
z2(i+1,2) = f(1,max_index2);
z(i+1,1) = abs(max);
z(i+1,2) = f(1,max_index2);
if(max_index2 > 4*count)
max_index2 = max_index2 - count;
end

% RECORD MOVIE
% semilogy(z2(:,2),z2(:,1),f,YY2,z(:,2),z(:,1),f,Y2);
% axis([0 20000 10^-3 Inf]);

end

semilogy(z2(:,2),z2(:,1),z(:,2),z(:,1));
axis([0 20000 0.00001 Inf]);

## Step Response Pole Effect Demo 2

April 1, 2011 Coded in Matlab
for omega = 0.5:0.5:20
sys = tf([omega^2],[1,omega,omega^2]);
p = [1 omega omega^2];
r = roots(p);
hold on

subplot(1,2,1),plot(real(r), imag(r), 'o'), title('Pole Locations')
hold on
subplot(1,2,2),step(sys)
pause(0.5)
end

## Step Response Pole Effect Demo

April 1, 2011 Coded in Matlab
for zeta = 0.1:0.1:2
sys = tf([9],[1,6*zeta,9]);
p = [1 6*zeta 9];
r = roots(p);
hold on
subplot(1,2,1),plot(real(r), imag(r),'o'), title('Pole Locations'), xlabel('Real'), ylabel('Imaginary')
hold on
subplot(1,2,2),step(sys)
pause(.5)
end

## WAV 2D Color Spectrogram

April 1, 20111 comment Coded in Matlab
%FFT window size
window = 267;

zz = zeros(0,window);
Y = Y.*1000;
for m = window+1:window:length(Y)
y = Y(m-window:m);
yy = fft(y);
zz = [zz, yy];
end
colormap(jet(280));
image(abs(zz));
xlabel('Time')
ylabel('Frequency')

## DFT/IDFT Demo

April 1, 2011 Coded in Matlab
clf;
L = 100;
for (L_factor = 1:1:6)
%L_factor = 3;
N_factor = 1;
L = round(L*L_factor);
n = [0:1:L-1];
x = sin(pi/1000*(n.*n));

subplot(3,1,1);
stem(n,x);
title ('x[n]');

subplot(3,1,2);
y = fft(x,L*N_factor);
plot(abs(y));
title ('DFT');

subplot(3,1,3);
xr = ifft(y);
stem(n,xr(1:L));
title ('IDFT');

pause;
end

## Waterfall Spectrograph

April 1, 2011 Coded in Matlab
function waterfallspect(s, fs, spectsize, spectnum)
%WATERFALLSPECT Display spectrum of signal s as a 3-D plot.
%   s - input signal, use wavload to load a WAV file
%   fs - Sampling frequency (Hz).
%   spectsize - FFT window size
%   spectnum - number of windows to analyze

frequencies = [0:fs/spectsize:fs/2];
offset = floor((length(s)-spectsize)/spectnum);
for i=0:(spectnum-1)
start = i*offset;
A = abs(fft(s((1+start):(start+spectsize))));
magnitude(:,(i+1)) = A(1:spectnum/2+1);
end
waterfall(frequencies, 0:(spectnum-1), magnitude');
xlabel('frequency');
ylabel('time');
zlabel('magnitude');

## Image Quantizer Demo

April 1, 2011 Coded in Matlab
clf;
for (truncate = 0:8)
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))));
pause(0.5);
end

## Signal-to-Distortion SDR Ratio

April 1, 2011 Coded in Matlab
[f, fs, nbits] = wavread('input.wav');
[f2, fs2, nbits2] = wavread('input-truncated.wav');

top = 1/size(f,1)*(sum(f.^2))
bottom = 1/size(f2,1)*(sum(f2.^2))

Y = 10*log10(top/bottom)

## RGB (JPEG) Image Separator

April 1, 2011 Coded in Matlab
clf;
red = x(:,:,1);
green = x(:,:,2);
blue = x(:,:,3);

subplot(2,2,1);
imshow(x);
title('Original');
subplot(2,2,2);
imshow(red);
title('Red');
subplot(2,2,3);
imshow(green);
title('Green');
subplot(2,2,4);
imshow(blue);
title('Blue');

## Image Quantizer

April 1, 2011 Coded in Matlab
%number of bits to truncate
tred = 4;
tgreen = 2;
tblue = 3;

clf;
red = x(:,:,1);
green = x(:,:,2);
blue = x(:,:,3);

red_t = (red/2^tred) * 2^tred;
green_t = (green/2^tgreen) * 2^tgreen;
blue_t = (blue/2^tblue) * 2^tblue;

x_t(:,:,1) = red_t;
x_t(:,:,2) = green_t;
x_t(:,:,3) = blue_t;

imshow(x);
pause;
imshow(x_t)

size = (24 - tred - tgreen - tblue) / 24