%  This script will load a sample file obtained from a frequency sweep
%  Labview program of a high-pass filter and estimate the cut-off from
%  all available data using the phase of the transfer function:
%
%    vout         s/wc
%   ------   = -----------------------
%     vin         1    + s/wc
%
%  where s is complex frequency and wc is the cut-off frequency in 
%  radian per second. From the cut-off with a known resistor value
%  the capacitor value is also estimated through 1/(RC) = fc

r = 10.4e3;  %  Measured value of resistor
m = dlmread('hpfs1.txt');   %  Read in data file create by sweep program

% The format of matrix m is 4 columns where column 1 is the frequency in Hz,
% column 2 is the input peak to peak amplitude, column 3 is the output peak
% to peak and column 4 is the phase in degrees of input phase minus output
% phase.


%  Create vector of possible cut-off frequencies
%  with range over suspected region in small increments 
wc = 2*pi*[10:.1:500];  %  In radians per second

% Frequency points tested
f = m(:,1);   %  Frequency in Hz
s = j*2*pi*f;  %  Complex Frequency in radians per second
%  Loop to compute squared error for every iteration of test values

%  The measured transfer function phase (output-input, just the opposite
%  of format in matrix m, so must negate)
ms = -m(:,4);
% Frequency points tested
f = m(:,1);   %  Frequency in Hz
s = j*2*pi*f;  %  Complex Frequency in radians per second
%  Loop to compute squared error for every iteration of test values
err = [];  %  clear/initalize array to store error values in
for k=1:length(wc)
    testtf = 180*angle( (s/wc(k)) ./ (1+(s/wc(k))) )/pi;  %  Convert radians to degrees
    err(k) = sum( (testtf- ms).^2);
end

%  Find point at which the minimum error occured
[mv, mp] = min(err);  %(mv is the minimum value and mp is the index where the minimumn occurs)
wcest = wc(mp(1));  %  Estimated cut-off frequency
disp([ 'Estimated Cut-off Frequency: ' num2str(wcest/(2*pi)) 'Hz' ])
c = 1/(r*wcest);
disp([ 'Estimated Capacitor Value: ' num2str(c) ])

%  Create a dense frequency axis to plot the parametric TF
%  at estimated cut-off
fplot = [min(f):.5:max(f)];
splot = j*pi*2*fplot;
tf = (splot/wcest) ./ (1 + (splot/wcest));   

%  Plot TF with estimated cut-off and data points and compare
%  to actual data points
figure(1)   %  Bode Plot
semilogx(fplot, 180*(angle(tf))/pi,'b',f,(ms),'xr')
title('Comarison of best fit TF phase to Data Points')
xlabel('Hz')
ylabel('Degrees')
figure(2)  %  % regular plot (no log scale)
plot(fplot, 180*angle(tf)/pi,'b',f,(ms),'xr')
title('Comarison of best fit TF phase to Data Points')
xlabel('Hz')
ylabel('Degrees')
figure(3)  %  Plot error curve
plot(wc/(2*pi),err)
title('Suared Error for each guess of a cut-off frequency')
xlabel('Test Cut-off Frequency in Hz')
ylabel('Squared Error')