% Lecture 35: Signal processing via FFT

x = rand(1,20);
y = fft(x);
iy = ifft(y);
x2 = real(iy);
norm(x2-x)

ans =

  4.7510e-016

% example of a signal 
et = 4;
dt = 0.01;
t = 0:dt:et;
y = 3*sin(4*pi*t)+4*cos(12*pi*t);
plot(t,y);
Y = fft(y);
n = size(y,2)/2;
spec = abs(Y)/n;
figure
subplot(2,1,1);
plot(t,y);
subplot(2,1,2);
freq = (0:79)/(2*n*dt);
plot(freq,spec(1:80));

% adding noise
noise = randn(1,size(y,2));
ey = y + noise;
eY = fft(ey);
espec = abs(eY)/n;
figure
subplot(2,1,1);
plot(t,ey);
subplot(2,1,2);
plot(freq,espec(1:80));


% clean up the spectrum and recover the original signal.
fY = fix(eY/100)*100;
figure
fspec = abs(fY)/n;
plot(freq, fspec(1:80));
ifY = ifft(fY);
cy = real(ifY);

% compare the original, perturbed and reconstructed signals
figure
plot(t,y);
hold on;
plot(t,ey,'g');
plot(t,cy,'r');
diary off