% This script will load a sample file obtained from a frequency sweep
% Labview program of a low pass filter and estimate the cut-off from
% all available data using the magnitude of the transfer function:
%
% vout 1
% ------ = -----------------------
% 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 = 1940; % Measured value of resistor
m = dlmread('1960.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*[1000:1:15000]; % In radians per second
% The measured transfer function amplitudes
ms = (m(:,3) ./ m(:,2)); % Ouput over input
% 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 = abs((s/wc(k)) ./ (1+(s/wc(k))));
testtf = abs(1./ (1+(s/wc(k))));
err(k) = sum( (testtf- ms).^2); % Compute sum of all deviations for total error
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)));
tf = 1 ./ (1 + (splot/(wcest)));
% Plot TF with estimated cut-off and data points and compare
% to actual data points
figure(1) % Bode Plot
semilogx(fplot, 20*log10(abs(tf)),'b',f,20*log10((ms)),'xr')
title('Comarison of best fit TF Magnitude to Data Points')
xlabel('Hz')
ylabel('dB')
figure(2) % % regular plot (no log scale)
plot(fplot, abs(tf),'b',f,(ms),'xr')
title('Comarison of best fit TF Magnitude to Data Points')
xlabel('Hz')
ylabel('Gain')
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')