% A basic demonstration of the FFT for translation from 
% time to frequency domains, with scaling for an arbitrary 
% time step and point set.  H. Motteler, May 1999.
%
% A typical FFT application has three parameters: period, 
% size of the time step, and number of points, two of which
% are independent; i.e., specifying any two fixes the third.
% In this example we specify time step and number of points,
% solve for period, and scale the frequency plot accordingly.
% This would be typical for the case where the sample time
% step was given and we wanted to do the transform with 2^k
% points.
%
% For this demo, you might want to start with a time step 
% of 1/64 (0.015625), 128 points, and any two frequencies 
% below Nyquist.

dt = input('time step > ');
fprintf(1, 'we need at least %g points\n', 2/dt)
np = input('number of points > ');

tt = dt * (np-1); % period is dt*np
tseq = 0:dt:tt;   % time step grid

fmax = 1 / (2*dt);  % max representable frequency
fprintf(1, 'nyquist frequency is %g\n', fmax);

sf1 = input('signal #1 frequency > ');
sf2 = input('signal #2 frequency > ');

sp1 = 1/sf1;
fprintf(1, 'signal #1 period is %g\n', sp1);
sp2 = 1/sf2;
fprintf(1, 'signal #2 period is %g\n', sp2);

% generate a 2-component signal
y1 = cos(2*pi*tseq/sp1);
y2 = cos(2*pi*tseq/sp2);
y = y1 + y2;

% plot time-domain signal y, with unit tics
tx = 0:1:tt;
ty = zeros(1,length(tx));

subplot(2,1,1)
plot(tx, ty, '+', tseq, y)
v = axis;
axis([0, tseq(length(tseq)), v(3), v(4)])
xlabel('Time')
grid

% calculate and plot the FFT of y
fy = fft(y);
fx = (0:np-1) / (dt*np); % scale X axis

subplot(2,1,2)
plot(fx, real(fy), fx, imag(fy));
v = axis;
axis([0, fx(length(fx)), v(3), v(4)])
legend('real', 'imag.', 4)
xlabel('Frequency')
grid