%Lecture_042312.m

dt = .01;  %step size
maxt = 60;  
t = 0:dt:maxt;
nt =length(t);
k =2;

n = 3;  %number of cascades

impulse = zeros(size(t));
impulse(1) = 1/dt;
h = impulse;
for i=1:n
    h = leakyIntegrator(h,k,t)/k;
end

figure(1)
clf
plot(t,h,'b-');


%% Uncomment one of these 
% random h
% h = randn(size(h));
% h(t>10) = 0;

% delay
% h = zeros(size(t));
% h(600) = 1/dt;

% difference of Gammas
% k1 = 2;
% n1 = 3;
% h1 = Gamma(n1,k1,t);
% 
% k2= 2;
% n2 = 5;
% h2 = Gamma(n2,k2,t);
% h = h1-h2;

% sin X exponential
% freq = .1;
% h = sin(2*pi*t*freq).*exp(-.1*t);

%Gabor!
freq = .1;
sig = 10;
t0 =   maxt/6;
h = sin(2*pi*t*freq).* exp( -((t-t0).^2/(2*sig)));

% Gamma
% h = Gamma(n,k,t);

clf
hold on
plot(t,h,'r-');

%% response to sinusoids

freq = 1/20;  %Hz

y = sin(2*pi*freq*t);

figure(2)
clf
subplot(2,1,1)
plot(t,y,'b-');


resp = conv(y,h)*dt;
resp = resp(1:length(t));

subplot(2,1,2)
plot(t,resp,'r-');
set(gca,'YLim',[-1,1]);

%% frequency response function

nCycles = [1:30];  

freqs = nCycles/maxt;

amps = zeros(1,length(nCycles));
for i=1:length(nCycles)
    y = sin(2*pi*nCycles(i)*t/maxt);
    resp = conv(y,h)*dt;
    resp = resp(1:length(t));
    hi = max(resp(t>10));
    lo = min(resp(t>10));
    amps(i) = (hi-lo)/2;
end

figure(3)
clf
stem(freqs,amps);
xlabel('frequency (Hz)');
ylabel('response');



%% FFt of the impulse response function

H = fft(h)*dt;

Hamp = abs(fft(h))*dt;

hold on
plot(freqs,Hamp(nCycles+1),'r.');



%% The Fourier Theorem

%random input

%y = randn(size(t));

y =double(t>20& t<30)/dt;



resp = conv(y,h)*dt;
resp = resp(1:length(t));
figure(4)
clf
subplot(2,1,1)
plot(t,y);
subplot(2,1,2)
plot(t,resp);


%% FFT of y and resp

Y = complex2real(fft(y),t);
R = complex2real(fft(resp),t);

figure(5)
clf
plot(Y.freq(nCycles),Y.amp(nCycles),'ro-');
hold on
plot(R.freq(nCycles),R.amp(nCycles),'go-');

%% linear filtering through the FFT


resp2 =        ifft(fft(h).*fft(y))*dt;

figure(4)
subplot(2,1,2)
hold on
plot(t,resp2,'r.');









