% Homework 2 Code - Dane Sabo
% ME 2046 Digital Control Systems
% Problem 4

syms omega T real positive
sympref('FloatingPointOutput',false);
%FOH = ((exp(1j*omega*T/2) - exp(-1j*omega*T/2))/(1j*omega*T))^2
FOH = ((1-exp(-T*1j*omega))/(1j*omega))^2 * (T*1j*omega + 1)/T
latex(FOH)
magnitude_FOH = sqrt(real(FOH)^2 + imag(FOH)^2)
magnitude_FOH = simplify(magnitude_FOH)
phase_FOH = simplify(atan2(imag(FOH),real(FOH)))
figure;
% Example parameters
T = 1;                          % A sample scalar value for T
omega_range = linspace(0.1, 10, 100000);  % Frequency array (avoiding 0 to prevent division-by-zero)

% Compute the magnitude element-wise:
%   sqrt(T.^2 .* omega_range.^2 + 1): element-wise operations for square and square root
%   cos(T .* omega_range): cosine computed for each element
%   Multiplications (.*) and divisions (./) are used to enforce element-wise operations.
magnitude_range = ((-2 .* sqrt(T.^2 .* omega_range.^2 + 1) .* cos(T .* omega_range) - 1) ./ T) ./ (omega_range.^2);

% Find the maximum magnitude and its corresponding index:
[maxMagnitude, idx] = max(magnitude_range);
% Find the frequency corresponding to the maximum magnitude:
maxFrequency = omega_range(idx);

% Print the maximum magnitude and frequency:
fprintf('The maximum magnitude is %f at frequency %f T.\n', maxMagnitude, maxFrequency*T);

% Optional: Plotting the magnitude vs frequency for visualization
figure;
plot(omega_range, magnitude_range, 'LineWidth', 2);
xlabel('\omega');           % Label for the x-axis
ylabel('Magnitude');        % Label for the y-axis
title('Magnitude Plot');    % Plot title
grid on;                    % Enable grid for better visualization
% Define parameters
T = 1;                           % Sample period (seconds)
omega_range = linspace(1e-3, 100*T, 10000); % Frequency vector (rad/s)
                                   % Starting at 1e-3 to avoid division by zero

% Compute the Zero-Order Hold (ZOH) frequency response:
% H_zoh(jω) = (1 - exp(-jωT)) / (jω)
% '1j' in MATLAB represents the imaginary unit.
H_zoh = (1 - exp(-1j * T * omega_range)) ./ (1j * omega_range);

% Compute the First-Order Hold (FOH) frequency response:
% H_foh(jω) = ((1 - exp(-T*1j*omega))/(1j*omega))^2 * ((T*1j*omega + 1)/T)
% Use element-wise operations for arrays:
H_foh = (((1 - exp(-T * 1j * omega_range)) ./ (1j * omega_range)).^2) .* ((T * 1j * omega_range + 1) / T);
% Compute magnitude and phase for each hold
mag_zoh = abs(H_zoh);
phase_zoh = angle(H_zoh);
mag_foh = abs(H_foh);
phase_foh = angle(H_foh);

% Plot the magnitude responses for both ZOH and FOH
figure;
plot(omega_range*T, mag_zoh, 'b-', 'LineWidth', 2);  % Blue solid line for ZOH
hold on;
plot(omega_range*T, mag_foh, 'r-', 'LineWidth', 2);  % Red solid line for FOH
xlabel('\omega (rad/s) \times T');
ylabel('Magnitude');
title('Magnitude Response of ZOH and FOH');
legend('ZOH', 'FOH');
grid on;
hold off;
% Plot the phase responses for both ZOH and FOH
figure;
plot(omega_range*T, phase_zoh, 'b-', 'LineWidth', 2);  % Blue solid line for ZOH
hold on;
plot(omega_range*T, phase_foh, 'r-', 'LineWidth', 2);  % Red solid line for FOH
xlabel('\omega (rad/s) \times T');
ylabel('Phase (radians)');
title('Phase Response of ZOH and FOH');
legend('ZOH', 'FOH');
grid on;
hold off;
%This follows some written work I was doing, but is incomplete.
% I may return to this later, but for your times sake I wouldn't bother reading.
syms s
syms T real positive
ZOH = (1-exp(-s*T))/s
FOH = (1-exp(-T*s))^2 * (T*s + 1)/(T*s^2)
G = (T*s+1)/(T*s) * (T*s*ZOH)/(T*s + ZOH)
collect(expand(G))
collect(expand(FOH))