%% The Fourier Transform
% Table of Contents
% Moving into the Frequency Domain
%     Introducing Frequency Space 1D Audio Example:
% The Frequency Domain
% Applications of the Fourier Transform
% Fourier representation of a Sin Wave
% Fourier Transfrom: A 2D Image Example
% Visualising Frequency Domain Transforms
%     Sinusoidal Decomposition
%         So What Does All This Mean?
% Fundamental Application: Filtering
%     The Basic Idea of Filtering in Frequency Space
% Visualising the Frequency Domain
%     So are we ready for the Fourier Transform?
% Fourier Transform: Introducing The Fourier Transform
%     1D Fourier Transform
%         Inverse 1D Fourier Transform
% Fourier Transform Example
% The Sinc Function: The Fourier Transform of the Top Hat Function
% The 2D Fourier Transform
% The Discrete Fourier Transform
%     1D Discrete Fourier transform
%     2D Discrete Fourier transform
%         Balancing the 2D DFT
% The Fourier Transform in MATLAB
%         FFT for Spectral Analysis
%         See also: MATLAB Central's FFT Easy Demo
% Visualising the Fourier Transform
%     Magnitude Spectrum
%     Phase Spectrum
%     Relating a Sample Point to a Frequency Point
% The Spectrogram
%     Another MATLAB spectrogram Example: Aphex Twin
% Properties of Fourier Transforms
%     Linear Operator Theorem
%         Fourier Transform Linear Operator MATLAB Demo: FFT_Linear.m
%     Shifting Operator Theorem
%         Fourier Transform Shifting Operator MATLAB Demo: FFT_Shift.m
%     Scaling Theroem
%         Fourier Transform Scaling Operator MATLAB Demo: FFT_Scale.m
%     2D Fourier Transform: Rotation
%         Fourier Transform Rotate Operator Demo: FFT_Rotate.m
%     Fourier Transform: Zeroth Component
%         MATLAB Illustration of zeroth element
% Centring the Frequency of a Fourier Transform
% Supporting functions
%
%
%% Moving into the Frequency Domain
%%
% So far we have consider the representation of a signal as amplitudes of sound (1D, Audio) or image intensities (2D image) or 3D video frames (2D images sampled over time).
% There is another way to represent data: the data's underlying frequency componenents or content.
%% Introducing Frequency Space 1D Audio Example:
% Let's consider a 1D (e.g. Audio) example to see what the different domains mean:
% Consider a complicated sound such as the a chord played on a piano or a guitar.
%         We can describe this sound in two related ways:
%                     Temporal Domain: Sample the amplitude of the sound many times a second, which gives an approximation to the sound as a functionof time.
%
%
%                       Frequency Domain: Analyse the sound in terms of the pitches of the notes, or frequencies, which make the sound up, recording the amplitude of each frequency.
%
% The human ear perceived 4 Fundamental Frequencies for this chord:
% D flat: 554.40Hz
% F: 698.48Hz
% A flat: 830.64Hz
% C:1046.56Hz
% People with a musical ear may be able to discern the notes by name or by comparison to pitches played on an instrument.
% No normal person can determine the precise frequency or know the precise frequency of not names!
%                         This is not how humans/musicians work
%                         No one remebers exact frequencies to this accuracy but humans like labels, the note names
% Note: harmonics/partial frequencies are not inculded in the above discussion:
%                         The wave form of the chord played on a piano will sound and look different to one played on a guitar or any other instument
%                         Humans percived the lowest Fundamental Frequencies well from all musical instruments --- for polyphoninc audio (chords and/or concurrenrt playing instruments) this is still a very difficult machine learning task.
%                         Harmonics/partial frequencies at certain discrete intervals higher that the Fundamental Frequencies give the discerning information to hear different instruments. All instruments are different, and it may depend on how they a played too (e.g. bow or plucked violin).
%                        Not all chords sound the same --- at basic level playing just Fundamental Frequencies would  all be playing sine waves which do not like or anywhere near as nice a piano (or other instrtument)!
%
%% The Frequency Domain
%
% The Frequency domain can be obtained through the transformation, via Fourier Transform (FT) (There are other transforms: for example, the DCT will be met later), from one  domain
% Temporal (Time) or  Spatial domain
% to the anoother
% Frequency Domain
% In the Frequency Domain, we do not think in terms of signal or pixel intensities but rather underlying sinusoidal waveforms of varying frequency, amplitude and phase representing the underlying data (now suitably transformed).
%
%% Applications of the Fourier Transform
%
% There many applications of Fourier Transform, these Applications include (to name but a few):
% Essential tool for Engineers, Physicists, Mathematicians and Computer Scientists in many areas.
% Fundamental tool for Digital Signal Processing and Image Processing
% Many types of Frequency Analysis: JPEG use related techniques.
% Many more \ldots\ldots
%
% FT of some simple shapes
%%
%
%% Fourier representation of a Sin Wave
% Back to Basics for a moment
% An 8 Hz Sine Wave A signal that consists of a sinusoidal wave at 8 Hz.
% 8 Hz means that wave is completing 8 cycles in 1 second
% The  frequency of that wave is 8 Hz.

