%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% LINEAR SYSTEMS TUTORIAL 
%%% EPS, 3/99.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% The "Impulse" vector, delta[n], is 1 only when n=0. Here we plot
% delta[n-16]:
impulse = mkImpulse([1,32],[1,16]);
lplot(impulse);

% Linear shift-invariant systems are computed via the operation of
% "convolution".  Matlab provides a convolution function called
% "conv":

sz = 16;
% Step vector, u[n], is 0 when n<0.  Here we plot u[n-8]:
signal = zeros(sz,1);
signal((sz/2+1):sz) = ones(sz/2,1);
subplot(2,1,1); lplot(signal);

blurfilter = [1 4 6 4 1]'/16;
bsignal = conv(signal, blurfilter);
subplot(2,1,2); lplot(bsignal);

% Matlab's default conv function produces a larger vector than than
% the input signal.  We have provided a function upConv that always
% returns a signal of the same size as the input.  The boundaries can
% be handled in many different ways (try "help upConv" to see the
% options).  This allows a more direct comparison of the input and
% output signals.

bsignal = upConv(signal, blurfilter, 'zero')
subplot(2,1,2); lplot(bsignal);

% Now we want to look at the matrix M corresponding to this linear
% operation.  Consider the response of the system to an impulse
% located in the first sample.  By the rules of matrix multiplication,
% this must be equal to the first column of the matrix M.  Similarly,
% the nth column of the matrix M must contain the response to an
% impulse in location n.  Since the system is shift-invariant, the
% responses to these impulses are just shifted copies of the same
% vector.  Thus, we can characterize the full system response by just
% knowing the response to an impulse, and we can construct the matrix
% by concatenating these columns:

M = zeros(sz,sz);
for n=1:sz
  M(:,n) = upConv(mkImpulse([sz,1], [n,1]), blurfilter, 'zero');
end
M

% Compare the matrix multiplication to the previous convolution
% result, using the "imStats" function:
msignal = M*signal;
subplot(2,1,1); lplot(bsignal);
subplot(2,1,2); lplot(msignal);
imStats(msignal, bsignal)

% We can view the matrix M as an image using our command showIm.  The
% matrix is square, and has its non-zero entries in a band along the
% diagonal.  This is known as a "Toeplitz" matrix.
clf; showIm(M)

% The previously considered convolution filter was a blurring or
% "lowpass" filter.  Now consider a filter that computes the
% derivative of a blurred signal.  The blurring function in this case
% is a Gaussian.  Since differentiation is also linear and
% shift-invariant, and since convolution is assocative, we can combine
% the two operations (blurring and differentiation) into a single
% convolution kernel:

Xvals = [-3:3]'; sig = 1.2;
derivfilter = Xvals .* exp( -Xvals.^2 ./ (2 * sig^2));
clf; lplot(derivfilter);

subplot(2,1,1); lplot(signal);
subplot(2,1,2); lplot(upConv(signal,derivfilter,'zero'));


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% Sinusoids

% Sinusoidal and cosinusoidal vectors play a particularly important
% role in representing signals, because complex exponential vectors
% (including sines and cosines) are the eigenfunctions of linear
% shift-invariant systems.  In other words, applying a linear
% shift-invariant system to a sinusoidal vector produces another
% sinusoidal vector of the same frequency.  Only the phase and
% amplitude of the output sinusoid may change.

% A sinusoidal vector with period 8:
nRange = [0:31]';
amplitude = 1; phase_off = 0; period = 8; freq = 2*pi/period;
sinusoid = amplitude .* sin(freq .* nRange + phase_off);
subplot(2,1,1); lplot(sinusoid);

% Notice that for discrete sinusoids, unlike continous sinusoids,
% adding 2*pi to the frequency gives the same sinusoid:
freq = 2*pi/period + 2*pi;
sinusoid_shift_2pi = amplitude .* sin(freq .* nRange + phase_off);
subplot(2,1,2); lplot(sinusoid_shift_2pi);

% The importance of this is that we need only consider frequencies in
% a frequency interval of length 2*pi such as -pi to pi.  

% Also notice that although continuous sinusoids with frequency w are
% periodic with period 2*pi/w, this is not necessarily true of
% discrete sinusoids.  For example, a discrete sinusoid with frequency
% w=1 is NOT periodic:
period = 2*pi; freq = 2*pi/period;
non_periodic_sinusoid = amplitude .* sin(freq .* nRange + phase_off);
figure(1), clf
plot(nRange,non_periodic_sinusoid,'r-',nRange,non_periodic_sinusoid,'go');
axis([0 31 -1.2 1.2]);

% Why isn't this sequence periodic?  Is it because we've plotted
% only 32 samples?  If we were to plot more samples, would it
% ever repeat?

