# program to demontrate filtering. We desire a low pass filter so we first
# create the ideal freq response being sure to symmetrize. The filter
# is the obtained by inverse fft-ing and retaining only the most significant
# coeffs (just to speed up calculations).
# we then run the result on white noise and observe the spectrum of the
# filtered data.
# then we filter actual audio
#
N = 1024 # number points in fft
cutoff = 100 # cutoff freqs higher than this (in cycles/fftlen)
retain = 40 # will keep this many twice this many filter coeffs
A = rep(0,N) # A will be the idealized filter response
ran = 2:cutoff # the freqs we want to "pass"
A[ran] = 1
A[N-ran+2] = 1 # symmetrize the freq response
a = Re(fft(A,inverse=TRUE)) # the inverse fft --- our ideal filter
#AA = fft(a) # a check to see if it has the right freq response
#plot(Mod(AA[1:(N/2)]))
#stopp
a = c(a[(N-retain+2):N] , a[1:retain])/N # take the center 2*retain coeffs of filter
#plot(a,type='b');
#stopp
#AA = fft(c(a,rep(0,N-length(a)))) # pad to have length N and FFT
#plot(Mod(AA[1:(N/2)]),type='l') # filter response is close to what we want
#stopp;
# try it on noise
y = runif(N,min=-.5,max=.5) # white noise
z = filter(y,a,method="conv",circular=T) # noise filtered by our filter a
Z = fft(z)
plot(Mod(Z[1:(N/2)])) # note that the desired characteristics are achieved
# only low freqs appear in filtered signal
#stopp;
# try it on a real sound
library(tuneR)
w = readWave("sound/bass_oboe.wav") # take actual sound data from file
play(w,"play") # listen to a bit
bits = w@bit
sr = w@samp.rate
y = w@left # just work with left channel
y = y[1:160000]/max(y) # take 20 secs of sound
y = filter(y,a,method="conv",circular=T) # filter the data with our lowpass filter
y = as.vector(y) # obnoxious kludge: Wave only deals with "vector" data
u = Wave(round(2^(bits-1)*y), samp.rate = sr, bit=bits) # make wave struct
play(u,"play")