4 changed files with 388 additions and 6 deletions
Binary file not shown.
@ -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 ========== |
||||
@ -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…
Reference in new issue