amp =1; phase = 0;
freq = 8;
samp_freq = 256;
dur = 256;
y = mysin(amp,freq,phase,dur, samp_freq);
figure;
plot(y)

% From the frequency domain we can see that the composition of our signal is
% one peak occurring with a frequency of 8 Hz — there is only one sine wave here.
% with a magnitude/fraction of 1.0 i.e. it is the whole signal.

y_fft = fft(y)*2/dur; % Normalised "2 two-sided" FFT
% Plot "zoom in" on  first half of FFT Spectra for this illustration's purpose
% See below for full explanation
figure
f=((0:dur-1)/dur)*samp_freq;
axis([-0.2*samp_freq/16 max(f(1:50)) -0.1 1.1]);


stem(f(1:50),abs(y_fft(1:50))); 
legend('Magnitude spectrum |y|'); 
xlabel('f in Hz \rightarrow');

% Note: the MATLAB fft() function computes the Fourier transform – we'll formally introduce shortly (see below).
%
%% Fourier Transfrom: A 2D Image Example
% What do frequencies in an image mean?
% Images are, really, no more complex than 1D data signals:
% Brightness along a line can be recorded as a set of values measured at equally spaced distances apart,
%  Or equivalently, at a set of spatial frequency values.
% Each of these frequency values is a frequency component.
% An image is a 2D array of pixel measurements.
% We form a 2D grid of spatial frequencies.
% A given frequency component now specifies what contribution is made by data which is changing with specified x and y direction spatial frequencies.
% Large values at high frequency components then the data is changing rapidly on a short distance scale.
%                     e.g. a page of text
%                       However, Noise contributes (very) High Frequencies also
% Large low frequency components then the large scale features of the picture are more important. e.g. a single fairly simple object which occupies most of the image.
%
%% Visualising Frequency Domain Transforms
%
%% Sinusoidal Decomposition
% Any  digital signal (function) can be decomposed into purely sinusoidal components
% Sine waves of different size/shape — varying amplitude, frequency and phase.
% When added back together they reconstitute the original signal.
% The Fourier transform is the tool that performs such an operation.
%
% Summing Sine Waves. Example: to give a Square(ish) Wave (E.g. Additive Synthesis)
% Digital signals are composite signals made up of many sinusoidal frequencies
% A 200Hz digital signal (square(ish) wave) may be a composed of 200, 600, 1000, etc. sinusoidal signals which sum to give:
%
%% So What Does All This Mean?
% Transforming a signal into the frequency domain allows us:
% To see what sine waves make up our underlying signal
%  E.g.
%                 One part sinusoidal wave at 50 Hz and
%                  Second part sinusoidal wave at 200 Hz.
%                  Etc.
% More complex signals will give more complex decompositions but the idea is exactly the same.
%
%% Fundamental Application: Filtering
% How is this all this sinisoidal decompostion useful?
%
% A fundemental task in Signal and Image Process is filtering.
%     We have met filtering in real space with IIR/FIR filters
%     We now develop filtering in frequency space
%                 I think filtering in frequecncy space is much more intuitive
%% The Basic Idea of Filtering in Frequency Space
%
% In frequency space, filtering  involves attenuating or removing certain frequencies — easily performed:
% Low pass filter —
%         Ignore high frequency noise components — make zero or a very low value.
%     Only store lower frequency components
% High Pass Filter — opposite of above
% Bandpass Filter — only allow frequencies in a certain range.
%%
%
%% Visualising the Frequency Domain
% Think Graphic Equaliser
% An easy way to visualise what is happening is to think of a graphic equaliser on a stereo system (or some software audio players, e.g. iTunes).
%
%
% As we will see later (See Digital Audio Effects Section), an graphic equaliser is an implementaion of a bank of bandpass equalisers, which are implemented via filters
%
%% So are we ready for the Fourier Transform?
% We have all the Tools....
% This lecture, so far, (hopefully) set the context for Frequency decomposition.
% Past CM2104 maths lectures:
%                     Odd/Even Functions: \sin(-x) = -\sin(x), \cos(-x) = \cos(x)
%                     Complex Numbers: Phasor Form r e^{i\phi} = r(\cos\phi + i \sin\phi)
%                     Calculus Integration:  \int e^{kx}\:dx = \frac{e^{kx}}{k}
% Digital Signal Processing:
%             Basic Waveform Theory. Sine Wave
% y = A. sin( 2\pi.n.F_w/F_s)
% where:  A = amplitude,  F_w = wave frequency,  F_s = sample frequency, n is the sample index.
%             Relationship between Amplitude, Frequency and Phase:
%
%             Cosine is a Sine wave 90^\circ out of phase
%             Impulse Responses
% DSP + Image Proc.: Filters and other processing, Convolution (Recapped later in this lecture)
%
%% Fourier Transform: Introducing The Fourier Transform
%
%
% 1. The tool which converts a spatial or temporal (real space) description of audio/image data, for example, into one in terms of its frequency components is called the Fourier transform
%
% 2. The new version is usually referred to as the Fourier space description of the data.
% We then essentially process the data:
%           E.g. for filtering basically this means attenuating or setting certain frequencies to zero
%
% 3. We then need to convert data back (or invert) to real audio/imagery to use in our applications.
% The corresponding inverse transformation which turns a Fourier space description back into a real space one is called the inverse Fourier transform.
%% 1D Fourier Transform
%
% Our typical example in the 1D Case is an Audio Signal, but any one 1D signal is treated similarly (e.g. time series data, such as financial stock trends, Heart ECG, Brain EEG)
% Consider a continuous function,f(x), of a single variable x representing distance (or time).
% The Fourier transform of that function is denoted F(u), where u represents spatial (or temporal) frequency is defined by:
% F(u) = \int_{-\infty}^\infty f(x) {\bf e^{-2\pi i x u}}\,dx.
% Note: In general, F(u) will be a complex quantity even though the original data is purely real.
% The meaning of this is that not only is the magnitude of each frequency present important, but that its phase relationship is too. Recall Phasors from Complex Number Lectures (CM2208).
% e^{-2\pi i x u} above is a Phasor.
%% Inverse 1D Fourier Transform
%
% Inverse 1D Fourier Transform The inverse Fourier transform for regenerating f(x) from F(u) is given by:
% f(x) = \int_{-\infty}^\infty F(u) e^{{\bf  2\pi i x u}}\,du
% which is rather similar to the (forward) Fourier transform
% except that the exponential term has the opposite sign – It is not negative
%% Fourier Transform Example
%
% Fourier Transform of a Top Hat Function Let's see how we compute a Fourier Transform: consider a particular function f(x) defined as:
% f(x)   =   1 if \vert x\vert \leq 1
%                   0 otherwise,
%
%
%% The Sinc Function: The Fourier Transform of the Top Hat Function
%
% We derive the Sinc function by taking the Fourier Transform of the Top Hat Function
% So its Fourier transform is defined by:
% F(u)  =  \int_{-\infty}^\infty f(x) e^{-2\pi i x u}\,dx
%      =  \int_{-1}^1 1\times e^{-2\pi i x u}\,dx
%            =  \frac{-1}{2\pi i u}(e^{2\pi i u} - e^{-2\pi i u})
% Now, recall that
%  \sin \theta = \frac{e^{i\theta} - e^{-i\theta}}{2i}
% So,
% F(u) = \frac{\sin 2 \pi u}{\pi u}.}
%
% Note: In this (special) case, F(u) is purely real, which is a consequence of the original data being symmetric in x and -x.
% f(x) is an even function.
% The Fourier transform of a top hat function, the Sinc function, looks like this:

