clear
close all

% We can also use Fourier transforms to clean up signals.
t = .01:.01:10;

% Lets make a function with two frequencies and some noise
y1 = sin(2*pi*t);
y2 = sin(2*pi*3*t + .1);
y = y1 + y2 + 3*(rand(length(t),1)')-1.5;
figure
plot(t,y)

% We can take the FFT in exaclty the same way
Y = fft(y);
Yabs = abs(Y);
f = (100/2)*linspace(0,1,length(y)/2);

figure
plot(f,Yabs(1:500))
xlabel('Frequency (Hz)')

% If we keep only frequencies less than 5 Hz, we should keep all of our
% signal.  I will make a mask of 1s and 0s to keep only the frequencies we
% want.
Fcutoff = find(f > 5);
mask = ones(500,1)';
mask(Fcutoff) = 0;

figure
plot(f,Yabs(1:500),f,200*mask)

% But you have to remember that the full FFT is 2x as big as what we're
% looking at.  So we have to apply our filter to that also.  I use a fliplr
% to generate the full mask.
fullmask = [mask fliplr(mask)];

figure
plot(Yabs)
hold all
plot(200*fullmask)
hold off

% So now we can multiply our FFT by the mask, and then take the IFFT to get
% our filtered signal back.  We have to take the real component because the
% transform back is not perfect, so there are some imaginary numbers.
Yfilt = Y.*fullmask;
yfilt = real(ifft(Yfilt));

figure
subplot(2,1,1)
plot(t,y,t,yfilt,'linewidth',1.5)
set(gco,'linewidth',4)
subplot(2,1,2)
tP = 101:150;
plot(t(tP),y(tP),t(tP),yfilt(tP))


% Now what if we make a high pass filter
Fcutoff = find(f < 5);
mask = ones(500,1)';
mask(Fcutoff) = 0;
fullmask = [mask fliplr(mask)];

figure
plot(f,Yabs(1:500),f,200*mask)

% Now we perform the same ifft
Yfilt = Y.*fullmask;
yfilt = real(ifft(Yfilt));

figure
subplot(2,1,1)
plot(t,y,t,yfilt,'linewidth',1.5)
set(gco,'linewidth',4)
subplot(2,1,2)
tP = 101:190;
plot(t(tP),y(tP),t(tP),yfilt(tP))

mo = y- yfilt;
figure
plot(mo)
hold all
plot(y)

%% 2D Fourier transform
clear
close all

img = imread('Football.jpg');

figure
imshow(img)

img2 = double(rgb2gray(img));
figure
imshow(img2,[])

% We can use the fft2 command to take a 2D fourier transform
IMG = fft2(img2 - mean(img2(:)));
IMG(1,1) = 0;
size(img2)
size(IMG)

% We need to look at the absolute value again, since there are imaginary
% phase components in the 2D spectrum

figure
imshow(log(abs(IMG)),[])

% When we take the 2D FFT, we need to run the circshift command to put the
% DC component in the center of the image
IMG2 = fftshift(IMG);
figure
subplot(1,2,1)
imshow(img2,[])
subplot(1,2,2)
imshow(log(abs(IMG2)),[])
