% EE3032 W6D3 Example, Winter 2019-'20, Dr. Durant
% Apply an averaging/integrating/rectangular impulse response to a periodic
% waveform using the transfer function and Fourier Series coefficients.
% Bug from class; now fixed: in the expressions for part, "2*pi" should have
% been used to convert from hertz to radians/second, but the "*pi" was missing.
% ***** Try different values of T, T0, n, etc.! (Adjust t and the first omega if needed.)
close all
T = 0.372; % s, duration of the averaging system
omega = -40:0.01:40; % omegas for plotting the system transfer function
Hmag = T * sinc((omega*T/2)/pi); % on the board, we separated magnitude...
Hang = -T*omega/2; % and phase components of H(omega).
% The "/pi" above converts from the engineer's definition of sinc to MATLAB's definition.
figure
subplot(2,1,1), plot(omega,Hmag), ylabel('Gain (V/V)'), title('H(\omega)')
subplot(2,1,2), plot(omega,Hang), ylabel('Phase shift (radians)')
xlabel('Frequency \omega (rad/s)')
% Now, calculate the Fourier Series coefficients for the even, 0 V average,
% triangle wave. These will later be put through the transfer function
% H(omega) to compute the output. These are based on Table 5-4 from page 210.
n = 1:2:5; % :2: means step by 2, so only use odd numbers; the table told us only the odd components are present in this signal
an = 8 ./ (n.^2*pi^2); % formula from the table
T0 = 1; % period of the input signal
t = -2:0.01:2; % time for plotting input and input
x = zeros(size(t));
for idx=1:length(n) % foreach n
part = an(idx) * cos(2*pi*n(idx)*t/T0); % FS coefficient * basis function/building block
x = x + part; % summation per the FS formula
end
%figure,plot(t,x) % plot just the input
omega = n * (1/T0) * 2*pi; % omegas (n * (f * 2pi)) for the FS frequencies of the input
Hmag = T * sinc((omega*T/2)/pi); % Same formulas as before; apply to the new omegas
Hang = -T*omega/2;
y = zeros(size(t));
for idx=1:length(n)
% same as the loop for x, except we've added gain Hmag and phase shift Hang
part = an(idx) * Hmag(idx) * ...
cos(2*pi*n(idx)*t/T0 + Hang(idx));
y = y + part;
end
figure,plot(t,x, t,y),legend('x(t)','y(t)')
xlabel('Time (s)'), ylabel('Voltage (V)')