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');
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
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
#define COMB_FILTER_ORDER 3 // Filter Order, i.e., M
#define DELAY_ARRAY_SIZE 9 // Number of (delayed) input samples to store
int filtersum; // Sum-of-products accumulator for calculating the filter output value
int delay[DELAY_ARRAY_SIZE]; // Array for storing (delayed) input samples
static int gain1 = 1;
static int gain2 = 4;
interrupt void isr() //Interrupt service routine
{
short i; // Loop counter
delay[0] = get_sample() / gain1; // Read filter input sample from ADC and scale the sample
filtersum = delay[0]; // Initialize the accumulator
for (i = COMB_FILTER_ORDER; i > 0; i--)
{
filtersum += delay[i]; // Accumulate array elements
delay[i] = delay[i-1]; // Shift delay elements by one sample
}
filtersum *= gain2; // Scale the sum-of-products
send_output(filtersum); // write filter output sample to DAC
return; // interrupt servicing complete
}
void main()
{
short i; // loop counter
for (i=0; i< DELAY_ARRAY_SIZE; i++) // Initialize delay elements with 0
{
delay[i] = 0;
}
init_all(); // global initialization
while(1); // infinite loop
// do nothing except wait for the interrupt
}
[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]);
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
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
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
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');
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