% For a finite length sequence, we have an even more stringent
% requirement.  By a periodic finite length sequence, we mean
% circularly periodic.  When you go off the end you start back at
% the beginning.  To be periodic, the sequence length must be a
% multiple of the period.

% Finally, consider the response of a random filter to a sinusoid (do this
% a few times), and notice what has changed:
filt = rand(5,1); 
result = upConv(sinusoid,filt,'circular'); 
subplot(2,1,2); lplot(result);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% Discrete Fourier Transforms

% Any signal can be expressed as a (weighted) linear sum of impulses.
% Likewise, a signal can be expressed as a (weighted) linear sum of
% sines and cosines.  More specifically, if we treat our signal vector
% (of length N) as one period of an infinite periodic sequence, then
% it may be written as a sum of a fixed set of N sinusoids that are
% periodic with period N.  These frequencies are: 2*pi*k/N for
% k=0,1,...,N-1.  In our examples, N=32, so the frequencies are: 0,
% pi/16, 2pi/16, 3pi/16,..., 31pi/16.

% Alternatively, remember that adding 2*pi to the frequency doesn't
% change the sinusoid, and thus we could just as well use a set of N
% sinusoids indexed as: 2*pi*k/N for k=-N/2,...,-1,0,1,...,N/2.  In
% our examples, these periods are:
% -pi,-15pi/16,...,-pi/16,0,pi/16,...,15pi/16.  Take a look at some of
% these sinusoids and cosinusoids to see that these frequencies are
% all distinct.

% As an example, we construct a Gaussian as a sum of cosines (where we
% assume the origin of the spatial axis is at sample 17):
gaussian = exp(-((nRange-16).^2)./(2*4^2));
subplot(2,1,1); lplot(gaussian);

% Gradually add in the sinusoidal components to construct this Gaussian:
gaussian_series = (10.0258 * cos(2*pi*0/32*(nRange-16))/32);
subplot(2,1,2); lplot(gaussian_series); imStats(gaussian, gaussian_series);
gaussian_series = gaussian_series + (7.3662 * cos(2*pi*1/32*(nRange-16))/16);
subplot(2,1,2); lplot(gaussian_series); imStats(gaussian, gaussian_series);
gaussian_series = gaussian_series + (2.9192 * cos(2*pi*2/32*(nRange-16))/16);
subplot(2,1,2); lplot(gaussian_series); imStats(gaussian, gaussian_series);
gaussian_series = gaussian_series + (0.6252 * cos(2*pi*3/32*(nRange-16))/16);
subplot(2,1,2); lplot(gaussian_series); imStats(gaussian, gaussian_series);
gaussian_series = gaussian_series + (0.0716 * cos(2*pi*4/32*(nRange-16))/16);
subplot(2,1,2); lplot(gaussian_series); imStats(gaussian, gaussian_series);

% The vector of weights [10.0258 7.3662 2.9192 0.6252 0.0716]
% corresponds to the first four coefficients of the Discrete Fourier
% Transform (DFT).  They can be efficiently computed using an
% algorithm known as the Fast Fourier Transform (FFT).

% MatLab provides a function "fft" to compute these coefficients.
% We have to first periodically shift the signal around so that the
% origin is located at the first sample.  The matlab function
% "fftshift" does this:
shift_gaussian = fftshift(gaussian);
lplot(shift_gaussian);

% The central sample now corresponds to the zero-frequency  (also
% known as the "DC") coefficient.

% Complex numbers provide a convenient bookkeeping notation for the
% DFT coefficients, and are the standard method of representing them.
% The magnitude of each complex number gives the amplitude of the
% associated sinusoid.  The angle of each complex number gives the
% phase.  Alternatively, the real part of each number corresponds to
% the cosine coefficient and and the imaginary part to the sine. 

% Since our gaussian signal is symmetric about the origin, it is
% represented with only cosine components (phase = 0), and thus the
% coefficients all have zero imaginary part:
dft = fft(shift_gaussian)
lplot(abs(dft))  % "abs" computes magnitude of complex values

% If we shift the input vector, this changes the phases of the
% associated Fourier components, but not their magnitudes.  Thus:
dft = fft(shift(gaussian,3))
lplot(abs(dft))

% Although Matlab computes the DFT with the origin (zero frequency) at
% the first sample, it is common to plot the DFT with the origin in
% the middle.  We can accomplish this by fftshift'ing the result:
lplot(abs(fftshift(dft)))

% Now look at some more example DFTs. First, an impulse:
sz = 32;
signal = mkImpulse([sz,1]);
dft = fft(signal);
subplot(2,1,1); lplot(signal);
subplot(2,1,2); lplot(abs(fftshift(dft)))

