From db79e41b301fbd7102ff2fdd0051f3e549bc9645 Mon Sep 17 00:00:00 2001 From: Apostolof Date: Tue, 14 Aug 2018 18:08:05 +0300 Subject: [PATCH] Band elimination filter init retry --- .../band_elimination_design.m | 369 ++++++++++++++++++ 1 file changed, 369 insertions(+) diff --git a/Band Elimination Chebyshev/band_elimination_design.m b/Band Elimination Chebyshev/band_elimination_design.m index e69de29..f814b93 100644 --- a/Band Elimination Chebyshev/band_elimination_design.m +++ b/Band Elimination Chebyshev/band_elimination_design.m @@ -0,0 +1,369 @@ +%% AUTHOR : Apostolos Fanakis +%% $DATE : 14-Aug-2018 15:32:12 $ +%% $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 = 2.5*1000; % Hz +specification_central_radial_frequency = 2*pi* ... + specification_central_frequency; % rad/s + +% Specification f1 +specification_low_pass_frequency = 1650+50*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 = (specification_central_frequency^2- ... + specification_low_pass_frequency^2)/ ... + (specification_low_pass_frequency*2.5); + +% 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 = 26+AEM(3)*5/9; % dB +specification_max_pass_attenuation = 0.5+AEM(4)/18; % dB + +clear design_param_D + +%{ +specification_low_pass_radial_frequency = 1000; % rad/s +specification_high_pass_radial_frequency = 3000; % rad/s +specification_low_stop_radial_frequency = 1400; % rad/s +specification_high_stop_radial_frequency = 2142; % rad/s +specification_min_stop_attenuation = 15; % dB +specification_max_pass_attenuation = 1; % dB +specification_central_radial_frequency = sqrt( ... + specification_low_pass_radial_frequency* ... + specification_high_pass_radial_frequency); +%} +% ========== DESIGN SPECIFICATIONS END ========== + +%% ========== PROTOTYPE LOW PASS DESIGN SPECIFICATIONS START ========== +% Calculates the specifications for the low pass design that will later be +% converted to a high pass filter + +% Calculates the specifications using the eq. 13-9 +% prototype_normalized_pass_radial_frequency = 1; % rad/s +prototype_normalized_stop_radial_frequency = ... + (specification_high_pass_radial_frequency- ... + specification_low_pass_radial_frequency)/ ... + (specification_high_stop_radial_frequency- ... + specification_low_stop_radial_frequency); % rad/s + +% Min and max attenuations remain the same + +% Calculates the geometric middle radial frequency using the eq. 13-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. 13-1 +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 +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-102 + low_pass_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 + low_pass_prototype_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-106 + low_pass_prototype_poles_radial_frequencies(1,1) = sqrt( ... + low_pass_prototype_poles_real_parts(1,1)^2+ ... + low_pass_prototype_poles_imaginary_parts(1,1)^2); + % Calculates the first pole's Q using the same equations (9-106) + low_pass_prototype_poles_Q(1,1) = ... + low_pass_prototype_poles_radial_frequencies(1,1)/ ... + (2*low_pass_prototype_poles_real_parts(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-102 + low_pass_prototype_poles_real_parts(1,i) = ... + -sinh(alpha_parameter)*cosd(design_butterworth_angles(1,i)); + % Pole's imaginary part, eq. 9-103 + low_pass_prototype_poles_imaginary_parts(1,i) = ... + cosh(alpha_parameter)*sind(design_butterworth_angles(1,i)); + % Pole's radial frequency, eq. 9-106 + low_pass_prototype_poles_radial_frequencies(1,i) = sqrt( ... + low_pass_prototype_poles_real_parts(1,i)^2+ ... + low_pass_prototype_poles_imaginary_parts(1,i)^2); + % Pole's Q, eq. 9-106 + low_pass_prototype_poles_Q(1,i) = ... + low_pass_prototype_poles_radial_frequencies(1,i)/ ... + (2*abs(low_pass_prototype_poles_real_parts(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-102 + low_pass_prototype_poles_real_parts(1,i) = ... + -sinh(alpha_parameter)*cosd(design_butterworth_angles(1,i)); + % Pole's imaginary part, eq. 9-103 + low_pass_prototype_poles_imaginary_parts(1,i) = ... + cosh(alpha_parameter)*sind(design_butterworth_angles(1,i)); + % Pole's radial frequency, eq. 9-106 + low_pass_prototype_poles_radial_frequencies(1,i) = sqrt( ... + low_pass_prototype_poles_real_parts(1,i)^2+ ... + low_pass_prototype_poles_imaginary_parts(1,i)^2); + % Pole's Q, eq. 9-106 + low_pass_prototype_poles_Q(1,i) = ... + low_pass_prototype_poles_radial_frequencies(1,i)/ ... + (2*abs(low_pass_prototype_poles_real_parts(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 ========== + +%% ========== LOW PASS TO HIGH PASS TRANSFORMATION START ========== +% Inverses the poles of the low pass prototype to transform it to a high +% pass filter that will later be converted to the desired band +% elimination filter + +% Initializes necessary variables +high_pass_prototype_number_of_poles = low_pass_prototype_number_of_poles; +% Creates five vector arrays of dimensions [1 * number_of_poles] filled +% with zeros to store: +% - the radial frequencies (Omega0) of the poles and +% - the Q's of the poles +% - the angles of the poles +% - the real parts of the poles, +% - the imaginary parts of the poles, +high_pass_prototype_poles_radial_frequencies = ... + zeros([1 high_pass_prototype_number_of_poles]); +high_pass_prototype_poles_Q = ... + zeros([1 high_pass_prototype_number_of_poles]); +high_pass_prototype_poles_angles = ... + zeros([1 high_pass_prototype_number_of_poles]); +high_pass_prototype_poles_real_parts = ... + zeros([1 high_pass_prototype_number_of_poles]); +high_pass_prototype_poles_imaginary_parts = ... + zeros([1 high_pass_prototype_number_of_poles]); + +for i=1:high_pass_prototype_number_of_poles + % Calculates the inversed pole's radial frequencies using the eq. 13-14 + high_pass_prototype_poles_radial_frequencies(1,i) = 1/ ... + low_pass_prototype_poles_radial_frequencies(1,i); + % The Q of the poles remain the same + high_pass_prototype_poles_Q(1,i) = low_pass_prototype_poles_Q(1,i); + % Calculates the inversed pole's angle using the eq. 9-113 + high_pass_prototype_poles_angles(1,i) = acosd(1/ ... + (2*high_pass_prototype_poles_Q(1,i))); + % Calculates the real and imaginary parts of the inverted poles + high_pass_prototype_poles_real_parts(1,i) = ... + -high_pass_prototype_poles_radial_frequencies(1,i)* ... + cosd(high_pass_prototype_poles_angles(1,i)); + high_pass_prototype_poles_imaginary_parts(1,i) = ... + high_pass_prototype_poles_radial_frequencies(1,i)* ... + sind(high_pass_prototype_poles_angles(1,i)); +end + +% Clears unneeded variables from workspace +% +clearVars = {'i', 'high_pass_prototype_poles_radial_frequencies', ... + 'high_pass_prototype_poles_Q', 'high_pass_prototype_poles_angles'}; +clear(clearVars{:}) +clear clearVars +% +clear -regexp ^low_pass_prototype_ + +% ========== LOW PASS TO HIGH PASS TRANSFORMATION 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_elimination_number_of_poles = 2*high_pass_prototype_number_of_poles-1; +else + band_elimination_number_of_poles = 2*high_pass_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 transfer function zeros +band_elimination_poles_Q = zeros([1 band_elimination_number_of_poles]); +band_elimination_poles_angle = zeros([1 band_elimination_number_of_poles]); +band_elimination_poles_radial_frequencies = zeros( ... + [1 band_elimination_number_of_poles]); +% Every pole transformation produces one transfer function zero at (0,0) for +% every new pole +band_elimination_transfer_function_zeros = zeros( ... + [1 2*design_filter_order]); +% temp_index is a helper variable +temp_index = 1; + +% The transformation from high pass to band elimination produces pairs of +% imaginary zeros. The number of the pairs is equal to the filter order. +for i=1:2:(2*design_filter_order) + band_elimination_transfer_function_zeros(1,i) = 1i* ... + design_geometric_central_radial_frequency; + band_elimination_transfer_function_zeros(1,i+1) = -1i* ... + design_geometric_central_radial_frequency; +end + +for i=1:high_pass_prototype_number_of_poles + if high_pass_prototype_poles_imaginary_parts(1,i)==0 % Real pole + transformation_sigma_1 = -high_pass_prototype_poles_real_parts(1,i); + + % Calculates the transformed pole's Q using the eq. 11-11 + band_elimination_poles_Q(1,temp_index) = ... + transformation_qc_parameter/transformation_sigma_1; + % Calculates the transformed pole's angle using the eq. 11-12 + band_elimination_poles_angle(1,temp_index) = acosd(1/ ... + (2*band_elimination_poles_Q(1,i))); + band_elimination_poles_radial_frequencies(1,temp_index) = ... + design_geometric_central_radial_frequency; + temp_index = temp_index+1; + else % Complex pole + geffe_sigma_2 = -high_pass_prototype_poles_real_parts(1,i); + geffe_Omega_2 = high_pass_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_elimination_poles_Q(1,temp_index) = geffe_Q; + band_elimination_poles_Q(1,temp_index+1) = geffe_Q; + % Calculates the transformed pole's angle using the eq. 11-37b + band_elimination_poles_angle(1,temp_index) = acosd(1/ ... + (2*band_elimination_poles_Q(1,temp_index))); + band_elimination_poles_angle(1,temp_index+1) = ... + band_elimination_poles_angle(1,temp_index); + + band_elimination_poles_radial_frequencies(1,temp_index) = ... + geffe_Omega_0_1; + band_elimination_poles_radial_frequencies(1,temp_index+1) = ... + geffe_Omega_0_2; + temp_index = temp_index+2; + end +end + +% Clears unneeded variables from workspace +clearVars = {'i', 'temp_index'}; +clear(clearVars{:}) +clear clearVars +clear -regexp ^high_pass_prototype_ +clear -regexp ^geffe_ +clear -regexp ^transformation_ +% ========== POLES TRANSFORMATION END ========== \ No newline at end of file