% 1 the fft function
% implements the fast Fourier transform
% first we see how to use it
x = rand(1,10);
y = fft(x);
iy = ifft(y);
x2 = real(iy);
norm(x-x2);
norm(x-x2)

ans =

  3.6926e-016

% 2. computing the spectrum of a signal
% let us build a signal
dt = 1/100;
et = 4;
t = 0:dt:et;
y = 3*sin(4*2*pi*t) + 5*sin(2*2*pi*t);
subplot(2,1,1)
plot(t,y); grid on
axis([0 et -8 8])
xlabel('Time (s)')
ylabel('Amplitude')
% the fft transforms the time into
% the frequency domain
Y = fft(y);
n = size(y,2)/2;
amp_spec = abs(Y)/n;
subplot(2,1,2)
freq = (0:79)/(2*n*dt);
plot(freq,amp_spec(1:80)); grid on
xlabel('frequency (Hz)')
ylabel('Amplitude')
% 3. Filtering noise from signals
noise = randn(1,size(y,2));
ey = y + noise;
figure
subplot(2,1,1)
plot(t,y)
plot(t,ey)
% to remove the noise, we compute first
% the frequency spectrum of (t,ey)
eY = fft(ey);
n = size(ey,2)/2;
amp_spec = abs(eY)/2;
xlabel('Time (s)')
ylabel('Amplitude')
subplot(2,1,2)
xlabel('Frequency (Hz)')
ylabel('Amplitude')
freq = (0:79)/(2*n*dt);
plot(freq,amp_spec(1:80)); grid on
% forgot to scale the amplitudes
amp_spec = abs(eY)/n;
plot(freq,amp_spec(1:80)); grid on
% we observe the noise to be of low
% amplitude, we recognize the spikes
% for the real important frequencies
fY = fix(eY/100)*100;
% we removed all Fourier data less than 100
ifY = ifft(fY);
cy = real(ifY);
figure
subplot(2,1,1);
plot(t,y);
subplot(2,1,2);
plot(t,cy);
diary off