warning off;
fplot(@(u) sin(2*pi*u)/(pi*u), [-50 50]);

%
% MATLAB built in sinc() function:

showDoc = 0;
if showDoc
    doc sinc;
end

% Note: the sinc function has some important theoritical properties for filters (a top hat is essentially a band bass "brick wall" filter function).
% Further investiagtion is not that relevant to, and beyond the scope of, this module though.
%
% Create a sinc function

Fs=42;Ts=1/Fs;
t=-1:Ts:40*Ts-Ts;
f=5;
y=sinc(pi*t*f);
figure;
plot(t,y);
xlabel('x');
ylabel(' magnitude');

% Compute FFT of sinc, hey presto a top hat:

 N=512;
fy=(fft(y,N));
figure;
fr=(0:N-1)*Fs/N;
 %fy=fy(1:floor(end/2));
 %fr=fr(1:floor(end/2));
 plot(fr,fftshift(abs(fy)));
 xlabel(' frequency');
 ylabel(' magnitude');

%
%% The 2D Fourier Transform
%
% For the purposes of the modules, the 2D Case is most typically represented by Image data.
%
% If f(x,y) is a function, for example intensities in an image, its Fourier transform is given by
% F(u, v) = \int_{{-\infty}}^{\infty} \int_{{-\infty}}^{\infty}f (x, y) e^{-2\pi i (x u+y v)} dx dy
% and the inverse transform, as might be expected, is
%  f (x, y) = \int_{{-\infty}}^{\infty}\int_{{-\infty}}^{\infty} F(u, v)e^{2\pi i (x u+y v)} du dv
%
%% The Discrete Fourier Transform
%
% All our Audio and Image data are digitised so we need to move to the Discrete Fourier Transform.
% Now, we need a discrete formulation of the Fourier transform. This:
% Assumes regularly spaced data values, and
% Returns the value of the Fourier transform for a set of values in frequency space which are equally spaced.
% This is done quite naturally by replacing the integral by a summation, to give the  discrete Fourier transform or DFT for short.
%% 1D Discrete Fourier transform
%
% Consider the 1D Case: In 1D it is convenient now to assume that x goes up in steps of 1, and that there are N samples, at values of x from 0 to N-1.
%
% So the DFT takes the form
% F(u) = \frac{1}{N} \sum_{{x=0}}^{{N-1}} f (x)e^{-2\pi ixu/N}
%
% while the inverse DFT is
% f (x) = \sum_{{x=0}}^{{N-1}}F(u)e^{2\pi ixu/N}.
%
% NOTE: Minor changes from the continuous case are a factor of 1/N in the exponential terms, and also the factor 1/N in front of the forward transform which does not appear in the inverse transform.
%
%% 2D Discrete Fourier transform
%
% The 2D DFT works in a similar fashion:
% So for an N\times M grid in x and y we have
% F(u,v) = \frac{1}{N M}\sum_{{x=0}}^{{N-1}} \sum_{{y=0}}^{{M-1}}f (x, y)e^{-2\pi i(xu/N + yv/M)
% and
% f (x, y) =\sum_{{u=0}}^{{N-1}} \sum_{{v=0}}^{{M-1}}F(u, v)e^{2\pi i(xu/N+yv/M)}
% .
%
%% Balancing the 2D DFT
% Most Images are Square Often N=M, and it is then it is more convenient to redefine F(u,v) by multiplying it by a factor of N, so that the forward and inverse transforms are more symmetric:
%  F(u, v) = \frac{1}{N} \sum_{{x=0}}^{{N-1}}\sum_{{y=0}}^{{N-1}} f (x, y)e^{-2\pi i(xu+yv)/N}
%  and
%  f (x, y) = \frac{1}{N} \sum_{{u=0}}^{{N-1}} \sum_{{v=0}}^{{M-1}}F(u, v)e^{2\pi i(xu/N+yv/ M)}.
%
%% The Fourier Transform in MATLAB
% fft() and fft2() MATLAB provides functions for 1D and 2D Discrete Fourier Transforms (DFT):
% fft(X) is the 1D discrete Fourier transform (DFT) of vector X. For matrices, the FFT operation is applied to each column— NOT a 2D DFT transform.

showDoc = 1;
if showDoc
    doc fft;
end

%
% fft2(X) returns the 2D Fourier transform of matrix X. If X is a vector, the result will have the same orientation.

showDoc = 0;
if showDoc
    doc fft2;
end

%
% fftn(X) returns the N-D discrete Fourier transform of the N-D array X.

showDoc = 0;
if showDoc
 doc fftn;
end

%
% Inverse DFT: ifft(), ifft2(), ifftn() perform the inverse DFT. See appropriate MATLAB help/doc pages above for full details.

showDoc = 0;
if showDoc
    doc ifft;
end

%

showDoc = 0;
if showDoc
    doc ifft2;
end

%

showDoc = 0;
if showDoc
    doc ifftn;
end

%
% Plenty of examples to Follow.
%
%% FFT for Spectral Analysis
%
% This example shows (Borrowed from MATLAB's built in fftdemo.m ('doc FFT for Spectral Analysis') the use of the FFT function for spectral analysis. A common use of FFT's is to find the frequency components of a signal buried in a noisy time domain signal.
%
% First create some data. Consider data sampled at 1000 Hz. Start by forming a time axis for our data, running from t=0 until t=.25 in steps of 1 millisecond. Then form a signal, x, containing sine waves at 50 Hz and 120 Hz.

t = 0:.001:.25;
x = sin(2*pi*50*t) + sin(2*pi*120*t);

%
% Add some random noise with a standard deviation of 2 to produce a noisy signal y. Take a look at this noisy signal y by plotting it.

y = x + 2*randn(size(t));
plot(y(1:50))
title('Noisy time domain signal')

%
% Clearly, it is difficult to identify the frequency components from looking at this signal; that's why spectral analysis is so popular.
% Finding the discrete Fourier transform of the noisy signal y is easy; just take the fast-Fourier transform (FFT).

Y = fft(y,251);

%
% Compute the power spectral density, a measurement of the energy at various frequencies, using the complex conjugate (CONJ). Form a frequency axis for the first 127 points and use it to plot the result. (The remainder of the points are symmetric.)

Pyy = Y.*conj(Y)/251;
f = 1000/251*(0:127);
plot(f,Pyy(1:128))
title('Power spectral density')
xlabel('Frequency (Hz)')

%
% Zoom in and plot only up to 200 Hz. Notice the peaks at 50 Hz and 120 Hz. These are the frequencies of the original signal.

plot(f(1:50),Pyy(1:50))
title('Power spectral density')
xlabel('Frequency (Hz)')

%
%% See also: MATLAB Central's FFT Easy Demo
% MATLAB Central's FFT Easy: https://uk.mathworks.com/matlabcentral/fileexchange/71001-fft-easy-with-live-editor
%
%% Visualising the Fourier Transform
% Visualising the Fourier Transform Having computed a DFT it might be useful to visualise its result:
% It's useful to visualise the Fourier Transform
% Standard tools
% Easily plotted in MATLAB
% fourierspectraeg.m:
%
% Create an N-point cosine:

N=16;
Fs = 16;
Freq = 2;

x=cos(2*pi*Freq*(0:1:N-1)/N)';

figure
stem(0:N-1,x,'.');
axis([-0.2 N -1.2 1.2]);
legend('N-pt Sample Cosine signal x(n)');
xlabel('n \rightarrow');

%
% Take its Fourier transform and compute normalised Magnitude Spectrum:
%

X=abs(fft(x))/N;
stem(0:N-1,X,'.');
axis([-0.2 N -0.1 1.1]);
legend('Magnitude spectrum |X(k)|');
xlabel('k \rightarrow')

% Note: Axes are still labelled a the Npoints (same dimension as input signal)
%
% Let's relabel the x axes with frequencies, :

figure
X=abs(fft(x))/N;
f=((0:N-1)/N)*Fs;

stem(f,X,'.');
axis([0 (N-1)*Fs/N -0.1 1.1]);
legend('Magnitude spectrum |X(f)|');

xlabel('f in Hz \rightarrow')

%
% Let's relabel the x axes with frequencies, but now lets centre the FFT with fftshift() :

figure
X=abs(fftshift(fft(x)))/N;

stem(-N/2:N/2-1,X,'.');
axis([-N/2 N/2 -0.1 1.1]);
legend('Magnitude spectrum |X(f)|');
xlabel('f in Hz\rightarrow')

%
%
% Now let's increase the sample rate, and change the number of points, and frequency:

Fs=1024;
N=2*Fs;
Freq = 200;

x=cos(2*pi*(Freq/Fs)*(0:1:N-1))';
f=((0:N-1)/N)*Fs;
X=abs(fft(x))/N;
stem(f,X,'.');
axis([-0.2*Fs/N max(f) -0.1 1.1]);
legend('Magnitude spectrum |X(f)|');
xlabel('f in Hz \rightarrow')

%
% And centred:

figure
fshift=((-N/2:N/2-1)/N)*Fs;
X=abs(fftshift(fft(x)))/N;

stem(fshift,X,'.');
axis([-0.5*(N-1)*Fs/N 0.5*(N-1)*Fs/N -0.1 1.1]);
legend('Magnitude spectrum |X(f)|');
xlabel('f in Hz\rightarrow')



%
%% Magnitude Spectrum
% The Magnitude Spectrum of Fourier Transform (as shown above)
% Recall that the Fourier Transform of our real audio/image data is always complex
% Phasors: This is how we encode the phase of the underlying signal's Fourier Components.
% How can we visualise a complex data array?
% Back to Complex Numbers:
% Magnitude spectrum: Compute the absolute value of the complex data:
% where F_R(k) is the real part and F_I(k) is the imaginary part of the N sampled Fourier Transform, F(k).
%  Recall MATLAB: Sp = abs(fft(X,N))/N; (Normalised form)
%
%% Phase Spectrum
%
% The Fourier Transform also represent phase, the phase spectrum is given by:
%
% Recall MATLAB: phi = angle(fft(X,N))
%
%% Relating a Sample Point to a Frequency Point
%
% When plotting graphs of Fourier Spectra and doing other DFT processing we may wish to plot the x-axis in Hz (Frequency) rather than sample point number k =0,1,\ldots,N-1
%
% There is a simple relation between the two:
% The sample points go in steps k =0,1,\ldots,N-1
% For a given sample point k the frequency relating to this is given by: f_k = k\frac{f_s}{N} where f_s is the sampling frequency and N the number of samples.
% Thus we have equidistant frequency steps of \frac{f_s}{N} ranging from 0 Hz to \frac{N-1}{N}f_s Hz
% See the code in fourierspectraeg.m (above) for example plot:
% Fs=44100;
% f=((0:N-1)/N)*Fs;
% X=abs(fft(x,N))/N;
% stem(f,X);
% axis([-0.2*FS/N max(f) -0.1 1.1]);
% legend('Magnitude spectrum |X(f)|');
% xlabel('f in Hz \rightarrow')
%
%
%%
%
%% The Spectrogram
% A Spectrogram gives a time-Frequency representation of the data.
% It is often useful, for many forms of signal analyis, to look at the frequency distribution over a short-time:
% Split signal into N segments
% Do a windowed Fourier Transform — Short-Time Fourier Transform (STFT)
% Window needed to reduce leakage effect of doing a shorter sample SFFT.
% Apply a Blackman, Hamming or Hanning Window
% MATLAB function does the job:  Spectrogram — see doc/help spectrogram

showDoc = 0;
if showDoc
    doc spectrogram;
end

%

load('handel')
[N M] = size(y);
figure;
spectrogram(y,512,20,1024,Fs);

%
%
%% Another MATLAB spectrogram Example: Aphex Twin
% Previous slide use the free and excellent Sonic Visualiser.
% We of course know how to display the spectrogram image in MATLAB:

aphex = webread('http://users.cs.cf.ac.uk/Dave.Marshall/CM2208/MATLAB/Fourier_Transform/FormulaSnippet.wav');

% aphex = audioread('FormulaSnippet.wav'); % Local file version

mono = (aphex(:,1) + aphex(:,2))/2;
spectrogram(mono,1024,120,2048,'power','yaxis'); 
set(gca,'YScale','log'); 
colormap('winter'); 
xlabel('Time');
ylabel('Frequency (Log Scale)');

%
% Note: we change the display of the spectrogram to a log scale, which looks better.
% The audio clip used above is available here: FormulaSnippet.wav
%
% A MATLAB emulation of the basic idea:
% Use the Daphex.m  code  and  im2sound.m supporting code (listed below) ,
% To run code below you need to download this file:  Dave_Frame0001.jpg and put it the same folder as this script or change path or use the loader switch in code:

% Download image in current directory?
downloadIMInCurrentDirectory=1;
if downloadIMInCurrentDirectory
    Im = webread('http://users.cs.cf.ac.uk/Dave.Marshall/CM2208/MATLAB/Fourier_Transform/Dave_Frame0001.jpg');
    imwrite(Im,'Dave_Frame0001.jpg');
end

figure;

imshow(imread('Dave_Frame0001.jpg'));
dave_im2snd = im2sound('Dave_Frame0001', 'jpg', 44100, 40,6000,0.00002, 10);

sound(dave_im2snd,44100);

figure;

spectrogram(dave_im2snd,1024,120,2048,'power','yaxis');
set(gca,'YScale','log');

colormap('winter');

%
%% Properties of Fourier Transforms
% Here are just a few of the other properties of Fourier transforms which are useful when reasoning, or computing, with Fourier transforms:
% Linear Operator
% Shifting
% Scaling
% Rotation
% Zeroth component
% Convolution Theorem (see next section     applications)
%
%% Linear Operator Theorem
%
% The Fourier transform is a linear operator. This means that:
% If f(x) and g(x) are two functions with Fourier transforms F(u) and G(u), then the Fourier transform of a f(x)+b g(x) where a and b are
%     constants is simply a F(u)+b G(u).
%
%
%% Fourier Transform Linear Operator MATLAB Demo: FFT_Linear.m
%

% Fourier Transform Linear Operator Demo
% Create a white box on a black background image
M= 256; N = 256;
image1 = zeros(M,N);
box = ones(64,64);
%box at centre
image1(97:160,97:160) = box;

% Create a single white box in an image:


figure
imshow(image1)
title('Image 1: White Box 1');

% Create another  white box on a black background image (at different image location):


% Create another  white box on a black background image
M= 256; N = 256;
image2 = zeros(M,N);
box = ones(32,32);
%box at centre
image2(37:68,37:68) = box;


figure
imshow(image2)
title('Image 2: Another white box');

% Make a composite image by adding these two images together:

% Make composite image
image3 = image1 + image2;

figure
imshow(image3)
title('Image 3: Two white boxes (Image 1 + Image 2)');

% Now compute the Fourier transfgorm of the respective images

% Compute Fourier Transforms.

imft1 = fft2(double(image1));
imft2 = fft2(double(image2));
imft3 = fft2(double(image3));

% Lets add the Fourier Spectra of Image 1 + Image 2 together and display the Spectra:

figure
imagesc(abs(fftshift(imft1) + fftshift(imft2)))
title('Image 1 + Image 2 Spectra');

% As the linear operator theorem predicts:  the Image 1 spectra + Image 2  spectra is the same as the spectra of Image 1 + Image 2 (Image 3):

figure
imagesc(abs(fftshift(imft3)))
title('imft3 Spectra');

% To prove this there is no residual if we subtrat the Image 1 spectra + Image 2  spectra form Image 3's spectra:

figure
imagesc(abs(fftshift(imft1) + fftshift(imft2) - fftshift(imft3)) > eps('single'))


title('Difference (Image 1 + Image 2) - imft3 Spectra');

% Also if we Fourier transforminvert the Image 1 spectra + Image 2  spectra we get Image 3 as expected:

figure(7)
imshow(ifft2(imft1 + imft2));
title('Inverse FT of imft1 + imft2');

%
%% Shifting Operator Theorem
%
% Shifting the real space data through a fixed distance has the effect that:
% If f(x) is a function with Fourier transform F(u), then the Fourier transform of f(x-x_0) is given by e^{-2\pi i x_0 u}F(u).
%
% Note: that the exponential term here has unit magnitude for all values of u.
%                        Thus, the magnitude of the resulting Fourier transform is left unchanged, although the relative sizes of the real and imaginary parts are altered.
%
%% Fourier Transform Shifting Operator MATLAB Demo: FFT_Shift.m
% Create a white box on a black background image:

% Fourier Transform Shifting Operator Demo

% Define Initial box
x0 = 97; 
y0 = 97;
xwidth = 64;
ywidth = 64 ;

% Create a white box on a black background image
M= 256; N = 256;
image1 = zeros(M,N);
box = ones(xwidth,ywidth);
%box at centre
image1(x0:x0+xwidth-1,y0:y0+ywidth-1) = box;

figure
imshow(image1)
title('Image 1');

% Now let's shift the box (You could experiment with different values for dx and dy)

% Define shift
dx= 40; 
dy= 0;


% Shift Image 1
image2 = zeros(M,N);
image2(x0+dx:x0+dx+xwidth-1,y0+dy:y0+dy+ywidth-1) = box;

figure
imshow(image2)
title('Image 2: Shifted box by (dx,dy)');

% Now let's compute the Fourier transfoms and display Spectra:

% Compute  Fourier Transfoms.

imft1 = fft2(double(image1));
imft2 = fft2(double(image2));


figure
imagesc(abs(fftshift(imft1)));
title('FT of Image 1 (imft1)');
figure
imagesc(abs(fftshift(imft2)));
title('FT of Image 2 (Shofted imft1, imft2)');

% Now let's apply the shifty in the frequency domain, we need to compute and apply: e^{-2\pi i x_0 u}F(u) in 2D to FT of image 1:

% Define shift in frequency domain
[yF,xF] = meshgrid(-M/2:M/2-1,-N/2:N/2-1);

% Perform the shift
imft3=imft1.*exp(-1i*2*pi.*(xF*dx+yF*dy)/256);

figure
imagesc(abs(fftshift(imft3)));
title('FT of imft3 (imft1 shifted in Fourier Space');

% The inverse FT of FT image, imft3, is:

figure
imshow(abs(ifft2(imft3)));
title('Inverse FT of imft3');

% As expected this is the same as image 2, the original image shifted in real space.


figure
imagesc(abs(ifft2(imft3)-image2) > eps('single'));
title('Residual of IFT of imft3 and image2');

%
%
%% Scaling Theroem
% If we scale the spacing of the real space data in distance, we have that:
% If f(x) is a function with Fourier transform F(u), then the Fourier transform of f(a x) where a is a real constant is given by \frac{1}{|a|}F(\frac{u}{a}).
%
% In other words:
% if we spread out the data in real space, the data is compressed into more compact region of Fourier space.
% This is intuitive — if we double the spacing of a grid pattern, its spatial frequency is halved.
% Note that the magnitude of the Fourier representation is also affected.
%
%% Fourier Transform Scaling Operator MATLAB Demo: FFT_Scale.m
% Create nother white box to start with:

% Fourier Transform Rotate Operator Demo

% Define Initial box

xwidth = 40;
ywidth = 40;

% Create a white box on a black background image
M= 256; N = 256;

x0 = M/2 - xwidth/2;
y0 = N/2 - ywidth/2;


image1 = zeros(M,N);
box = ones(xwidth,ywidth);
%box at centre
image1(x0:x0+xwidth-1,y0:y0+ywidth-1) = box;

figure
imshow(image1);
title('Image 1: A white box');

% Now let's scale this box, we define a scaling in x and y: (scalex, scaley) – expeiment with different values.

% Define scale
scalex = 0.5;
scaley = 0.5;


% Image 2 = Scaled Image 1
image2 = zeros(M,N);
[px,qx] = rat(scalex);
[py,qy] = rat(scaley);

in_image2 = resample(resample(double(image1),px,qx)',py,qy)';

xrange = int16(M/2-M*scalex/2:M/2+M*scalex/2-1);
yrange = int16(N/2-N*scaley/2:N/2+N*scaley/2-1);
image2(xrange,yrange)= in_image2; 

figure
imshow(image2);
title('Image 2: Scaled Image 1');

% Now compute Fourier transforms and display, note the scaling:

% Compute  Fourier Transfoms.
imft1 = fft2(double(image1));
imft2 = fft2(double(image2));

figure
imagesc(abs(fftshift(imft1)))
title('Fourier Transform of Image 1 (imft1)');
figure
imagesc(abs(fftshift(imft2)))
title('Fourier Transform of Image 2 (imft2)');

%
%% 2D Fourier Transform: Rotation
%
% One useful property in two dimensions is that:
% if we rotate the real space data, its Fourier transform is rotated by the same angle, or more exactly:
%                         If f(x,y) is a function with Fourier transform F(u,v), on expressing these functions in terms of polar coordinates r, \theta, \rho, \phi where x=r\cos\theta, y=r\sin\theta, u=\rho\cos\phi, v=\rho\sin\phi so that f(x,y) and F(u,v) become f(r,\theta) and F(\rho,\phi) respectively, the Fourier transform of f(r,\theta+\omega) where \omega is a constant is given by F(\rho,\phi+\omega).
%
% Basically,  this means the two-dimensional Fourier transform is an intrinsic property of the data which is independent of our choice of axis directions.
%
%% Fourier Transform Rotate Operator Demo: FFT_Rotate.m
% Create a box again:

% Fourier Transform Rotate Operator Demo

% Define Initial box

xwidth = 64;
ywidth = 32;

% Create a white box on a black background image
M= 256; N = 256;
image1 = zeros(M,N);

x0 = M/2 - xwidth/2;
y0 = N/2 - ywidth/2;

box = ones(xwidth,ywidth);
%box at centre
image1(x0:x0+xwidth-1,y0:y0+ywidth-1) = box;

figure
imshow(image1);
title('Image 1');

% Rotate the image, defined (in degress) by rot:

% Define Rotation in degrees
rot = 45; 
image2 = imrotate(image1, rot,'bilinear','crop');

figure
imshow(image2);
title('Image 2: Image 1 rotated');

% Compute the Fourier transforms and display, note the rotated spectra:

% Compute  Fourier Transfoms.

imft1 = fftshift(fft2(double(image1)));
imft2 = fftshift(fft2(double(image2)));


figure
imagesc(abs(imft1))
title('Fourier Transform of Image 1 (imft1)');

figure
imagesc(abs(imft2))
title('Fourier Transform of Image 2 (imft2)');

%
%
%
%
%
%
%
%
%% Fourier Transform: Zeroth Component
%
% A final property is that the zeroth component of the Fourier space representation is just the average data value (apart from a factor of 1/N, in two dimensions – assuming balanced FT)).
% This is illustrated below for the two-dimensional DFT case:
% F(u,v) = \frac{1}{N}\sum_{x=0}^{N-1}\sum_{y=0}^{N-1} f(x,y) e^{-2\pi i (x u+y v)/N},
%  so F(0,0) = \frac{1}{N}\sum_{x=0}^{N-1}\sum_{y=0}^{N-1} f(x,y) = N \overline f(x,y).
%
%% MATLAB Illustration of zeroth element
%
% So F(1) (1D), or F(1,1) (2D), is  the mean of the data plus a factor N: FFT_zeroth.m

% Illustration of zeroth element, F(1) (1D) or F(1,1) (2D), being the mean of the data plus a factor N

% 1D Example
y = 1:100
yft = fft(y);
yft(1)/100  % Zeroth Element of y
mean(y) % Mean of Data

% 2D Example
im = imread('cameraman.tif');

[N, M] = size(im);
imft = fft2(im);

% The zeroth element is:

imft(1,1)/(N*M) % Zeroth Element of imft

% The mean of the actual image is:

mean(mean(im)) % Mean of Data

% Note:  The FT is not balanced in MATLAB so we divide by N*M
%
%
%% Centring the Frequency of a Fourier Transform
% Most computations of FFT represent the frequency from 0 --- (N-1) samples (similarly in 2D, 3D etc.) with corresponding frequencies ordered accordingly — the 0  frequency is not really the centre.
% We frequently like to visualise the FFT as the centre of the spectrum.
% In 1D (Audio/Vector): swaps the left and right halves of the vector
% Similarly in 2D (Image/Matrix) we swap the first quadrant with the third and the second quadrant with the fourth:
%
% This is possible due the invariant shift property of the Fourier Transform.

% See doc fftshift for more information.
% Also see lab sheet exercises.

showDoc = 0;
if showDoc
    doc fftshift;
end

%
%% Supporting functions
% mysin():

function s = mysin(amp,freq,phase,dur, samp_freq)
% example function to so show how amplitude,frequency and phase
% are changed in a sin function
% Inputs: amp - amplitude of the wave
%         freq - frequency of the wave
%         phase - phase of the wave in degree
%         dur - duration in number of samples
%         samp_freq - sample frequncy


x = 0:dur-1;
phase = phase*pi/180;

s = amp*sin(2*pi*x*freq/samp_freq + phase);

end

% im2sound.m :

function [final_sound] = im2sound(filename, ext, f_sample, f_low,f_high, amp_mod, sample_t)

%INPUTS:
%'filename' - Name of the image to be encoded (not including extension
%ext' - Extension of the image (not including "." at the beginning).  
%'f_sample' - Sampling frequency (Hz)
%'f_low' - Lowest frequency (Hz) (e.g. 40)
%'f_high' - Highest frequency (Hz) (e.g. 6000)
%'amp_mod' - Multiplication factor for the amplitude.  Decrease until 
%image is clear. Too high and the waveform clips.  Too low and the image 
%is very dark (e.g. 0.00002)
%'sample_t' - Duration of the sample in seconds.  Longer samples have
%better quality (e.g. 10)

%OUTPUTS:
%'final_sound' - the final sound containing the image.  This is
%automatically saved to a .wav file with the original image filename


%INITIALISING VARIABLES:
%The waveform at each time point.  This is reset at the beginning of each
%time point
temp_sound = 0; 
%The final waveform
final_sound = 0; 


%MAIN BODY
%Loading the sample image and calculating the image size
raw_im = imread(strcat(filename,'.',ext));
size_raw_im = size(raw_im);

%Making a frequency table for the height of the image.  Each row of the
%image is assigned a particular frequency from the corresponding row of 
%this table.  The frequencies are linearly distributed between the highest and
%lowest user-definied frequencies.  "f_step" is the increment between each
%adjacent frequency
f_step = (f_high - f_low)/size_raw_im(1,1);
f_table = (f_high:-f_step:f_low);

%The final sound will dwell on each column of the image for a specific
%time.  This time is defined by "t_start" and "t_end".  It depends on how
%long the user determined the sound-clip should be and how wide (how many
%columns) the image is.  
t_step = (sample_t/size_raw_im(1,2));

%Initial values for the start and end times.  These will be increased at
%the end of each loop iteration (when the script moves onto the next column
%of the image).
t_start = 0;
t_end = t_step;

%The loop which generates the sound file.  At each iteration it generates a
%segment of the final sound file, which is temporarily saved to 
%"temp_sound".  This segment is built up of frequencies from that
%particular column of the image.
for j = 1:size_raw_im(1,2)
    %Initialising the variable (the sound for each frequency (row) is added
    %to the existing sound)
    temp_sound = 0;
    
    %Setting the time in matrix format
    t = t_start:1/f_sample:(t_end);
    
    %For each iteration of this loop, the script goes down the current
    %column of the image and generates a waveform of the frequency
    %specified in "f_table".  The amplitude of the waveform is determined
    %by the pixel intensity.  This generated waveform is added to all the
    %previously generated waveforms in that particular column
    for i = size_raw_im(1,1):-1:1
        temp_sound = temp_sound+ sin(2*pi*t*f_table(i))*...
            double(raw_im(i,j))*amp_mod;
        
    end
       
    %At the end of each column the segment of sound generated is added to
    %the end of the existing sound file ("final_sound").
    final_sound = cat(2,final_sound,temp_sound);
    
    %The temporary sound is cleared ready for the start of the next column
    clear temp_sound
    
    %Moving to the next time frame
    t_start = t_start + t_step;
    t_end = t_end + t_step;
    
end

%This saves "final_sound" to the '.wav' file of the same name as the input
%fil

audiowrite(strcat(filename, '.wav'), final_sound, f_sample);

end