diff --git a/Band Pass Chebyshev/band_pass_design.m b/Band Pass Chebyshev/band_pass_design.m new file mode 100644 index 0000000..429a684 --- /dev/null +++ b/Band Pass Chebyshev/band_pass_design.m @@ -0,0 +1,388 @@ +%% AUTHOR : Apostolos Fanakis +%% $DATE : 03-Aug-2018 12:32:39 $ +%% $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]; + +% Specification f0 +specification_central_frequency = 0.9*1000; % Hz +specification_central_radial_frequency = 2*pi* ... + specification_central_frequency; % rad/s + +% Specification f1 +specification_low_pass_frequency = 650+25*AEM(3); % Hz +specification_low_pass_radial_frequency = 2*pi* ... + specification_low_pass_frequency; % rad/s + +% Specification f2 +specification_high_pass_frequency = specification_central_frequency^2/ ... + specification_low_pass_frequency; % Hz +specification_high_pass_radial_frequency = 2*pi* ... + specification_high_pass_frequency; % rad/s + +% Specification D +design_param_D = 2.2*((specification_central_frequency^2- ... + specification_low_pass_frequency^2)/specification_low_pass_frequency); + +% Specification f3 +specification_low_stop_frequency = (-design_param_D+ ... + sqrt(design_param_D^2+4*specification_central_frequency^2))/2; % Hz +specification_low_stop_radial_frequency = 2*pi* ... + specification_low_stop_frequency; % rad/s + +% Specification f4 +specification_high_stop_frequency = specification_central_frequency^2/ ... + specification_low_stop_frequency; % Hz +specification_high_stop_radial_frequency = 2*pi* ... + specification_high_stop_frequency; % rad/s + +specification_min_stop_attenuation = 28+AEM(4)*5/9; % dB +specification_max_pass_attenuation = 0.5+AEM(3)/36; % dB + +clear design_param_D + +%{ +specification_low_pass_radial_frequency = 500; % rad/s +specification_high_pass_radial_frequency = 800; % rad/s +specification_low_stop_radial_frequency = 400; % rad/s +specification_high_stop_radial_frequency = 1000; % rad/s +specification_min_stop_attenuation = 18; % 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 band pass filter + +% Calculates the specifications using the eq. 11-56 +% prototype_normalized_pass_radial_frequency = 1; % rad/s +prototype_normalized_stop_radial_frequency = ... + (specification_high_stop_radial_frequency- ... + specification_low_stop_radial_frequency)/ ... + (specification_high_pass_radial_frequency- ... + specification_low_pass_radial_frequency); % rad/s + +% Min and max attenuations remain the same + +% Calculates the geometric middle radial frequency using the eq. 11-2 +design_geometric_central_radial_frequency = ... + sqrt(specification_low_pass_radial_frequency* ... + specification_high_pass_radial_frequency); % rad/s + +% Calculates the pass bandwidth using the eq. 11-52 +design_filter_bandwidth = specification_high_pass_radial_frequency- ... + specification_low_pass_radial_frequency; % rad/s + +% ========== PROTOTYPE LOW PASS DESIGN SPECIFICATIONS END ========== + +%% ========== PROTOTYPE LOW PASS DESIGN START ========== +% The calculated low pass design specifications have a form fit for a low +% pass Chebyshev filter design (the pass radial frequency is normalized to +% one). +% Designs the prototype normalized filter. + +% Calculates the filter's order using the eq. 9-83 +design_filter_order = ceil(acosh(((10^ ... + (specification_min_stop_attenuation/10)-1)/(10^ ... + (specification_max_pass_attenuation/10)-1))^(1/2))/ ... + acosh(prototype_normalized_stop_radial_frequency)); + +% Calculates epsilon parameter using the eq. 9-76 +epsilon_parameter = sqrt(10^(specification_max_pass_attenuation/10)-1); +% Calculates alpha parameter using the eq. 9-92 +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))/design_filter_order); % rad/s + +% ----- +% Calculates stable poles, zeros, angles and other characteristic sizes +% using the Guillemin algorithm +% ----- + +% Initializes necessary variables +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 angles (Qk) of the poles +design_butterworth_angles = zeros([1 prototype_number_of_poles]); +prototype_poles_real_parts = zeros([1 prototype_number_of_poles]); +prototype_poles_imaginary_parts = zeros([1 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-102 + prototype_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 + prototype_poles_imaginary_parts(1,1) = cosh(alpha_parameter)* ... + sind(design_butterworth_angles(1,1)); + + % Calculates the rest of the poles in the same way + for i=2:prototype_number_of_poles + design_butterworth_angles(1,i)=double((i-1)*theta); + % Pole's real part, eq. 9-102 + prototype_poles_real_parts(1,i) = -sinh(alpha_parameter)* ... + cosd(design_butterworth_angles(1,i)); + % Pole's imaginary part, eq. 9-103 + prototype_poles_imaginary_parts(1,i) = cosh(alpha_parameter)* ... + sind(design_butterworth_angles(1,i)); + end +else % Even number of poles + % Theta is a helper parameter + theta=90/prototype_number_of_poles; + + for i=1: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-102 + prototype_poles_real_parts(1,i) = -sinh(alpha_parameter)* ... + cosd(design_butterworth_angles(1,i)); + % Pole's imaginary part, eq. 9-103 + prototype_poles_imaginary_parts(1,i) = cosh(alpha_parameter)* ... + sind(design_butterworth_angles(1,i)); + end +end + +% Clears unneeded variables from workspace +clearVars = {'prototype_normalized_stop_radial_frequency', ... + 'epsilon_parameter', 'alpha_parameter', 'theta'}; +clear(clearVars{:}) +clear clearVars + +% ========== PROTOTYPE LOW PASS DESIGN END ========== + +%% ========== POLES TRANSFORMATION START ========== +% Transforms the prototype's poles according to the Geffe algorithm + +% Initializes necessary variables +% Calculates the parameter qc, required for the transformation of the +% poles, using the eq. 11-6 +transformation_qc_parameter = design_geometric_central_radial_frequency/ ... + design_filter_bandwidth; +% Calculates the number of poles that will occur after the transformation +if mod(design_filter_order,2)~=0 + band_pass_number_of_poles = 2*prototype_number_of_poles-1; +else + band_pass_number_of_poles = 2*prototype_number_of_poles; +end +% Creates four vector arrays of dimensions [1 * 4] filled with zeros to +% store: +% - the Q's of the transformed poles +% - the angles of the transformed poles +% - the radial frequencies of the transformed poles +% - the tranfer function zeros +band_pass_poles_Q = zeros([1 band_pass_number_of_poles]); +band_pass_poles_angle = zeros([1 band_pass_number_of_poles]); +band_pass_poles_radial_frequencies = zeros([1 band_pass_number_of_poles]); +% Every pole transormation produces one transfer function zero at (0,0) for +% every new pole +band_pass_transfer_function_zeros = zeros([1 band_pass_number_of_poles]); +% temp_index is a helper variable +temp_index = 1; + +for i=1:prototype_number_of_poles + if prototype_poles_imaginary_parts(1,i)==0 % Real pole + transformation_sigma_1 = -prototype_poles_real_parts(1,i); + + % Calculates the transformed pole's Q using the eq. 11-11 + band_pass_poles_Q(1,temp_index) = transformation_qc_parameter/ ... + transformation_sigma_1; + % Calculates the transformed pole's angle using the eq. 11-12 + band_pass_poles_angle(1,temp_index) = acosd(1/(2*band_pass_poles_Q(1,i))); + band_pass_poles_radial_frequencies(1,temp_index) = ... + design_geometric_central_radial_frequency; + temp_index = temp_index+1; + else % Complex pole + geffe_sigma_2 = -prototype_poles_real_parts(1,i); + geffe_Omega_2 = prototype_poles_imaginary_parts(1,i); + + % Calculates the parameter C using the eq. 11-28 + geffe_C = geffe_sigma_2^2+geffe_Omega_2^2; + % Calculates the parameter D using the eq. 11-29 + geffe_D = (2*geffe_sigma_2)/transformation_qc_parameter; + % Calculates the parameter E using the eq. 11-30 + geffe_E = 4+geffe_C/transformation_qc_parameter^2; + % Calculates the parameter G using the eq. 11-31 + geffe_G = sqrt(geffe_E^2-4*geffe_D^2); + + % Calculates the parameter Q using the eq. 11-32 + geffe_Q = sqrt((geffe_E+geffe_G)/2)/geffe_D; + % Calculates the parameter k using the eq. 11-33 + geffe_k = (geffe_sigma_2*geffe_Q)/transformation_qc_parameter; + % Calculates the parameter W using the eq. 11-34 + geffe_W = geffe_k+sqrt(geffe_k^2-1); + + % Calculates the radius of the circles upon which the two poles + % reside using the eq. 11-15 + geffe_Omega_0_1 = design_geometric_central_radial_frequency* ... + geffe_W; + geffe_Omega_0_2 = design_geometric_central_radial_frequency/ ... + geffe_W; + + % The two poles have the same Q + band_pass_poles_Q(1,temp_index) = geffe_Q; + band_pass_poles_Q(1,temp_index+1) = geffe_Q; + % Calculates the transformed pole's angle using the eq. 11-37b + band_pass_poles_angle(1,temp_index) = acosd(1/ ... + (2*band_pass_poles_Q(1,temp_index))); + band_pass_poles_angle(1,temp_index+1) = ... + band_pass_poles_angle(1,temp_index); + + band_pass_poles_radial_frequencies(1,temp_index) = geffe_Omega_0_1; + band_pass_poles_radial_frequencies(1,temp_index+1) = geffe_Omega_0_2; + temp_index = temp_index+2; + end +end + +% Clears unneeded variables from workspace +clearVars = {'prototype_number_of_poles', 'i', 'temp_index', ... + 'prototype_poles_imaginary_parts', 'prototype_poles_real_parts'}; +clear(clearVars{:}) +clear clearVars +clear -regexp ^geffe_ +clear -regexp ^transformation_ + +% ========== POLES TRANSFORMATION END ========== + +%% ========== UNITS IMPLEMENTATION START ========== +% AEM(3) = 6, so the first design strategy is going to be used in the +% Delyiannis-Fried circuits. + +% ------------------------------------------------------------------------- +% Unit 1 has a pole pair with ±87.38 degrees of angle, Q equal to 10.954 +% and lies on a circle with a radius equal to 5938.94. +% ------------------------------------------------------------------------- +% Unit 2 has a pole pair with ±87.38 degrees of angle, Q equal to 10.954 +% and lies on a circle with a radius equal to 5384.37. +% ------------------------------------------------------------------------- +% Unit 3 has a pole pair with ±88.92 degrees of angle, Q equal to 26.599 +% and lies on a circle with a radius equal to 6363.12. +% ------------------------------------------------------------------------- +% Unit 4 has a pole pair with ±88.92 degrees of angle, Q equal to 26.599 +% and lies on a circle with a radius equal to 5025.44. +% ------------------------------------------------------------------------- + +% Initializes necessary arrays, each array is 1X4, the first element (1,1) +% corresponds to the first unit and the second element (1,2) to second unit +% etc. +units_BW = zeros([1 band_pass_number_of_poles]); +units_C21 = zeros([1 band_pass_number_of_poles]); +units_C22 = zeros([1 band_pass_number_of_poles]); +units_R21 = zeros([1 band_pass_number_of_poles]); +units_R22 = zeros([1 band_pass_number_of_poles]); +units_frequency_scale_factors = zeros([1 band_pass_number_of_poles]); +units_amplitude_scale_factors = zeros([1 band_pass_number_of_poles]); +units_central_frequency_gain = zeros([1 band_pass_number_of_poles]); +units_alpha = zeros([1 band_pass_number_of_poles]); +units_Z22 = zeros([1 band_pass_number_of_poles]); +units_Z23 = zeros([1 band_pass_number_of_poles]); +unit_transfer_function = [tf(1) tf(1) tf(1) tf(1)]; + +for i=1:band_pass_number_of_poles + % Calculates BW using the eq. 7-62 + units_BW(1,i) = band_pass_poles_radial_frequencies(1,i)/ ... + band_pass_poles_Q(1,i); + % Calculates C21 (=C22=C) using the eq. 7-87 + units_C21(1,i) = 1/(2*band_pass_poles_Q(1,i)); + units_C22(1,i) = units_C21(1,i); + % Using the eq. 7-86 + units_R21(1,i) = 1; + % Calculates R12 using the eq. 7-87 + units_R22(1,i) = 4*band_pass_poles_Q(1,i)^2; + + % Selects the appropriate frequency scale factor to transfer the + % normalized radial frequency back to the original + units_frequency_scale_factors(1,i) = ... + band_pass_poles_radial_frequencies(1,i); + % AEM(4) = 1, 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_C21(1,i)/ ... + (units_frequency_scale_factors(1,i)*0.1*10^(-6)); + + % Scales the circuit elements + units_C21(1,i) = 0.1*10^(-6); + units_C22(1,i) = 0.1*10^(-6); + units_R21(1,i) = units_amplitude_scale_factors(1,i); + units_R22(1,i) = units_R22(1,i)*units_amplitude_scale_factors(1,i); + + % Calculates the gain at the central radial frequency and the alpha + % parameter using the eq. 7-89 + units_central_frequency_gain(1,i) = 2*band_pass_poles_Q(1,i)^2; + units_alpha(1,i) = 1/units_central_frequency_gain(1,i); + + % Calculates the values of the resistors used to diminish the entry + % using the eq. 7-90 (to include the scaling already done the equations + % are used in the form presented at example 7.2) + units_Z22(1,i) = units_R21(1,i)/units_alpha(1,i); + units_Z23(1,i) = units_R21(1,i)/(1-units_alpha(1,i)); + + %TODO: build the tfs + % + unit_numerator = [-1/(units_R21(1,i)*units_C21(1,i)) ... + 0]; + unit_denominator = [1 ... + (2/units_C21(1,i))/units_R22(1,i) ... + 1/(units_R22(1,i)*units_R21(1,i)*units_C21(1,i)^2)]; + % + %{ + unit_numerator = [-2*band_pass_poles_Q(1,i)^2*units_BW(1,i) ... + 0]; + unit_denominator = [1 ... + units_BW(1,i) ... + band_pass_poles_radial_frequencies(1,i)^2]; + %} + unit_transfer_function(i) = tf(unit_numerator, unit_denominator); +end + +total_transfer_function = series(series(series(unit_transfer_function(1), ... + unit_transfer_function(2)), unit_transfer_function(3)), ... + unit_transfer_function(4)); +%total_transfer_function = total_transfer_function*38.2; + +%{ +plot_transfer_function(total_transfer_function, ... + [1 ... + specification_low_stop_frequency ... + specification_low_pass_frequency ... + 900 ... + specification_high_pass_frequency ... + specification_high_stop_frequency]); +%} +%plot_transfer_function(unit_transfer_function(2), [1 10]); +%ltiview(unit_transfer_function(1), unit_transfer_function(2), ... +% unit_transfer_function(3), unit_transfer_function(4)); + +%ltiview(total_transfer_function); + +% Clears unneeded variables from workspace +clearVars = {'units_central_frequency_gain', 'i', 'units_alpha', ... + 'unit_denominator', 'unit_numerator', ... + 'units_amplitude_scale_factors', 'units_frequency_scale_factors'}; +clear(clearVars{:}) +clear clearVars +clear -regexp _transfer_function$ + +% ========== UNITS IMPLEMENTATION END ========== \ No newline at end of file diff --git a/Band Pass Chebyshev/plot_transfer_function.m b/Band Pass Chebyshev/plot_transfer_function.m new file mode 100644 index 0000000..b56e2d2 --- /dev/null +++ b/Band Pass 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