Browse Source

Band pass init

master
Apostolos Fanakis 6 years ago
parent
commit
899dfa047f
  1. 388
      Band Pass Chebyshev/band_pass_design.m
  2. 44
      Band Pass Chebyshev/plot_transfer_function.m

388
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 ==========

44
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
Loading…
Cancel
Save