Interpolation does not generally yield the ideal result. One way of looking at interpolation errors is with Fourier spectra.
In Fourier analysis one regards signals as superpositions of plane waves (i.e. sines and cosines). Given the sampling rate and the length of the signal, we can only compute the amplitudes of a finite sequence of frequencies (and we constructed those frequencies in class for an audio signal).
Consider the following function and arrays:
In [1]:
import numpy as np
def func(val):
return np.sin(4*np.pi*val) + np.cos(6*np.pi*val)
samples = np.linspace(0, 1, 17, endpoint = False)
x = np.linspace(0, 1, 51, endpoint = False)
Compute the value of this function at every point in samples.
Create an interpolator object (using the class we defined in our lectures) from samples and the values you have obtained.
Afterwards, create 2 more arrays.
One which will contain the interpolated values at the points in x, and one which will contain the exact values.
Compute the Fourier transforms of these 3 arrays (the coarse grained version, the interpolated version, and the exact fine-grained version).
Construct the frequencies corresponding to samples and x, and then plot the power spectra for the 3 arrays (i.e. $|\hat{f}(k)|^2$ for the various $\hat{f}$ you get).
Use logarithmic scale on the $y$ axis.
Discuss the plot.
5 points for creating the interpolator object correctly
5 points for creating the arrays with interpolated values and exact values
5 points for the 3 Fourier transforms
5 points for constructing the frequencies.
5 points for saying that, additionally to there being much larger errors in the off-peak frequencies, there are also two fake peaks in the spectrum
In [2]:
%matplotlib nbagg
import matplotlib.pyplot as plt
class interpolator:
# members:
# t -- array of function arguments where we know the value of the function
# x -- array of corresponding values of the function
# methods:
# interpolate -- linear interpolation
def __init__(self, t, x):
self.t = t
self.x = x
return None
def interpolate(self, tt):
for i in range(self.t.shape[0]):
if tt < self.t[i]:
i0 = i
break
else:
i0 = i
y = (self.x[i0] +
((self.x[i0] - self.x[i0-1]) /
(self.t[i0] - self.t[i0-1])) * (tt - self.t[i0]))
return y
fvalues = func(samples)
i = interpolator(samples, fvalues)
tfine = np.linspace(0, 1, 51, endpoint = False)
yfine = np.zeros(tfine.shape[0], tfine.dtype)
for ii in range(tfine.shape[0]):
yfine[ii] = i.interpolate(tfine[ii])
ax = plt.figure(figsize=(6, 3)).add_subplot(111)
ax.plot(samples, fvalues, marker = '.', label = 'coarse data')
ax.plot(tfine, yfine, marker = '+', linestyle = 'none', label = 'interpolated data')
ax.legend(loc = 'best')
Out[2]:
In [3]:
sfull0 = np.fft.rfft(fvalues)
sfull1 = np.fft.rfft(yfine)
sfull2 = np.fft.rfft(func(tfine))
k0 = np.array(range(sfull0.shape[0])).astype(samples.dtype) / samples[-1]
k1 = np.array(range(sfull1.shape[0])).astype(tfine.dtype) / tfine[-1]
ax = plt.figure().add_subplot(111)
ax.plot(k0, np.abs(sfull0)**2, label = 'orig data')
ax.plot(k1, np.abs(sfull1)**2, label = 'interpolated')
ax.plot(k1, np.abs(sfull2)**2, label = 'exact refinement')
ax.set_yscale('log')
ax.legend(loc = 'best')
Out[3]:
Follow a similar procedure for different numbers of samples: 15, 18, 21, ..., 30. Compare the spectra you obtain by interpolating on the differently sampled data.
10 points for doing this correctly
5 points for pointing out that the erroneous peaks move to the right, and they get less and less prominent.
In [4]:
def get_spec(nsamples):
xtmp = np.linspace(0, 1, nsamples, endpoint = False)
i = interpolator(xtmp, func(xtmp))
ftmp = np.zeros(tfine.shape[0], tfine.dtype)
for ii in range(tfine.shape[0]):
ftmp[ii] = i.interpolate(tfine[ii])
sfull = np.fft.rfft(ftmp)
return np.abs(sfull)**2
spec_list = []
nlist = range(15, 32, 3)
for n in nlist:
spec_list.append(get_spec(n))
ax = plt.figure().add_subplot(111)
for counter in range(len(nlist)):
ax.plot(k1,
spec_list[counter],
label = '${0}$'.format(nlist[counter]))
ax.set_yscale('log')
ax.legend(loc = 'best')
Out[4]: