# Program to identify the freqs amps and phases of a function
# that is a sum of sines each of which oscillates an integer number
# of times for N points
N = 1024 # FFT works best with powers of 2
t = seq(from=0,by=2*pi/N,length=N) # N evenly spaced pts 0 -- 2*pi
f = 100 # number of cycles per N pts
ph = 1 # the phase at f
amp = 3 # the amplitude at f
#y = amp*cos(f*t+ph); # our "signal"
y = 1*cos(f*t+ph) +
2*sin(2*f*t + 10) +
3*cos(3*f*t -1);
Y = fft(y) # take the fft
plot(Mod(Y[1:(N/2)]) / (N/2)) # plot of Modulus = amplitudes
# Note how only one (at f) is non-zero
print("amp at f is ")
print(Mod(Y[f+1])/(N/2)) # remember sine oscillating f times
# appears at f+1
# also must divide by N/2 to get correct amp
print("phase at f is")
print(Arg(Y[f+1])) # the phase