Low pass and band pass response to input signal Matlab and Multisim

Band Pass Chebyshev/band_pass_design.m

@@ -480,28 +480,20 @@ clear clearVars
 
 % ========== UNITS IMPLEMENTATION END ==========
 
-%% ========== TRANSFER FUNCTIONS START ==========
-
-%{
-total_gain = units_filter_central_frequency_gain(1,1)* ...
-    units_filter_central_frequency_gain(1,2)* ...
-    units_filter_central_frequency_gain(1,3)* ...
-    units_filter_central_frequency_gain(1,4);
-%}
+%% ========== TRANSFER FUNCTIONS STUDY START ==========
 
 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*(1/total_gain);
-
-%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), total_transfer_function);
-
-%ltiview(total_transfer_function);
-
 %{
+% Plots each unit's frequency response
+for i=1:band_pass_number_of_poles
+    plot_transfer_function(unit_transfer_function(1,i), ...
+        design_geometric_central_radial_frequency);
+end
+
+% Plots the total filter frequency response
 plot_transfer_function(total_transfer_function, ...
     [specification_low_stop_frequency ...
     specification_low_pass_frequency ...
@@ -509,9 +501,90 @@ plot_transfer_function(total_transfer_function, ...
     941.83 ...
     specification_high_pass_frequency ...
     specification_high_stop_frequency]);
+
+% Plots the total filter attenuation function
+plot_transfer_function(inv(total_transfer_function), ...
+    [specification_low_stop_frequency ...
+    specification_low_pass_frequency ...
+    specification_central_frequency ...
+    941.83 ...
+    specification_high_pass_frequency ...
+    specification_high_stop_frequency]);
 %}
+%{
+ltiview(unit_transfer_function(1,1));
+ltiview(unit_transfer_function(1,2));
+ltiview(unit_transfer_function(1,3));
+ltiview(unit_transfer_function(1,4));
+ltiview(total_transfer_function);
+ltiview(unit_transfer_function(1,1), unit_transfer_function(1,2), ...
+    unit_transfer_function(1,3), unit_transfer_function(1,4), ...
+    total_transfer_function);
+%}
+
+hold off
+
+sampling_time_seconds = 60; % s
+sampling_frequency_Fs = 80000; % Hz
+sampling_period_T = 1/sampling_frequency_Fs; % s
+sampling_time_vector = 0:sampling_period_T: ...
+    sampling_time_seconds-sampling_period_T;
+sampling_length_L = length(sampling_time_vector);
+frequency_vector = sampling_frequency_Fs/sampling_length_L* ...
+    (0:(sampling_length_L/2));
+
+input_cos_signal_frequency_1 = specification_central_radial_frequency- ...
+    (specification_central_radial_frequency- ...
+    specification_low_pass_radial_frequency)/2;
+input_cos_signal_frequency_2 = specification_central_radial_frequency+ ...
+    (specification_central_radial_frequency+ ...
+    specification_low_pass_radial_frequency)/3;
+input_cos_signal_frequency_3 = 0.4*specification_low_stop_radial_frequency;
+input_cos_signal_frequency_4 = 2.5*specification_high_stop_radial_frequency;
+input_cos_signal_frequency_5 = 3*specification_high_stop_radial_frequency;
+
+input_signal = cos(input_cos_signal_frequency_1*sampling_time_vector)+ ...
+    0.8*cos(input_cos_signal_frequency_2*sampling_time_vector)+ ...
+    0.8*cos(input_cos_signal_frequency_3*sampling_time_vector)+ ...
+    0.6*cos(input_cos_signal_frequency_4*sampling_time_vector)+ ...
+    0.5*cos(input_cos_signal_frequency_5*sampling_time_vector);
+
+[gain,phase] = ...
+        bode(total_transfer_function, ...
+        850*2*pi);
+
+system_output = lsim(total_transfer_function, input_signal, ...
+    sampling_time_vector);
+
+% Plots only the first 2000 samples
+figure(1)
+plot(sampling_time_vector(1:2000), input_signal(1:2000), ...
+    sampling_time_vector(1:2000), system_output(1:2000));
+grid on
+
+
+% Plots the power spectrum of the input signal
+input_square_wave_fft = fft(input_signal);
+Pyy = input_square_wave_fft.*conj(input_square_wave_fft)/sampling_length_L;
+figure(2)
+semilogx(frequency_vector,Pyy(1:sampling_length_L/2+1))
+grid on
+
+
+% Plots the power spectrum of the output signal
+system_output_fft = fft(system_output, length(sampling_time_vector));
+Pyy = system_output_fft.*conj(system_output_fft)/sampling_length_L;
+figure(3)
+semilogx(frequency_vector,Pyy(1:sampling_length_L/2+1))
+grid on
 
 % Clears unneeded variables from workspace
+clearVars = {'temp', 'Pyy', 'frequency_vector', ...
+    'system_output', 'system_output_fft'};
+clear(clearVars{:})
+clear clearVars
 clear -regexp _transfer_function$
+clear -regexp ^sampling_
+%clear -regexp ^input_
 
-% ========== TRANSFER FUNCTIONS END ==========
+% ========== TRANSFER FUNCTIONS STUDY END ==========

Low Pass Inverse Chebyshev/low_pass_design.m

@@ -59,7 +59,7 @@ design_filter_order = ceil(temp_filter_order);
 % Calculates epsilon parameter using the eq. 9-123
 epsilon_parameter = 1/(10^(specification_min_stop_attenuation/10)-1)^(1/2);
 
-% Calculates alpha using the eq. ?? wtf is this?  =========^^^^^^^^^^^^^^^^^^&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&==========
+% Calculates alpha 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-139
@@ -91,12 +91,12 @@ inverse_chebyshev_poles_Q = zeros([1 design_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
-% inverse Chebyshev poles
+% 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
+    % The rest of the poles are scattered in the left half pane in equal
+    % angle arcs
     % Theta is a helper parameter
     theta=180/design_filter_order;
     
@@ -281,6 +281,7 @@ for i=1:design_number_of_poles
     unit_low_pass_notch_resistors_5(1,i) = ...
         (4*inverse_chebyshev_poles_Q(1,i)^2)/ ...
         (normalized_transfer_function_zero^2-1); % Ohm
+    
     unit_low_pass_notch_resistors_4(1,i) = 1; % Ohm
     % Calculates the resistance of R3 using the eq. 7-155
     unit_low_pass_notch_resistors_3(1,i) = ...
@@ -449,51 +450,72 @@ ltiview(unit_transfer_function(1,1), unit_transfer_function(1,2), ...
 
 hold off
 
-Fs = 1000;            % Sampling frequency                    
-T = 1/Fs;             % Sampling period       
-L = 1500;             % Length of signal
-t = (0:L-1)*T;        % Time vector
-
 input_signal_frequency = 2000; % Hz
-input_signal_number_of_periods = 100000;
-input_signal_T = input_signal_number_of_periods*(1/input_signal_frequency);
 input_signal_duty_cycle = 20;
-input_signal_Fs = 80000;
-input_signal_dt = 1/input_signal_Fs;
-input_signal_t = 0:input_signal_dt:input_signal_T-input_signal_dt;
 
-input_square_wave = (1+square(2*pi*input_signal_frequency*input_signal_t, ...
+number_of_periods_sampled = 120000;
+sampling_time_seconds = number_of_periods_sampled* ...
+    (1/input_signal_frequency); % s
+
+sampling_frequency_Fs = 80000; % Hz
+sampling_period_T = 1/sampling_frequency_Fs; % s
+sampling_time_vector = 0:sampling_period_T: ...
+    sampling_time_seconds-sampling_period_T;
+sampling_length_L = length(sampling_time_vector);
+frequency_vector = sampling_frequency_Fs/sampling_length_L* ...
+    (0:(sampling_length_L/2));
+
+input_square_wave = (1+ ...
+    square(2*pi*input_signal_frequency*sampling_time_vector, ...
     input_signal_duty_cycle))/2;
 
-system_output = lsim(total_transfer_function, input_square_wave, input_signal_t);
+system_output = lsim(total_transfer_function, input_square_wave, ...
+    sampling_time_vector);
+
 % Plots only the first 10 periods
-temp = 10*input_signal_Fs/input_signal_frequency;
+temp = 10*sampling_frequency_Fs/input_signal_frequency;
 figure(1)
-plot(input_signal_t(1:temp), input_square_wave(1:temp), ...
-    input_signal_t(1:temp), system_output(1:temp));
+plot(sampling_time_vector(1:temp), input_square_wave(1:temp), ...
+    sampling_time_vector(1:temp), system_output(1:temp));
 grid on
 
+% Plots the power spectrum of the input signal
 input_square_wave_fft = fft(input_square_wave);
-L = length(input_signal_t);
-P2 = abs(input_square_wave_fft/L);
-P1 = P2(1:L/2+1);
-P1(2:end-1) = 2*P1(2:end-1);
-f = input_signal_Fs*(0:(L/2))/L;
+Pyy = input_square_wave_fft.*conj(input_square_wave_fft)/sampling_length_L;
 figure(2)
-%plot(f,P1)
-semilogx(f,P1)
+semilogx(frequency_vector,Pyy(1:sampling_length_L/2+1))
 grid on
 
-system_output_fft = fft(system_output, length(input_signal_t));
-P2 = abs(system_output_fft/L);
-P1 = P2(1:L/2+1);
-P1(2:end-1) = 2*P1(2:end-1);
+% Plots the power spectrum of the output signal
+system_output_fft = fft(system_output, length(sampling_time_vector));
+Pyy = system_output_fft.*conj(system_output_fft)/sampling_length_L;
 figure(3)
-%plot(f, P11)
-semilogx(f,P1)
+semilogx(frequency_vector,Pyy(1:sampling_length_L/2+1))
 grid on
 
-%filter_fft = fft(total_transfer_function);
+%{
+asdf = 1:length(frequency_vector);
+for i=1:length(frequency_vector)
+    if i<175000
+        asdf(i) = 1;
+    elseif mod(i,200)==0
+        [asdf(i),phase] = ...
+        bode(total_transfer_function, ...
+        frequency_vector(i)*2*pi);
+    else
+        asdf(i) = 0;
+    end
+    %{
+    if mod(i,10000)==0
+        disp(i)
+    end
+    %}
+end
+
+figure(4)
+semilogx(frequency_vector,asdf,'.',frequency_vector,P1,'-')
+grid on
+%}
 
 % Clears unneeded variable from workspace
 clear low_frequency