diff --git a/Band Elimination Chebyshev/band_elimination_design.m b/Band Elimination Chebyshev/band_elimination_design.m index b890f68..a8a8542 100644 --- a/Band Elimination Chebyshev/band_elimination_design.m +++ b/Band Elimination Chebyshev/band_elimination_design.m @@ -95,7 +95,7 @@ alpha_parameter = asinh(1/epsilon_parameter)/design_filter_order; % Calculates the frequency at which half power occurs using the eq. 9-80 % TODO: denormalize!! ====================%%%%%%%%%%%%%%%%%%%%%%%%============================ design_half_power_radial_frequency = cosh(acosh(( ... - 10^(specification_max_pass_attenuation/10-1))^(-1/2))/ ... + 10^(specification_max_pass_attenuation/10)-1)^(-1/2))/ ... design_filter_order); % rad/s % ----- @@ -639,25 +639,26 @@ ltiview(high_pass_notch_units_transfer_functions(1,1), ... low_pass_notch_units_transfer_functions(1,2)); %} -% +%{ ltiview(high_pass_notch_units_transfer_functions(1,1), ... high_pass_notch_units_transfer_functions(1,2), ... low_pass_notch_units_transfer_functions(1,1), ... low_pass_notch_units_transfer_functions(1,2), ... total_transfer_function); -% +%} %ltiview(total_transfer_function); -%{ +% plot_transfer_function(total_transfer_function, ... [specification_central_frequency ... - design_half_power_radial_frequency/(2*pi) ... + band_elimination_poles_radial_frequencies(1,1)/(2*pi) ... + design_half_power_radial_frequency*specification_low_pass_radial_frequency/(2*pi) ... specification_low_stop_frequency ... specification_low_pass_frequency ... specification_high_pass_frequency ... specification_high_stop_frequency]); -%} +% % Clears unneeded variable from workspace clearVars = {'total_transfer_function'}; diff --git a/High Pass Butterworth/Multisim/high_pass_butterworth.ms14 b/High Pass Butterworth/Multisim/high_pass_butterworth.ms14 new file mode 100644 index 0000000..e7f3559 Binary files /dev/null and b/High Pass Butterworth/Multisim/high_pass_butterworth.ms14 differ diff --git a/High Pass Butterworth/high_pass_design.m b/High Pass Butterworth/high_pass_design.m new file mode 100644 index 0000000..9d4c892 --- /dev/null +++ b/High Pass Butterworth/high_pass_design.m @@ -0,0 +1,337 @@ +%% AUTHOR : Apostolos Fanakis +%% $DATE : 16-Aug-2018 19:23:40 $ +%% $Revision : 1.00 $ +%% DEVELOPED : 9.0.0.341360 (R2016a) +%% FILENAME : band_pass_design.m +%% AEM : 8261 +%% +%% ========== DESIGN SPECIFICATIONS START ========== +% Figures out design specifications according to my AEM number + +AEM = [8 2 6 1]; + +if AEM(4)<4 + design_param_m = 1; +elseif AEM(3)<7 + design_param_m = 0; +else + design_param_m = 3; +end + +specification_pass_frequency = (4+design_param_m)*1000; % Hz +specification_pass_radial_frequency = ... + specification_pass_frequency*(2*pi); % rad/s + +specification_stop_frequency = specification_pass_frequency/2.6; % Hz +specification_stop_radial_frequency = ... + specification_stop_frequency*(2*pi); % rad/s + +specification_min_stop_attenuation = 24+AEM(4)*(6/9); % dB +specification_max_pass_attenuation = 0.5+AEM(3)/36; % dB + +clear design_param_m + + +%{ +specification_pass_radial_frequency = 15707.96; % rad/s +specification_stop_radial_frequency = 6283.2; % rad/s + +specification_min_stop_attenuation = 25; % dB +specification_max_pass_attenuation = 0.5; % dB +%} +% ========== DESIGN SPECIFICATIONS END ========== + +%% ========== PROTOTYPE LOW PASS DESIGN SPECIFICATIONS START ========== +% Calculates the specifications for the low pass design that will later be +% converted to the desired high pass filter + +% Calculates the specifications using the eq. 12-4 +% prototype_normalized_pass_radial_frequency = 1; % rad/s +prototype_normalized_stop_radial_frequency = ... + specification_pass_radial_frequency/ ... + specification_stop_radial_frequency; % rad/s + +% Min and max attenuations remain the same + +% ========== PROTOTYPE LOW PASS DESIGN SPECIFICATIONS END ========== + +%% ========== PROTOTYPE LOW PASS DESIGN START ========== +% Designs the prototype normalized filter. + +% Calculates the filter's order using the eq. 9-52 +design_filter_order = ceil(log10(((10^ ... + (specification_min_stop_attenuation/10)-1)/(10^ ... + (specification_max_pass_attenuation/10)-1)))/ ... + (2*log10(prototype_normalized_stop_radial_frequency))); + +% Calculates the frequency at which half power of the low pass prototype +% occurs using the eq. 9-48 +low_pass_prototype_half_power_radial_frequency = 1/ ... + (10^(specification_max_pass_attenuation/10)-1)^ ... + (1/(2*design_filter_order)); % rad/s + +% Transforms the result to get the corresponding frequency for the high +% pass using the eq. 12-12 +design_half_power_radial_frequency = specification_pass_radial_frequency/ ... + low_pass_prototype_half_power_radial_frequency; % rad/s + +% ----- +% Calculates stable poles, zeros, angles and other characteristic sizes +% using the Guillemin algorithm for a normalized low pass Butterworth +% filter. +% +% So for the time being we assume that the pass radial frequency is equal +% to unity (1). +% ----- + +% Initializes necessary variables +low_pass_prototype_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 low_pass_prototype_number_of_poles]); +low_pass_prototype_poles_real_parts = ... + zeros([1 low_pass_prototype_number_of_poles]); +low_pass_prototype_poles_imaginary_parts = ... + zeros([1 low_pass_prototype_number_of_poles]); +low_pass_prototype_poles_radial_frequencies = ... + zeros([1 low_pass_prototype_number_of_poles]); +low_pass_prototype_poles_Q = zeros([1 low_pass_prototype_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 +% 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-31 + low_pass_prototype_poles_real_parts(1,1) = ... + -cosd(design_butterworth_angles(1,1)); + % Calculates the first pole's imaginary part using the eq. 9-31 + low_pass_prototype_poles_imaginary_parts(1,1) = ... + sind(design_butterworth_angles(1,1)); + + low_pass_prototype_poles_radial_frequencies(1,1) = 1; + % Calculates the first pole's Q using the eq. 9-38 + low_pass_prototype_poles_Q(1,1) = 1/ ... + (2*cosd(design_butterworth_angles(1,1))); + + % Calculates the rest of the poles in the same way + for i=2:low_pass_prototype_number_of_poles + design_butterworth_angles(1,i) = double((i-1)*theta); + % Pole's real part, eq. 9-31 + low_pass_prototype_poles_real_parts(1,i) = ... + -cosd(design_butterworth_angles(1,i)); + % Pole's imaginary part, eq. 9-31 + low_pass_prototype_poles_imaginary_parts(1,i) = ... + sind(design_butterworth_angles(1,i)); + + low_pass_prototype_poles_radial_frequencies(1,i) = 1; + % Pole's Q, eq. 9-38 + low_pass_prototype_poles_Q(1,i) = 1/ ... + (2*cosd(design_butterworth_angles(1,i))); + end +else % Even number of poles + % Theta is a helper parameter + theta=90/low_pass_prototype_number_of_poles; + + for i=1:low_pass_prototype_number_of_poles + design_butterworth_angles(1,i) = double(90)/ ... + double(design_filter_order)+double((i-1)*theta); + % Pole's real part, eq. 9-31 + low_pass_prototype_poles_real_parts(1,i) = ... + -cosd(design_butterworth_angles(1,i)); + % Pole's imaginary part, eq. 9-31 + low_pass_prototype_poles_imaginary_parts(1,i) = ... + sind(design_butterworth_angles(1,i)); + + low_pass_prototype_poles_radial_frequencies(1,i) = 1; + % Pole's Q, eq. 9-106 + low_pass_prototype_poles_Q(1,i) = 1/ ... + (2*cosd(design_butterworth_angles(1,i))); + end +end + +% Clears unneeded variables from workspace +% +clearVars = {'i', 'prototype_normalized_stop_radial_frequency', ... + 'low_pass_prototype_half_power_radial_frequency', 'theta'}; +clear(clearVars{:}) +clear clearVars +% + +% ========== PROTOTYPE LOW PASS DESIGN END ========== + +%% ========== POLES TRANSFORMATION START ========== +% Transforms the prototype's poles + +% Initializes necessary variables +% Calculates the number of poles that will occur after the transformation +high_pass_number_of_poles = low_pass_prototype_number_of_poles; + +% According to the course notes (chapter 12, end of page 5) the +% transformation leaves the poles unchanged. +high_pass_poles_real_parts = low_pass_prototype_poles_real_parts; +high_pass_poles_imaginary_parts = ... + low_pass_prototype_poles_imaginary_parts; +high_pass_poles_radial_frequencies = ... + low_pass_prototype_poles_radial_frequencies; +high_pass_poles_Q = low_pass_prototype_poles_Q; + +% The transormation also produces a number of zeros at (0,0) equal to the +% filter order +high_pass_transfer_function_zeros = zeros([1 design_filter_order]); + +% Clears unneeded variables from workspace +clear low_pass_prototype_number_of_poles +clear -regexp ^low_pass_prototype_ + +% ========== POLES TRANSFORMATION END ========== + +%% ========== POLES DE-NORMALIZATION START ========== +% The high pass filter poles calculated above are those of a normalized +% filter. De-normalization is needed to get the actual poles for the +% desired high pass filter. + +for i=1:high_pass_number_of_poles + high_pass_poles_real_parts(1,i) = high_pass_poles_real_parts(1,i)* ... + design_half_power_radial_frequency; + high_pass_poles_imaginary_parts(1,i) = ... + high_pass_poles_imaginary_parts(1,i)* ... + design_half_power_radial_frequency; + high_pass_poles_radial_frequencies(1,i) = ... + design_half_power_radial_frequency; +end + +% Clears unneeded variables from workspace +clear i +clear -regexp ^geffe_ +clear -regexp ^transformation_ + +% ========== POLES DE-NORMALIZATION END ========== + +%% ========== UNITS IMPLEMENTATION START ========== +% AEM(2) = 2, so the first design strategy is going to be used in the +% Sallen-Key high pass circuits. + +% ------------------------------------------------------------------------- +% Unit 1 has a pole pair with Q equal to 0.5412 and lies on a circle with a +% radius equal to 25097.78. +% ------------------------------------------------------------------------- +% Unit 1 has a pole pair with Q equal to 1.3066 and lies on a circle with a +% radius equal to 25097.78. +% ------------------------------------------------------------------------- + +% 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. +units_R = ones([1 high_pass_number_of_poles]); +units_C = ones([1 high_pass_number_of_poles]); +units_r1 = ones([1 high_pass_number_of_poles]); +units_r2 = zeros([1 high_pass_number_of_poles]); +units_k = zeros([1 high_pass_number_of_poles]); +units_frequency_scale_factors = zeros([1 2]); +units_amplitude_scale_factors = zeros([1 2]); +units_transfer_functions = [tf(1) tf(1)]; + +for i=1:high_pass_number_of_poles + % Calculates k and r2 using the eq. 6-75 + units_r2(1,i) = 2-1/high_pass_poles_Q(1,i); + units_k(1,i) = 3-1/high_pass_poles_Q(1,i); + + % Selects the appropriate frequency scale factor to transfer the + % normalized radial frequency back to the original + units_frequency_scale_factors(1,i) = ... + high_pass_poles_radial_frequencies(1,i); + % AEM(3) = 6, so the magnitude scaling will be performed to achieve a + % capacitor value of 0.1uF using the eq. 6-33 + units_amplitude_scale_factors(1,i) = ... + units_C(1,i)/ ... + (units_frequency_scale_factors(1,i)*0.1*10^(-6)); + + % Performs scaling + units_R(1,i) = units_R(1,i)* ... + units_amplitude_scale_factors(1,i); + units_C(1,i) = 0.1*10^(-6); + units_r1(1,i) = units_r1(1,i)* ... + units_amplitude_scale_factors(1,i); + units_r2(1,i) = units_r2(1,i)* ... + units_amplitude_scale_factors(1,i); + + % Builds unit's transfer function + % Builds numerator and denominator of the transfer function using the + % eq. 6-68 + G = (units_R(1,i)+units_r2(1,i))/units_R(1,i); + unit_numerator = [G ... + 0 ... + 0]; + unit_denominator = [1 ... + 2/(units_C(1,i)*units_R(1,i))+(1-G)/(units_C(1,i)*units_R(1,i)) ... + 1/(units_C(1,i)^2*units_R(1,i)^2)]; + + units_transfer_functions(1,i) = ... + tf(unit_numerator, unit_denominator); +end + + + +% Clears unneeded variables from workspace +clearVars = {''}; +clear(clearVars{:}) +clear clearVars +clear -regexp _transfer_function$ + +% ========== UNITS IMPLEMENTATION END ========== + +%% ========== GAIN ADJUSTMENT START ========== + +% +total_gain_high = units_k(1,1)*units_k(1,2); +unit_adjustment_gain = 1/total_gain_high; +% We arbitrarily choose to use a 10KOhm series resistor in the adjustment +% unit +unit_adjustment_feedback_resistor = 10*10^3*unit_adjustment_gain; +% + +total_transfer_function = series(units_transfer_functions(1,1), ... + units_transfer_functions(1,2)); +total_transfer_function = total_transfer_function*unit_adjustment_gain; + +%{ +ltiview(units_transfer_functions(1,1), ... + units_transfer_functions(1,2)); +%} + +%{ +ltiview(units_transfer_functions(1,1), ... + units_transfer_functions(1,2), ... + total_transfer_function); +%} + +%ltiview(total_transfer_function); + +% +plot_transfer_function(total_transfer_function, ... + [design_half_power_radial_frequency/(2*pi) ... + specification_stop_frequency ... + specification_pass_frequency ... + 15000]); +% + +% Clears unneeded variable from workspace +clearVars = {'total_transfer_function'}; +clear(clearVars{:}) +clear clearVars +clear -regexp _transfer_functions$ + +% ========== GAIN ADJUSTMENT END ========== \ No newline at end of file diff --git a/High Pass Butterworth/plot_transfer_function.m b/High Pass Butterworth/plot_transfer_function.m new file mode 100644 index 0000000..b56e2d2 --- /dev/null +++ b/High Pass Butterworth/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