% Fourier transform of a constant function:
signal = ones([sz,1]);
dft = fft(signal);
subplot(2,1,1); lplot(signal);
subplot(2,1,2); lplot(abs(fftshift(dft)))

% A cosine (remember, the origin is being plotted at sample 17):
num_cycles = 4;
signal = cos(2*pi*num_cycles/sz*nRange);
dft = fft(signal);
subplot(2,1,1); lplot(signal);
subplot(2,1,2); lplot(abs(fftshift(dft)))

% We get a pair of impulses in the transform.  One of the impulses
% corresponds to the frequency of the signal (4 cycles per image) at
% position 4 in the transform.  The other shows up at position -4.
% Why is there a second impulse?   [Hint: how do you write a cos(phi) in
% terms of e^{i phi} and e^{-i phi}?]

% A sine of the same frequency:
signal = sin(2*pi*num_cycles/sz*nRange);
dft = fft(signal);
subplot(2,1,1); lplot(signal);
subplot(2,1,2); lplot(abs(fftshift(dft)))

% A higher-frequency cosine:
num_cycles = 6;
signal = cos(2*pi*num_cycles/sz*nRange);
dft = fft(signal);
subplot(2,1,1); lplot(signal);
subplot(2,1,2); lplot(abs(fftshift(dft)))

% Convolution in the signal domain is the same as multiplication in the
% frequency domain, and vice versa.  This "convolution theorem" is
% extremely useful.  Some filters are simple to characterize in the
% frequency domain, but complicated in the signal domain.  For example,
% the filter may be very compact in the frequency domain (i.e., zero
% nearly everywhere), but very big (i.e., lots of samples needed) in
% the signal domain.  In such cases, you can do the filtering by Fourier
% transforming the signal, multiplying in the frequency domain, and
% then Fourier transforming back.

% Besides the issue of characterization, another reason to compute
% convolutions this way is efficiency.  The FFT requires a number of
% operations proportional to N log( N), where N is the size of the
% input signal.  Computing a convolution directly might cost N^2
% operations.  As N grows, the Fourier transform has a larger
% computational advantage.  If the convolution kernel is small,
% however, it is usually more efficient to compute the convolution
% directly.

% The convolution theorem is also helpful in understanding the Fourier
% transforms of different functions.  Consider a "Gabor" function
% (product of a Gaussian and a sinusoid).  The DFT is the convolution
% of the Gaussian and sinusoid DFTs:
nRange = [0:31]';
gaussian = exp(-((nRange-16).^2)./(2*4^2));
sinusoid = mkSine(size(gaussian),32/6,pi/2);
signal = gaussian .* sinusoid;
dft = fftshift(fft(signal));
conv_dft = cconv2(fft(gaussian),fft(sinusoid),1)/32;
subplot(3,1,1); lplot(signal);
subplot(3,1,2); lplot(abs(dft),[-16 15])
subplot(3,1,3); lplot(abs(conv_dft),[-16 15])

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%% Symmetry properties of Fourier Transform:

% For any real-valued, antisymmetric (odd) function, in which
% f(x) = -f(-x), the real part of the FT is zero, and the
% imaginary part of the FT is antisymmetric (odd).  For any
% real-valued, symmetric (even) function, in which f(x) = f(-x),
% the imaginary part of the FT is zero and the real part is
% symmetric (even).

sz = 64;
random_signal = 0.5-rand(1,sz);

% Decompose random_signal into even and odd components:
%   even_signal(x) = 0.5*(random_signal(x)+random_signal(-x))
%   odd_signal(x)  = 0.5*(random_signal(x)-random_signal(-x))
%
% One can easily verify that
%   random_signal(x) = even_signal(x)+odd_signal(x).
%
% Note that the origin (i.e. x=0) is at location 33 
% (i.e. random_signal(33)).

figure(1);
even_signal = random_signal + random_signal(1+mod([sz:-1:1],sz));
subplot(2,1,1); plot(-32:31,even_signal,'r-',-32:31,even_signal,'go');
axis([-32 31 -0.5 0.5]);

odd_signal =  random_signal - random_signal(1+mod([sz:-1:1],sz));
subplot(2,1,2); plot(-32:31,odd_signal,'r-',-32:31,odd_signal,'go');
axis([-32 31 -0.5 0.5]);

figure(2);
% Real part of Fourier transform of an even signal is
% even-symmetric; imaginary part is zero.
ft_even_signal = fftshift(fft(fftshift(even_signal)));
real_ft_even = real(ft_even_signal);
subplot(2,1,1); plot(-32:31,real_ft_even,'r-',-32:31,real_ft_even,'go');
axis([-32 31 -8 8]);

imag_ft_even = imag(ft_even_signal);
subplot(2,1,2); plot(-32:31,imag_ft_even,'r-',-32:31,imag_ft_even,'go');
axis([-32 31 -8 8]);

