% This file demonstrates different implementations of numerical
% convolutions of two functions, as well as how to convert between time and
% frequency domain using the Fast Fourier Transform (FFT). You are
% encouraged to try out the different parts of the file, changing the
% parameters and see how it affects the outcome. Understanding these
% operations makes it easier to write your own codes.
%
% Daniel Sjoberg
% 2016-08-30
% We first demonstrate how two general functions f1 and f2 stored in row
% vectors can be convolved with each other, creating a new function f3 of
% the same length. This can be used for filtering, where the filter has
% bounded support (the function f2 below). In the case below, the resulting
% function f3 is smooth, that is, the fast oscillations in f1 are averaged
% out when f1 is convolved with the filter f2.
Nt = 1000;
t = linspace(0, 10, Nt);
dt = t(2) - t(1);
f1 = sin(t) + 0.1*sin(2*pi*10*t); % The sum of a slow oscillation and a fast oscillation
f2 = exp(-(t-5).^2/(2*0.1^2))/(0.1*sqrt(2*pi)); % Unit area Gaussian pulse
f3 = dt*conv(f1, f2, 'same'); % The keyword 'same' makes the result the same length as the input, keeping the central part
figure()
plot(t, f1, t, f2, t, f3);
xlabel('t')
legend(['f1'; 'f2'; 'f3'])
grid on
% Then we show how two causal signals can be convolved, producing another
% causal signal. A causal signal is zero until a time t=0 (corresponding
% to index 1).
alpha = 1.5;
tau = 1.5;
t0 = 4.80;
omega_a = 1;
Nt = 1000;
t = linspace(0, 30, Nt);
dt = t(2) - t(1);
f1 = alpha*exp(-t/tau);
f2 = exp(-t/t0).*sin(omega_a*t);
f3 = dt*conv(f1, f2);
f3 = f3(1:Nt); % Strip off the remainder of the convolution, making
% the output the same length as f1 and f2.
figure()
plot(t, f1, t, f2, t, f3);
xlabel('t')
legend(['f1'; 'f2'; 'f3'])
grid on
% The remainder of this file demonstrates the FFT technique to convert
% signals between time and frequency domain.
% Material parameters (Debye model)
alpha = 1;
tau = 1;
% Material parameters (Lorentz model)
omegap = 1;
omega0 = 2;
nu = 1;
nu0 = sqrt(omega0^2 - nu^2/4);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Time domain
Nt = 1000; % Number of time steps
T = 10; % Stop time
t = linspace(0, T, Nt); % A vector of Nt values uniformly spaced from 0 to T
dt = t(2) - t(1); % The step between two time values
chi_debye_t = alpha*exp(-t/tau);
chi_lorentz_t = omegap^2/nu0*exp(-nu*t/2).*sin(nu0*t);
% Plot the time domain susceptibilities
figure()
clf
subplot(2,1,1)
plot(t, chi_debye_t)
xlabel('t')
ylabel('\chi(t)')
grid on
title('Debye model')
subplot(2,1,2)
plot(t, chi_lorentz_t)
xlabel('t')
ylabel('\chi(t)')
grid on
title('Lorentz model')
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Frequency domain
% First compute the result by FFT. To change the frequency resolution,
% you can change the factor 10 in Nfft = Nt*10.
Nfft = Nt*10;
chi_debye_f_doublesided = fft(chi_debye_t, Nfft)*dt;
chi_lorentz_f_doublesided = fft(chi_lorentz_t, Nfft)*dt;
% The FFT computes also results for negative frequencies, but here we only
% need the positive ones, stored in the first half of the vector.
chi_debye_f = chi_debye_f_doublesided(1:Nfft/2+1);
chi_lorentz_f = chi_lorentz_f_doublesided(1:Nfft/2+1);
% The line below is how to compute the frequency points for an FFT
f = 1/(2*dt)*linspace(0, 1, Nfft/2+1);
omega = 2*pi*f;
% Compute the analytical result to compare with
chi_debye_f_analytical = alpha*tau./(1 + 1j*omega*tau);
chi_lorentz_f_analytical = omegap^2./(-omega.^2 + omega0^2 + 1j*omega*nu);
% Plot the frequency domain susceptibilities
omega_max = max([1/tau, omegap, omega0, nu])*3; % Estimate the largest frequency of interest
figure()
clf
subplot(2,1,1)
hold on
plot(omega, real(chi_debye_f), 'b-')
plot(omega, imag(chi_debye_f), 'r-')
plot(omega, real(chi_debye_f_analytical), 'b--')
plot(omega, imag(chi_debye_f_analytical), 'r--')
xlim([0, omega_max])
xlabel('\omega')
ylabel('\chi(\omega)')
title('Debye model')
grid on
hold off
subplot(2,1,2)
hold on
plot(omega, real(chi_lorentz_f), 'b-')
plot(omega, imag(chi_lorentz_f), 'r-')
plot(omega, real(chi_lorentz_f_analytical), 'b--')
plot(omega, imag(chi_lorentz_f_analytical), 'r--')
xlim([0, omega_max])
xlabel('\omega')
ylabel('\chi(\omega)')
title('Lorentz model')
grid on
hold off
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%