```
import numpy as np
import scipy as sp
import scipy.signal
import matplotlib
import matplotlib.pyplot as plt
%matplotlib inline
```

The linear interpolation is a convolution of the samples with a triangular pulse $h_l(t)$. The triangular pulse is the convolution of two rectangular pulses $ h_l(t) = (p_{\tau} \ast p_{\tau}) (t) $. The frequency response $ H_l(f) $ is therefore 
\begin{equation}
H_l(f) = (\tau \cdot \text{sinc}(f \tau))^2.
\end{equation}

For the linear interpolator $ \tau = T_S $ where $ T_S $ is the sampling time of the input discrete signal.


```
# sampling period & frequency ofthe input discrete signal
Ts = 1
Fs = 1/Ts

# analytic frequency response of the interpolating filter
N = 40
res = 100
f = np.linspace(-N*Fs,N*Fs, 2*N*res)
tau = Ts
H_l = tau**2 * np.sinc(f * tau)**2
```


```
fig, ax = plt.subplots(1, 1, figsize=(13,4))
ax.plot(f, 20 * np.log10( np.abs(H_l) ))
ax.set_ylabel('Amplitude [dB]', fontsize=12)
ax.set_xlabel('Frequency [Hz]', fontsize=12)

n = np.arange(0,9,1)
ax.set_xticks(Fs * n)
labels = ["$%d F_S$" % y for y in n[1:]]
labels.insert(0,"0")
ax.set_xticklabels(labels, fontsize=12)
ax.grid()

ax.axis('tight');
ax.set_ylim([-100, 10]); # semicolon suppresses the output
ax.set_xlim([0, n[-1]*Fs]);
```

With this operation I get a continuous signal from a discrete one. If I want to output a discrete signal, I must sample the output. This is the same as both:
 - sampling the impulse response of the linear interpolation filter, 
 - interpolating zeros to the input signal

with sampling time 
\begin{equation} T_{S1} = L \cdot T_S \qquad L \in \mathbb{N}. \end{equation}

This produces spectral copies of the original frequency response, at $ n F_{S1} $



```
L = 6 # Fs1 is L times Fs
Fs1 = L * Fs

# view
f_min = 0
f_max = 8*Fs

fig, ax = plt.subplots(1, 1, figsize=(12,6))

# plot some spectral copies
f = np.linspace(f_min, f_max, ((f_max-f_min)/Fs)*res)
tau = Ts
H_l_tot = np.zeros((len(f),))
for i in range(-10,10):
    H_l = tau**2 * np.sinc((f + i*Fs1) * tau)**2
    ax.plot(f, 20 * np.log10( np.abs(H_l) ), '--b', alpha=0.4)
    H_l_tot = H_l_tot + H_l

# plot the sum of all the spectral copies
ax.plot(f, 20 * np.log10( np.abs(H_l_tot)), 'b')

ax.set_ylabel('Amplitude [dB]', fontsize=12)
ax.set_xlabel('Frequency [Hz]', fontsize=12)

n = np.arange(0,9,1)
ax.set_xticks(Fs * n)
labels = ["$%d F_S$" % y for y in n[1:]]
labels.insert(0,"0")
ax.set_xticklabels(labels, fontsize=12)
ax.grid()

ax.axis('tight');
ax.set_ylim([-100, 10]); # semicolon suppresses the output
ax.set_xlim([f_min, f_max]);
```

We show that by fourier transforming the impulse response of the linear interpolator, sampled at $F_{S1}$, we obtain the same result.


```
nfft = 1024
h_l_sampled = np.hstack((np.linspace(0,1-1/L,L), np.linspace(1,0+1/L,L)))
H_l_sampled = np.fft.rfft(h_l_sampled, nfft, axis=-1)
f_sampled = Fs1 * np.linspace(0,0.5,(nfft/2)+1)
```


```
fig, ax = plt.subplots(1, 1, figsize=(6,4))

#plt.title('Digital filter frequency response')

# highest frequency (pi) corresponds to M * Fs/2
ax.plot(f_sampled, 20 * np.log10(abs(H_l_sampled)), 'b',  label=r"$ALIAS_{10}$")
ax.set_ylabel('Amplitude [dB]', fontsize=12)
ax.set_xlabel('Frequency [Hz]', fontsize=12)
ax.grid()

n = np.arange(0,9,1)
ax.set_xticks(Fs * n)
labels = ["$%d F_S$" % y for y in n[1:]]
labels.insert(0,"0")
ax.set_xticklabels(labels, fontsize=12)
ax.set_xlim([0, L*Fs/2]);
```

Lets normalize them to compare the frequency response of the sampled version with respect to the continuous version


```
# frequency response of the sampled (at L*Fs) interpolator 
L = 6
Fs1 = L * Fs

# fft computation
nfft = 1024
h_l_sampled = np.hstack((np.linspace(0,1-1/L,L), np.linspace(1,0+1/L,L)))
H_l_sampled = np.fft.rfft(h_l_sampled, nfft, axis=-1)
f_sampled = Fs1 * np.linspace(0,0.5,(nfft/2)+1) # [0, .5 Fs1]

#Normalization
H_l_sampled = H_l_sampled/np.max(np.abs(H_l_sampled))

# frequency response of the continuous-time interpolator
f = np.linspace(0*Fs,(L/2)*Fs, (L/2)*100) # we care only up to (L/2)*Fs
tau = Ts
H_l = tau**2 * np.sinc(f * tau)**2

#Normalization
H_l = H_l / np.max(H_l)
```


```
fig, ax = plt.subplots(1, 1, figsize=(12,4))

# highest frequency (pi) corresponds to M * Fs/2
ax.plot(f_sampled, 20 * np.log10(abs(H_l_sampled)), 'b',  label=r"$ALIAS_{%d}$" % L)
ax.plot(f, 20 * np.log10(abs(H_l)), 'g',  label=r"$sinc$")
ax.set_ylabel('Amplitude [dB]', fontsize=12)
ax.set_xlabel('Frequency [Hz]', fontsize=12)
ax.grid()
ax.legend(loc=0)

n = np.arange(0,9,1)
ax.set_xticks(Fs * n)
labels = ["$%d F_S$" % y for y in n[1:]]
labels.insert(0,"0")
ax.set_xticklabels(labels, fontsize=12)
ax.set_xlim([0, L*Fs/2]);
ax.set_ylim([-100, 10]);
```

096ab30778f022e72905d7c191d8dec9b899a329
ipynb
Jupyter Notebook
audio/LinearInterpolation.ipynb
brunodigiorgi/ipn-notes
c8840a45989f25442c1d800ef8acdf8c630cdafc
CC-BY-3.0
2018-03-07T13:46:17.000Z
2018-03-07T13:46:17.000Z
audio/LinearInterpolation.ipynb
brunodigiorgi/ipn-notes
c8840a45989f25442c1d800ef8acdf8c630cdafc
CC-BY-3.0
audio/LinearInterpolation.ipynb
brunodigiorgi/ipn-notes
c8840a45989f25442c1d800ef8acdf8c630cdafc
CC-BY-3.0
Qwen/Qwen-72B
1. YES
2. YES