% Imaginary part of Fourier transform of an odd signal
% is odd-symmetric; real part is zero.
ft_odd_signal = fftshift(fft(fftshift(odd_signal)));
real_ft_odd = real(ft_odd_signal);
subplot(2,1,1); plot(-32:31,real_ft_odd,'r-',-32:31,real_ft_odd,'go');
axis([-32 31 -8 8]);

imag_ft_odd = imag(ft_odd_signal);
subplot(2,1,2); plot(-32:31,imag_ft_odd,'r-',-32:31,imag_ft_odd,'go');
axis([-32 31 -8 8]);

% For any real-valued signal, the real part of the FT is
% even-symmetric and the imaginary part of the FT is
% odd-symmetric.
ft_random_signal = fftshift(fft(random_signal));
real_ft = real(ft_random_signal);
subplot(2,1,1); plot(-32:31,real_ft,'r-',-32:31,real_ft,'go');
axis([-32 31 -5 5]);

imag_ft = imag(ft_random_signal);
subplot(2,1,2); plot(-32:31,imag_ft,'r-',-32:31,imag_ft,'go');
axis([-32 31 -5 5]);

% Taken together, these symmetry properties mean that there is
% quite a lot of redundancy in the FT of a real signal.  A simple
% way to count the amount of redundancy is to compare the number
% of samples.  Take a real-valued input signal with 64 samples.
% Computing its fft gives a total of 128 samples (half in the
% real part and half in the imaginary part), a factor of 2
% redundant.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%% Circular Shifting: 

% If we shift the signal as if it were periodic (i.e., translate
% the signal, wrapping around at the edges), this does not affect
% the Fourier transform magnitude:

fft_gauss = fftshift(fft(fftshift(gaussian)));
fft_shift_gauss = fftshift(fft(fftshift(shift(gaussian,3))));

% Magnitudes are the same:
mag_fft_gauss = abs(fft_gauss);
subplot(2,1,1); plot(-16:15,mag_fft_gauss,'r-',-16:15,mag_fft_gauss,'go');
axis([-16 15 -.2 12]);

mag_fft_shift_gauss = abs(fft_shift_gauss);
subplot(2,1,2); plot(-16:15,mag_fft_shift_gauss,'r-',-16:15,mag_fft_shift_gauss,'go');
axis([-16 15 -.2 12]);

% Should be zero:
imStats(mag_fft_shift_gauss,mag_fft_gauss)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%% Differentiation:

% Taking a derivative of a signal in time is the same as multiplying
% by an imaginary ramp in frequency. In particular, 
%                  Fourier{d/dx f(x)} = -i w Fourier{f(x)}
% where i = sqrt(-1) and w is normalized frequency.
% For an intuition for the derivative property, recall that the
% d/dx[cos(wx)] = -w sin(wx).

% For example, let's consider a Gaussian and the first derivative
% of a Gaussian.

figure(1);
gaussian = exp(-(((nRange-16).^2)./(2*4^2)));
subplot(2,1,1); plot(-16:15,gaussian,'r-',-16:15,gaussian,'go');
axis([-16 15 -.2 1.2]);

gaussian_deriv = -2/(2*4^2).*(nRange-16).*gaussian;
subplot(2,1,2); plot(-16:15,gaussian_deriv,'r-',-16:15,gaussian_deriv,'go');
axis([-16 15 -.2 .2]);

figure(2);
ft_gaussian = fftshift(fft(fftshift(gaussian)));
real_ft_gaussian = real(ft_gaussian);
subplot(2,1,1); plot(-16:15,real_ft_gaussian,'r-',-16:15,real_ft_gaussian,'go');
axis([-16 15 -.2 12]);

ft_gaussian_deriv = fftshift(fft(fftshift(gaussian_deriv)));
imag_ft_gaussian_deriv = imag(ft_gaussian_deriv);
subplot(2,1,2); plot(-16:15,imag_ft_gaussian_deriv,'r-',-16:15,imag_ft_gaussian_deriv,'go');
axis([-16 15 -2 2]);

% ramp := i w :
ramp = 2*pi/32*sqrt(-1).*(nRange-16);

% Compute the Fourier transform of the derivative of a gaussian
% by multiplying the Fourier transform of a gaussian by the ramp.
ft_gaussian_mul_ramp = ramp .* ft_gaussian;
imag_ft_gaussian_mul_ramp = imag(ft_gaussian_mul_ramp);
subplot(2,1,1);
plot(-16:15,imag_ft_gaussian_mul_ramp,'r-',-16:15,imag_ft_gaussian_mul_ramp,'go');
axis([-16 15 -2 2]);

% Should be zero:
max(abs(ft_gaussian_deriv - ft_gaussian_mul_ramp))