From b68a326754bf9765cdd22a97e0dadb6a4aeeabc0 Mon Sep 17 00:00:00 2001 From: Apostolof Date: Sun, 12 Aug 2018 13:51:47 +0300 Subject: [PATCH] Init --- .../low_pass_design_single.m | 397 ++++++++++++++++++ .../plot_transfer_function.m | 44 ++ 2 files changed, 441 insertions(+) create mode 100644 Low Pass Inverse Chebyshev/low_pass_design_single.m create mode 100644 Low Pass Inverse Chebyshev/plot_transfer_function.m diff --git a/Low Pass Inverse Chebyshev/low_pass_design_single.m b/Low Pass Inverse Chebyshev/low_pass_design_single.m new file mode 100644 index 0000000..876047d --- /dev/null +++ b/Low Pass Inverse Chebyshev/low_pass_design_single.m @@ -0,0 +1,397 @@ +% AUTHOR : Apostolos Fanakis +% $DATE : 02-May-2018 17:21:18 $ +% $Revision : 1.00 $ +% DEVELOPED : 9.0.0.341360 (R2016a) +% FILENAME : low_pass_design_single.m +% AEM : 8261 + +% ========== DESIGN SPECIFICATIONS START ========== +% Figures out design specifications according to my AEM number + +AEM = [8 2 6 1]; + +if AEM(3)<3 + design_param_m = 1; +elseif AEM(3)<7 + design_param_m = 2; +else + design_param_m = 3; +end + +specification_pass_frequency = (1.1*(3+design_param_m))*1000; % Hz +specification_pass_radial_frequency = specification_pass_frequency*(2*pi); % rad/s + +specification_stop_frequency = 2.1*specification_pass_frequency; % Hz +specification_stop_radial_frequency = specification_stop_frequency*(2*pi); % rad/s + +specification_min_stop_attenuation = 23+(max(1,AEM(3))-5)*(3/4); % dB +specification_max_pass_attenuation = 0.6+((max(1,AEM(4))-5)/16); % dB + +clear design_param_m + +specification_pass_radial_frequency = 1000; % rad/s +specification_stop_radial_frequency = 1400; % rad/s +specification_min_stop_attenuation = 18; % dB +specification_max_pass_attenuation = 0.25; % dB +% ========== DESIGN SPECIFICATIONS END ========== + +% ========== NORMALIZED DESIGN START ========== +% Calculates normalized design specifications and designs a normalized +% filter using them + +normalized_pass_radial_frequency = specification_pass_radial_frequency/ ... + specification_stop_radial_frequency; % rad/s +% normalized_stop_radial_frequency = 1; % Hz (stop_frequency/stop_frequency) + +% Calculates the filter's order using the eq. 9-137 +design_filter_order = ceil( ... + acosh(((10^(specification_min_stop_attenuation/10)-1)/ ... + (10^(specification_max_pass_attenuation/10)-1))^(1/2)) ... + /acosh(1/normalized_pass_radial_frequency)); +% Calculates epsilon parameter using the eq. 9-123 +epsilon_parameter = 1/(10^(specification_min_stop_attenuation/10)-1)^(1/2); + +% Calculates alpha using the eq. ?? wtf is this? =========^^^^^^^^^^^^^^^^^^&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&========== +alpha_parameter = asinh(1/epsilon_parameter)/design_filter_order; + +% Calculates the frequency at which half power occurs using the eq. 9-139 +design_half_power_radial_frequency = specification_stop_radial_frequency/ ... + (cosh(acosh(1/epsilon_parameter)/design_filter_order)); % rad/s + +% ----- +% Calculates stable poles, zeros, angles and other characteristic sizes +% using the Guillemin algorithm +% ----- + +% Initializes necessary variables +design_number_of_poles = idivide(design_filter_order,int32(2),'ceil'); +% Creates five vector arrays of dimensions [1 * number_of_poles] filled +% with zeros to store: +% - the Butterworth angles with reference to the negative horizontal axes, +% - the real parts of the poles, +% - the imaginary parts of the poles, +% - the radial frequencies (Omega0) of the poles and +% - the Q's of the poles +design_butterworth_angles = zeros([1 design_number_of_poles]); +poles_real_parts = zeros([1 design_number_of_poles]); +poles_imaginary_parts = zeros([1 design_number_of_poles]); +poles_radial_frequencies = zeros([1 design_number_of_poles]); +inverse_poles_Q = zeros([1 design_number_of_poles]); + +% Calculates the Butterworth angles using the method suggested in chapter +% 9 (page 10) of the course notes and then uses them to calculate the +% inverse Chebyshev poles +if mod(design_filter_order,2)~=0 % Odd number of poles + % First pole has a zero angle + design_butterworth_angles(1,1)=0; + % The rest of the poles are scattered in the left half pane with + % equal angles + % Theta is a helper parameter + theta=180/design_filter_order; + + % Calculates the first pole's real part using the eq. 9-102 + poles_real_parts(1,1) = -sinh(alpha_parameter)* ... + cosd(design_butterworth_angles(1,1)); + % Calculates the first pole's imaginary part using the eq. 9-103 + poles_imaginary_parts(1,1) = cosh(alpha_parameter)* ... + sind(design_butterworth_angles(1,1)); + % Calculates the first pole's radial frequency using the eq. 9-150 + poles_radial_frequencies(1,1) = (poles_real_parts(1,1)^2+ ... + poles_imaginary_parts(1,1)^2)^(1/2); + % Calculates the first pole's Q using the eq. 9-151 + inverse_poles_Q(1,1) = 1/ ... + (2*cos(atan(poles_imaginary_parts(1,1)/poles_real_parts(1,1)))); + + % Calculates the rest of the poles in the same way + for i=2:design_number_of_poles + design_butterworth_angles(1,i)=double((i-1)*theta); + % Pole's real part, eq. 9-102 + poles_real_parts(1,i) = -sinh(alpha_parameter)* ... + cosd(design_butterworth_angles(1,i)); + % Pole's imaginary part, eq. 9-103 + poles_imaginary_parts(1,i) = cosh(alpha_parameter)* ... + sind(design_butterworth_angles(1,i)); + % Pole's radial frequency, eq. 9-150 + poles_radial_frequencies(1,i) = (poles_real_parts(1,i)^2+ ... + poles_imaginary_parts(1,i)^2)^(1/2); + % Pole's Q, eq. 9-151 + inverse_poles_Q(1,i) = 1/ ... + (2*cos(atan(poles_imaginary_parts(1,i)/poles_real_parts(1,i)))); + end +else % Even number of poles + % Theta is a helper parameter + theta=90/design_number_of_poles; + + for i=1:design_number_of_poles + design_butterworth_angles(1,i)=double(90)/ ... + double(design_filter_order)+double((i-1)*theta); + % Pole's real part, eq. 9-102 + poles_real_parts(1,i) = -sinh(alpha_parameter)* ... + cosd(design_butterworth_angles(1,i)); + % Pole's imaginary part, eq. 9-103 + poles_imaginary_parts(1,i) = cosh(alpha_parameter)* ... + sind(design_butterworth_angles(1,i)); + % Pole's radial frequency, eq. 9-150 + poles_radial_frequencies(1,i) = (poles_real_parts(1,i)^2+ ... + poles_imaginary_parts(1,i)^2)^(1/2); + % Pole's Q, eq. 9-151 + inverse_poles_Q(1,i) = 1/ ... + (2*cos(atan(poles_imaginary_parts(1,i)/poles_real_parts(1,i)))); + end +end + +% Initializes array to hold the inversed poles +inverse_poles_radial_frequencies = zeros([1 design_number_of_poles]); +% Calculates inverse Chebyshev poles by inversing the Chebyshev poles +% using the eq. 9-146 +for i=1:design_number_of_poles + inverse_poles_radial_frequencies(1,i) = 1/poles_radial_frequencies(1,i); +end + +% Initializes array to hold the transfer function zeros and a temporary +% helper variable +inverse_transfer_function_zeros = zeros([1 length(design_butterworth_angles)]); +temp_index = 1; +% Calculates the transfer function's zeros using the eq. 9-143 +for i=1:2:design_filter_order + inverse_transfer_function_zeros(1,temp_index) = sec((i*pi)/ ... + (2*design_filter_order)); + temp_index = temp_index + 1; +end + +% Clears unneeded variable from workspace +clearVars = {'theta', 'i', 'temp_index', 'alpha_parameter' ... + 'epsilon_parameter', 'normalized_pass_radial_frequency', ... + 'poles_real_parts', 'poles_imaginary_parts', 'poles_radial_frequencies'}; +clear(clearVars{:}) +clear clearVars + +% ========== NORMALIZED DESIGN END ========== + +% ========== ZEROS-POLES GROUPING START ========== + +% Grouping is done "by hand". The first pole, that has a radial frequency +% equal to 0.9631, is grouped with the first zero (1.0824) and the second +% pole, that has a radial frequency equal to 0.7484, is grouped with the +% second zero (2.6131). +% Two low pass notch units are required to implement the desired low pass filter. + +% ========== ZEROS-POLES GROUPING END ========== + +% ========== UNITS IMPLEMENTATION START ========== +% AEM(3) = 6, so the circuit shown in 7.23 is going to be used for the low +% pass notch units. + +% Initializes necessary arrays, each array is 1X2, the first element (1,1) +% corresponds to the first unit and the second element (1,2) to second unit. +unit_low_pass_notch_resistors_1 = zeros([1 2]); +unit_low_pass_notch_resistors_2 = zeros([1 2]); +unit_low_pass_notch_resistors_3 = zeros([1 2]); +unit_low_pass_notch_resistors_4 = zeros([1 2]); +unit_low_pass_notch_resistors_5 = zeros([1 2]); +unit_low_pass_notch_capacitors = zeros([1 2]); +unit_low_pass_notch_gains_high = zeros([1 2]); +unit_low_pass_notch_gains_low = zeros([1 2]); + +% The specifications for the first unit are: +% - inverse pole radial frequency = 0.9631 +% - zero at 1.0824 (zero > inverse pole radial frequency) +% - pole angle = 0.5822 + +% normalized_inverse_pole_radial_frequency = 1; +normalized_transfer_function_zero = inverse_transfer_function_zeros(1,1)/ ... + inverse_poles_radial_frequencies(1,3); + +% According to the design method outlined in 7.6-B, at page 35 +unit_low_pass_notch_resistors_1(1,1) = 1; % Ohm +% Calculates the capacity of the normalized circuit capacitors using the +% eq. 7-150 +unit_low_pass_notch_capacitors(1,1) = 1/(2*inverse_poles_Q(1,3)); % Farad +% Calculates the resistance of R2 using the same equations (7-150) +unit_low_pass_notch_resistors_2(1,1) = 4*inverse_poles_Q(1,3)^2; % Ohm +% Calculates the resistance of R5 using the eq. 7-152 +unit_low_pass_notch_resistors_5(1,1) = (4*inverse_poles_Q(1,3)^2)/ ... + (normalized_transfer_function_zero^2-1); % Ohm +unit_low_pass_notch_resistors_4(1,1) = 1; % Ohm +% Calculates the resistance of R3 using the eq. 7-155 +unit_low_pass_notch_resistors_3(1,1) = (normalized_transfer_function_zero^2)/ ... + (2*inverse_poles_Q(1,3)^2); % Ohm +% Calculates the gain of this unit in high frequencies using the eq. 7-143 +unit_low_pass_notch_gains_high(1,1) = (unit_low_pass_notch_resistors_4(1,1))/ ... + (unit_low_pass_notch_resistors_3(1,1)+unit_low_pass_notch_resistors_4(1,1)); +% Calculates the gain of this unit in low frequencies using the +% eq. 7-146, 7-147, 7-148, setting s = 0 +unit_low_pass_notch_gains_low(1,1) = unit_low_pass_notch_gains_high(1,1)* ... + normalized_transfer_function_zero^2; + +% The specifications for the first unit are: +% - inverse pole radial frequency = $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ +% - zero at $$$$$$$$$$$$$$$$$$$$ (zero > inverse pole radial frequency) +% - pole angle = $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ + +normalized_transfer_function_zero = inverse_transfer_function_zeros(1,2)/ ... + inverse_poles_radial_frequencies(1,2); + +% According to the design method outlined in 7.6-B, at page 35 +unit_low_pass_notch_resistors_1(1,2) = 1; % Ohm +% Calculates the capacity of the normalized circuit capacitors using the +% eq. 7-150 +unit_low_pass_notch_capacitors(1,2) = 1/(2*inverse_poles_Q(1,2)); % Farad +% Calculates the resistance of R2 using the same equations (7-150) +unit_low_pass_notch_resistors_2(1,2) = 4*inverse_poles_Q(1,2)^2; % Ohm +% Calculates the resistance of R5 using the eq. 7-152 +unit_low_pass_notch_resistors_5(1,2) = (4*inverse_poles_Q(1,2)^2)/ ... + (normalized_transfer_function_zero^2-1); % Ohm +unit_low_pass_notch_resistors_4(1,2) = 1; % Ohm +% Calculates the resistance of R3 using the eq. 7-155 +unit_low_pass_notch_resistors_3(1,2) = (normalized_transfer_function_zero^2)/ ... + (2*inverse_poles_Q(1,2)^2); % Ohm +% Calculates the gain of this unit in high frequencies using the eq. 7-143 +unit_low_pass_notch_gains_high(1,2) = (unit_low_pass_notch_resistors_4(1,2))/ ... + (unit_low_pass_notch_resistors_3(1,2)+unit_low_pass_notch_resistors_4(1,2)); +% Calculates the gain of this unit in low frequencies using the +% eq. 7-146, 7-147, 7-148, setting s = 0 +unit_low_pass_notch_gains_low(1,2) = unit_low_pass_notch_gains_high(1,2)* ... + normalized_transfer_function_zero^2; + +% Clears unneeded variable from workspace +clear normalized_transfer_function_zero + +% ========== UNITS IMPLEMENTATION END ========== + +% ========== DENORMALIZATION START ========== + +% Unit sizes rescale +unit_frequency_scale_factors = zeros([1 2]); +unit_amplitude_scale_factors = zeros([1 2]); + +% For unit 1 +unit_frequency_scale_factors(1,1) = specification_stop_radial_frequency* ... + inverse_poles_radial_frequencies(1,3); +% AEM(4) = 1, so the scaling will be performed to achieve a capacitor value +% of 0.1uF using the eq. 6-33 +unit_amplitude_scale_factors(1,1) = unit_low_pass_notch_capacitors(1,1)/ ... + (unit_frequency_scale_factors(1,1)*0.1*10^(-6)); + +unit_low_pass_notch_resistors_1(1,1) = unit_low_pass_notch_resistors_1(1,1)* ... + unit_amplitude_scale_factors(1,1); % Ohm +unit_low_pass_notch_resistors_2(1,1) = unit_low_pass_notch_resistors_2(1,1)* ... + unit_amplitude_scale_factors(1,1); % Ohm +unit_low_pass_notch_resistors_3(1,1) = unit_low_pass_notch_resistors_3(1,1)* ... + unit_amplitude_scale_factors(1,1); % Ohm +unit_low_pass_notch_resistors_4(1,1) = unit_low_pass_notch_resistors_4(1,1)* ... + unit_amplitude_scale_factors(1,1); % Ohm +unit_low_pass_notch_resistors_5(1,1) = unit_low_pass_notch_resistors_5(1,1)* ... + unit_amplitude_scale_factors(1,1); % Ohm +unit_low_pass_notch_capacitors(1,1) = 0.1*10^(-6); % Farad + +% For unit 2 +unit_frequency_scale_factors(1,2) = specification_stop_radial_frequency* ... + inverse_poles_radial_frequencies(1,2); +% AEM(4) = 1, so the scaling will be performed to achieve a capacitor value +% of 0.1uF using the eq. 6-33 +unit_amplitude_scale_factors(1,2) = unit_low_pass_notch_capacitors(1,2)/ ... + (unit_frequency_scale_factors(1,2)*0.1*10^(-6)); + +unit_low_pass_notch_resistors_1(1,2) = unit_low_pass_notch_resistors_1(1,2)* ... + unit_amplitude_scale_factors(1,2); % Ohm +unit_low_pass_notch_resistors_2(1,2) = unit_low_pass_notch_resistors_2(1,2)* ... + unit_amplitude_scale_factors(1,2); % Ohm +unit_low_pass_notch_resistors_3(1,2) = unit_low_pass_notch_resistors_3(1,2)* ... + unit_amplitude_scale_factors(1,2); % Ohm +unit_low_pass_notch_resistors_4(1,2) = unit_low_pass_notch_resistors_4(1,2)* ... + unit_amplitude_scale_factors(1,2); % Ohm +unit_low_pass_notch_resistors_5(1,2) = unit_low_pass_notch_resistors_5(1,2)* ... + unit_amplitude_scale_factors(1,2); % Ohm +unit_low_pass_notch_capacitors(1,2) = 0.1*10^(-6); % Farad + +% ========== DENORMALIZATION END ========== + +% ========== GAIN ADJUSTMENT START ========== + +total_fried_units_attenuation = unit_low_pass_notch_gains_low(1,1)*unit_low_pass_notch_gains_low(1,2); +unit_adjustment_gain = 1/total_fried_units_attenuation; +% We arbitrarily choose to use a 10KOhm series resistor in the adjustment +% unit +unit_adjustment_feedback_resistor = 10*10^3*unit_adjustment_gain; + +% ========== GAIN ADJUSTMENT END ========== + +% ========== TRANSFER FUNCTIONS START ========== + +% Builds numerator and denominator of the transfer function of each unit +% using the eq. 7-146, 7-147 & 7-148 +unit_1_numerator = [1 ... + (((unit_low_pass_notch_gains_high(1,1)-1)/ ... + (unit_low_pass_notch_gains_high(1,1)*unit_low_pass_notch_resistors_1(1,1)* ... + unit_low_pass_notch_capacitors(1,1)))+ ... + (2/(unit_low_pass_notch_resistors_2(1,1)* ... + unit_low_pass_notch_capacitors(1,1)))+ ... + (2/(unit_low_pass_notch_resistors_5(1,1)* ... + unit_low_pass_notch_capacitors(1,1)))) ... + (1/(unit_low_pass_notch_resistors_1(1,1)* ... + unit_low_pass_notch_resistors_5(1,1)*unit_low_pass_notch_capacitors(1,1)^2)+ ... + 1/(unit_low_pass_notch_resistors_1(1,1)* ... + unit_low_pass_notch_resistors_2(1,1)*unit_low_pass_notch_capacitors(1,1)^2))]; +unit_1_denominator = [1 ... + 2/(unit_low_pass_notch_resistors_2(1,1)* ... + unit_low_pass_notch_capacitors(1,1)) ... + 1/(unit_low_pass_notch_resistors_1(1,1)* ... + unit_low_pass_notch_resistors_2(1,1)*unit_low_pass_notch_capacitors(1,1)^2)]; + + +unit_1_transfer_function = tf(unit_1_numerator, unit_1_denominator); +unit_1_transfer_function = unit_1_transfer_function*unit_low_pass_notch_gains_high(1,1); + +unit_2_numerator = [1 ... + ((unit_low_pass_notch_gains_high(1,2)-1)/ ... + (unit_low_pass_notch_gains_high(1,2)*unit_low_pass_notch_resistors_1(1,2)* ... + unit_low_pass_notch_capacitors(1,2))+ ... + 2/(unit_low_pass_notch_resistors_2(1,2)* ... + unit_low_pass_notch_capacitors(1,2))+ ... + 2/(unit_low_pass_notch_resistors_5(1,2)* ... + unit_low_pass_notch_capacitors(1,2))) ... + (1/(unit_low_pass_notch_resistors_1(1,2)* ... + unit_low_pass_notch_resistors_5(1,2)*unit_low_pass_notch_capacitors(1,2)^2)+ ... + 1/(unit_low_pass_notch_resistors_1(1,2)* ... + unit_low_pass_notch_resistors_2(1,2)*unit_low_pass_notch_capacitors(1,2)^2))]; +unit_2_denominator = [1 ... + 2/(unit_low_pass_notch_resistors_2(1,2)* ... + unit_low_pass_notch_capacitors(1,2)) ... + 1/(unit_low_pass_notch_resistors_1(1,2)* ... + unit_low_pass_notch_resistors_2(1,2)*unit_low_pass_notch_capacitors(1,2)^2)]; + +unit_2_transfer_function = tf(unit_2_numerator, unit_2_denominator); +unit_2_transfer_function = unit_2_transfer_function*unit_low_pass_notch_gains_high(1,2); + +tmpHpRC = double(1)/double(0.1*10^(-6)*4.14*10^3); +unit_3_transfer_function = tf([0 tmpHpRC], [1 tmpHpRC]); + +total_transfer_function = series(series(unit_1_transfer_function,unit_2_transfer_function), unit_3_transfer_function); +total_transfer_function = total_transfer_function*unit_adjustment_gain; +%plot_transfer_function(unit_1_transfer_function, [1 10]); +%plot_transfer_function(unit_2_transfer_function, [1 192.82 222.817 234 238.73]); +%plot_transfer_function(0.656*series(unit_1_transfer_function,unit_2_transfer_function), ... +% [1 159.155 192.82 222.817 234 238.73]); +%{ +plot_transfer_function(total_transfer_function, ... + [1 ... + specification_pass_radial_frequency/(2*pi) ... + design_half_power_radial_frequency/(2*pi) ... + specification_stop_radial_frequency/(2*pi)]); +%} + +%ltiview(unit_1_transfer_function); +%ltiview(unit_2_transfer_function); +%ltiview(unit_3_transfer_function); +%ltiview(total_transfer_function); +%ltiview(unit_1_transfer_function, unit_2_transfer_function, ... +% unit_3_transfer_function, total_transfer_function); + +% Clears unneeded variable from workspace +clear tmpHpRC +clear -regexp _numerator$ +clear -regexp _denominator$ +clear -regexp _transfer_function$ + +% ========== TRANSFER FUNCTIONS END ========== \ No newline at end of file diff --git a/Low Pass Inverse Chebyshev/plot_transfer_function.m b/Low Pass Inverse Chebyshev/plot_transfer_function.m new file mode 100644 index 0000000..b56e2d2 --- /dev/null +++ b/Low Pass Inverse Chebyshev/plot_transfer_function.m @@ -0,0 +1,44 @@ +function plot_transfer_function( tf, frequency_markers ) +%PLOT_TRANSFER_FUNCTION Plots bode of a transfer function with markers +% +% tf - The transfer function (created using tf) +% frequency_markers - A matrix of frequencies in Hz +% +% Example: +% plot_transfer_function( tf([1000], [1 1000]), [10 1000 10000] ); + +figure; +x_space = logspace(1,5,1000); % 1000 points between 10^1 and 10^5 +x_space = 2 * pi * x_space; % to rad / sec +[mag,~,wout] = bode(tf,x_space); +mag = squeeze(mag); +wout = squeeze(wout); +mag = 20*log10(mag); +wout = wout/2/pi; +semilogx(wout,mag,'-b'); +axis([min(wout) max(wout) min(mag)-10 max(mag)+10]); +[num,den] = tfdata(tf); +syms s; +d1 = digits(5); +ltx = latex(vpa(poly2sym(cell2mat(num),s)/poly2sym(cell2mat(den),s))); +digits(d1); +title(strcat('$',ltx,'$'), 'Interpreter','latex', 'FontSize', 24); +xlabel('Frequency (Hz)', 'FontSize', 18); +ylabel('Magnitude (dB)', 'FontSize', 18); +grid on; +hold on; +[dbMarks,~,frequency_markers] = bode(tf,2 * pi * frequency_markers); +dbMarks = squeeze(dbMarks); +frequency_markers = squeeze(frequency_markers); +dbMarks = 20*log10(dbMarks); +frequency_markers = frequency_markers/2/pi; +Aw = cell(size(frequency_markers, 1) + 1, 1); +Aw{1} = 'Transfer function'; +for i = 1 : size(frequency_markers, 1) + semilogx(frequency_markers(i),dbMarks(i),'o'); + Aw{i+1} = sprintf('Attenuation at %.2f Hz is %.2f dB', ... + frequency_markers(i), dbMarks(i)); +end +legend(Aw,'Location','best','FontSize',12); +set(gca,'FontSize',14); +end \ No newline at end of file