% Computes the transfer function (gain and phase angle) for a 1-stage HPF
% You may wish to modify or adapt this code for similar problems.
%
% Code author: JE. W&L Phys-Engn, Sep 28 2018.


%% Define angular frequencies, omega. For computations speed I space these
%essentially logarithmically, but you can choose whichever values you like.

w = logspace(1,6,1000);  %to learn more about the logspace() function, type at the matlab command prompt >> doc logspace  

%% Define circuits component values and 1/RC cutoff frequencies, for convenience
% 1st stage parameters
R1 = 1e3; %Ohm
C1 = 0.68e-6; %F
w1 = 1/(R1*C1); %define omega_1

%2nd (cascaded) elements 
R2 = 1e4; %Ohm
C2 = 0.068e-6; %F
w2 = 1/(R2*C2); %define omega_2



%% initialize output.  This is only to speed up the execution of code by
%declaring memory space.
absH = zeros(size(w)); %the magnitude of the transfer function
Phi = zeros(size(w)); %phase shift angle (in radians)

%% Define transfer function
% Below is a 1-stage HPF

Numer = 1j*w/w1;  %numerator of Vout/Vin = H
Denom =(1+ 1j*w/w1); %denominator 

% IF we added another stage this would be our transfer function for a
% 2-stage HPF
% Numer =  (1j*w/w1)  .*  (1j*w/w2);
% Denom = (1+ 1i*w/w1).*(1+ 1i*w/w2);

%% Compute the magnitude of the transfer function
absH = abs(Numer./Denom);  %abs(...) computes magnitude. 
% Note: the ./ is a matlab's way to divide two vectors element-wise (instead of matrix opertion).

%% For phase angle, compute angle of numerator and denominator, then subtract
phi_numer = atan2(imag(Numer), real(Numer));
phi_denom = atan2(imag(Denom), real(Denom));
Phi = phi_numer - phi_denom;



%% Plot the results

%---------- 1. Plot magnitude response Vout/Vin vs log(f)
figure(1); hold on;
set(gcf, 'color', 'w')
set(gca, 'fontsize', 18)

plot(log10(w/(2*pi)), absH, 'b', 'linewidth', 2)

xlabel('log_{10}f (Hz)', 'fontsize', 20, 'interpreter', 'tex'); %the interpreter specifies to use fancy looking math font. 
ylabel('|H(f)|', 'fontsize', 20,  'interpreter', 'tex') %note greek letters are backslash then the letter, e.g \omega

grid on
axis tight

%---------2. Decibel Gain
figure(2);hold on;
set(gcf, 'color', 'w')
set(gca, 'fontsize', 18)

plot(log10(w/(2*pi)), 20*log10(absH), 'b', 'linewidth', 2)
xlabel('log_{10}f (Hz)', 'fontsize', 20,  'interpreter', 'tex')
ylabel('G(f) (dB)', 'fontsize', 20,  'interpreter', 'tex')

grid on
axis tight

%------ 3. Plot Phase Response
figure(3); hold on;
set(gcf, 'color', 'w')
set(gca, 'fontsize', 18)

plot(log10(w/(2*pi)), rad2deg(Phi), 'b', 'linewidth', 2)
xlabel('log_{10}f (Hz)','fontsize', 20, 'interpreter', 'tex')
ylabel('\phi (deg)', 'fontsize', 20, 'interpreter', 'tex')

grid on
axis tight