% Practice run with FFT

Fs = 1000;      % Sampling frequency
T  = 1/Fs;      % Sample time
L  = 1000;      % Length of signal in units of T
t  = (0:L-1)*T; % Time vector

% Time domain
y = 0.7*sin(2*pi*50*t) + 0.2*sin(2*pi*440*t);
% y = zeros(1,L);
subplot(221);
plot(t,y);
xlim([0,0.5]);
xlabel('time (milliseconds)');

% Frequency domain
subplot(222);
NFFT = 2^nextpow2(L);
Y = fft(y,NFFT)/L;
f = Fs/2*linspace(0,1,NFFT/2);

plot(f,2*abs(Y(1:NFFT/2)));

% Filter
F = zeros(1,length(Y));
a1 = floor(435/(Fs/NFFT));
a2 = floor(445/(Fs/NFFT));
F(a1:a2) = 1;
for i = NFFT/2+1:NFFT-1
    F(i) = F(NFFT-i);
end
Y = Y.*F;

subplot(223);
plot(f,2*abs(Y(1:NFFT/2)));

% inverse transform 
z = L*ifft(Y,NFFT);
subplot(224);
t = (0:length(z)-1)*T;
plot(t,real(z));
xlim([0,1]);
hold on;
plot(t,imag(z),'r');