Frequency Response Creator

Ron 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

Ron 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

Ron 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

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

[Y,FS,NBITS]=wavread('input.wav');
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

Ron 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

Ron 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

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

Signal-to-Distortion SDR Ratio

Ron 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

Ron April 1, 2011 Coded in Matlab
clf;
x = imread('input.jpg');
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

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

clf;
x = imread('input.jpg');
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