%
%
%
%
% Sample rates in excess of 10 times the highest frequency in the signal are typically chosen so that the frequencies of interest are minimally affected by the non-ideal zero order hold. The larger the sampling frequency relative to the frequencies of interest, the less amplitutde will be changed by the zero-order hold.
%
% Fix above. M(t) should blow up immediately
% e(t) has a spectrum amplitude of infinity at zero, but because of sampling effects, this will also be seen at every harmonic of the sampling frequency. The Zero-Order-Hold will then have a much larger signal per sample than the real signal e(t), as e*(t) will have aliasing included. c(t) will be slightly smaller thanks to the low pass filter, but will still capture several multiples of the sampling frequency, and have a much larger signal than e(t)
%
%
%

s = tf('s')
B = 10;
T = 1/1;
G = 1/(1 + (s/2*pi*B));
E = 1/s^2

%% Define the frequency range
omega1 = 1/T;       % Placeholder: maximum frequency (rad/s). Replace as needed.
nPoints = 10000;    % Number of frequency points
omega = linspace(-3*omega1, 3*omega1, nPoints);  % Linearly spaced frequency vector

%% Compute the frequency response
% Evaluate G(j*omega) over the specified range. 
% freqresp returns a multidimensional array; use squeeze to get a vector.
H = squeeze(freqresp(E, 1j*omega));

%% Plot the Magnitude and Phase Spectrum

% Create a figure with two subplots for magnitude and phase
figure;

% Magnitude Plot
subplot(2,1,1);
plot(omega, abs(H), 'b', 'LineWidth', 1.5);
xlabel('\omega (rad/s)');
ylabel('Magnitude');
title('Magnitude Spectrum of G(s)');
grid on;

% Phase Plot
subplot(2,1,2);
plot(omega, angle(H), 'r', 'LineWidth', 1.5);
xlabel('\omega (rad/s)');
ylabel('Phase (radians)');
title('Phase Spectrum of G(s)');
grid on;
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;
