% 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 phase 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 = 983; % Measured value of resistor m = dlmread('1kout.lvm'); % 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 % 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 testtf = 180*angle(1 ./ (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)); 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, 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')