%EA3
% 
%%---------------------------------------------------------
%%Example: Handel through a tone control
%%---------------------------------------------------------


clear all;
close all;

global A F;			%share these with function frate_fn which computes Xdot=AX+F

%%
%%System parameters
%more treble
L	=	0.00005;				%element 1 - mass
C	=	0.0001;				%element 2 - capacitor
R	=	.1;			%element 3 - resistor
%more bass
%L	=	0.0005;				%element 1 - mass
%C	=	0.0001;				%element 2 - capacitor
%R	=	.1
%
%
% Voltage source: Input signal is Handel's Messiah

%load handel;
load laughter;
Vin	=	y(1:16000);
sound(Vin/max(Vin),8192)
plot(Vin)
pause
%%
%State Variables:
%	X(1,:) = voltage across capacitor element 2
%	X(2,:) = current through inductor element 1

% Initial conditions quiescent.
X(1,1) = 0;
X(2,1) = 0;
%time parameters

%State matrix
A = [0,		1/C;
	-1/L,	-R/L];  % for use in equation Xdot = AX+F

%

% Integrate equations of motion using the Runge-Kutta Method
dt = 1/8192;
for i = 1:1:16000

		F = [0,	Vin(i)/L]';				% votage input source value;
        k1 = frate_fn(X(:,i));			% present rate
        k2 = frate_fn(X(:,i) + dt*k1/2);	% estimated mid-point rate
        k3 = frate_fn(X(:,i) + dt*k2/2);	% even better mid-point rate
        k4 = frate_fn(X(:,i) + dt*k3);	% excellent end-point rate

        X(:,i+1) = X(:,i) + dt*(k1/6 + k2/3 + k3/3 + k4/6);	% use weighted average of rates		

end
Vout	= R*X(2,:);
sound(Vout/max(Vout),8192)
figure
plot(Vout)
pause
%Show frequency content of input and output signals
figure
fVin = abs(fft(Vin));
subplot(2,1,1), plot(fVin(1:4000)/max(fVin))
fVout = abs(fft(Vout(1:16000)));
subplot(2,1,2), plot(fVout(1:4000)/max(fVout))
pause
%play input and output over again
for (kplay=1:1:2)
sound(Vin/max(Vin),8192)
sound(Vout/max(Vout),8192)
end