1 changed files with 369 additions and 0 deletions
@ -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 ========== |
||||
Loading…
Reference in new issue