%Computes the transfer function (gain and phase angle) for a 2-stage
%IDEAL 1-pole LPF

% The last part of this code shows you how to overlay actual data points on
% the theoretical plots (very handy for lab reports!)

% You may wish to modify or adapt this code for similar problems.
% Code author: JE. W&L Phys-Engn, Oct 10 2017.

define angular frequencies, omega. For computations speed I space these

use logarithmically spaced values since we are always doing plots vs. log10(omega) but you can choose whichever values you like.

w = logspace(-5,8,10000);

Enter circuit parameter values here

%first stage
R1 = 50e3; %Ohm
C1 = 10e-6; %F
w1 = 1/(R1*C1); %define omega_1

%high pass elements
R2 = 500e3; %Ohm
C2 = 1e-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)

choose a value for input. I choose 1 for convenience, since this makes Vout/Vin = Vout

Vin = 1;

% Enter transfer function numerator and denominator (divide and conquer)
Numer = 1;  %numerator
Denom =(1+ 1i*w/w1); %1-stage lpf.

absH = abs(Numer./Denom);  %abs(...) computes magnitude.

% the ./ is a trick to divide two vectors element-wise. You could choose
% instead just to loop over each value of omega above.

%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. Transfer function H(w) = Vout/Vin vs log(w)
figure(1); hold on;
plot(log10(w), absH, 'b', 'linewidth', 2)

xlabel('log_{10}(\omega)', 'fontsize', 20)
 ylabel('|H(\omega)|', 'fontsize', 20)
xlim([-5 4])

%2. Decibel Gain: G(w) = 20 log10 (|H(w)|) vs. log w
figure(2);hold on;
plot(log10(w), 20*log10(absH), 'b', 'linewidth', 2)
xlabel('log_{10}(\omega)', 'fontsize', 20)
ylabel('G(\omega) (dB)', 'fontsize', 20)
ylim([-120 3])
xlim([-5 4])
grid on

% 3. Phase angle: phi(omega) vs log(w)
figure(3); hold on;
plot(log10(w), rad2deg(Phi), 'b', 'linewidth', 2)
xlabel('log_{10}(\omega)','fontsize', 20)
ylabel('\phi (deg)', 'fontsize', 20)
xlim([-5 4])
grid on

Plot your data points overlaid the theoretical curves

% Step 1. Call up the figure you want to plot into.  In this case I'll plot decibel gain data
% This is contained in figure 2


%here's some sample data points to plot.
% You can IMPORT these directly from Excel using the Import Wizard in the
% Matlab main command window (toward upper left).  Take note of what your variable names.
% I'll assume they are G and logw here, like so:
logw = [-4.5, -4, -3.5, -2, -1, -0.6, -0.1, 0.2, 0.5, 0.7, 1, 2, 3];
G =    [-0.3, -0.2, -0.35, 0.5, -1, -3, -5, -4, -12, -14, -17 -31, -55];

% plot with fancy options
% 'rs' means, Red Squares.  the last part says fill in the face with
% Red also. To learn more type >> help plot  or >> doc plot
plot(logw, G, 'ks', 'markerfacecolor', 'k', 'markersize', 8);