# this program demonstrates how the fft can be used to construct an arbitrary
# function as a sum of cosine waves. we first take any function and then
# plot out better and better approxs to it using more and more sine waves.
# when the user exits from the loop by hitting any other than return,
# program plays both the original sound and the approx to it.
bits =16
N = 256 # number of points in function (usually power of 2 for fft)
#y = c(c(seq(0,1,length=N/4),seq(1,-1,length=N/2)),seq(-1,0,length=N/4)) # my sawtooth
y = cumsum(rnorm(N)) # a more exotic possibility
plot(y,type='l') # plot the function
Y = 2*fft(y)/N # the scaled fft of function
Y[1] = Y[1]/2 # some strange normalization needed on 1st and last components
# will not use these much in practice
Y[N/2+1] = Y[N/2+1]/2
yhat = rep(0,N) # our approximation using cosine waves ... start at 0
t = seq(from=0,by=2*pi/N,length=N) # time points for cosine waves
for (i in 1:(1+N/2)) { # loop through fft components
yhat = yhat + Mod(Y[i])*cos((i-1)*t + Arg(Y[i])) # add in new cosine with amp and phase given
# by fft. ith fft value oscillates i-1
# times per N pts.
error = yhat-y # the error at present
print(sum(error*error)) # sum of squared errors (always decreasing)
plot(y,type='l') # plot the true function and our approximation
lines(yhat,lty=2)
c = scan(what="") # wait for keyboard input
if (length(c) > 0) break # if get a character, break, ow continue
}
# now play out the sound of the true function followed by approx
scale = max (abs(max(y)) , abs(min(y)))
y = y / scale
scale = max (abs(max(yhat)) , abs(min(yhat)))
yhat = yhat / scale
periods = 500
library(tuneR) # need to like with the sound library
sr = 48000 # sampling rate
s = c(rep(y,periods),rep(yhat,periods)) # replicate approx periods times
u = Wave(round(2^(bits-5)*s), samp.rate = sr, bit=16) # make wave struct
play(u,"play") # play it