<jupyter_start><jupyter_text><NAME>,Professorship Signal Theory and Digital Signal Processing,Institute of Communications Engineering (INT),Faculty of Computer Science and Electrical Engineering (IEF),University of Rostock,Germany Tutorial Digital Signal Processing**DFT to DTFT Interpolation with the Periodic Sinc Function**,Winter Semester 2021/22 (Course 24505)- lecture: https://github.com/spatialaudio/digital-signal-processing-lecture- tutorial: https://github.com/spatialaudio/digital-signal-processing-exercisesFeel free to contact lecturer <EMAIL><jupyter_code>import numpy as np
from matplotlib import pyplot as plt
from numpy.fft import fft, ifft
#from scipy.fft import fft, ifft
from scipy.special import diric<jupyter_output><empty_output><jupyter_text>DFT to DTFT Interpolation with the Periodic Sinc FunctionThe DFT spectrum $X[\mu]$ can be interpolated towards the DTFT spectrum $X(\Omega)$ with\begin{equation}X(\Omega)=\sum_{\mu=0}^{N-1}X[\mu]\cdot\mathrm{e}^{-\mathrm{j}\frac{\left(\Omega-\frac{2\pi}{N}\mu\right)(N-1)}{2}}\cdot\text{psinc}_N\left(\Omega-\frac{2\pi}{N}\mu\right).\end{equation}This interpolation implies:- the DFT $X[\mu]$ stems from a signal $x[k]$ for which periodicity of $N$ is inherent, we assume the first period at $0\leq k\leq N-1$ - the DFT spectrum is discrete and $N$ periodic- the DTFT spectrum stems from the signal $x[k]$ for $0\leq k\leq N-1$, for all other $k$ is $x[k]=0$, i.e. considering only one period- the DTFT spectrum is continuous and $2\pi$ periodic, which means that the signal $x[k]$ can not be periodicThe interpolation kernel utilises the so-called **periodic sinc function**\begin{align}\text{psinc}_N(\Omega)=\begin{cases}\frac{1}{N}\cdot\frac{\sin\left(\frac{N}{2}\Omega\right)}{\sin\left(\frac{1}{2}\Omega\right)}&\text{for }\Omega\neq2\pi m\\(-1)^{m(N-1)}&\text{for }\Omega=2\pi m\end{cases},\,\,m\in\mathbb{Z},\end{align}which is also known as aliased sinc and Dirichlet function.Below, we give an example graph for $\text{psinc}_N(\Omega)$.Note that the orange dots indicate the DFT eigenfrequencies.<jupyter_code>def dft2dtft(X, W):
    N = len(X)
    Xi = np.zeros(len(W), dtype='complex')
    for i, Omega in enumerate(W):
        for mu_dft in range(N):
            Xi[i] += X[mu_dft] *\
                np.exp(-1j/2*(Omega - 2*np.pi/N*mu_dft)*(N-1)) *\
                diric(Omega - 2*np.pi/N*mu_dft, N)
    return Xi
K = 3
N_DFT = 2**K
Omega_DFT = np.arange(-2*N_DFT, 2*N_DFT) * 2*np.pi/N_DFT
N = 2**(K+4)
Omega = np.arange(-2*N, 2*N) * 2*np.pi/N
psinc_DFT = diric(Omega_DFT, N_DFT)
psinc = diric(Omega, N_DFT)

plt.plot(Omega, psinc, label='psinc function')
plt.plot(Omega_DFT, psinc_DFT, 'o', label='DFT bins')
plt.xlim(-4*np.pi, 4*np.pi)
plt.xticks(np.arange(-4, 5, 1)*np.pi)
plt.xlabel(r'$\Omega$')
plt.ylabel(r'psinc($\Omega$)')
plt.title('psinc for '+str(N_DFT)+'-point DFT')
plt.legend()
plt.grid(True)<jupyter_output><empty_output><jupyter_text>Task: DFT Analysis Using a Rectangular WindowA complex signal- $x[k]=\mathrm{exp}(\mathrm{j}(\Omega k + \frac{\pi}{\pi}))$ with- $\Omega=4\cdot\frac{2\pi}{N}$, - $N=8$is to be analysed with the DFT \begin{align}X[\mu]=\sum_{k=0}^{N-1}x[k]\cdot\mathrm{e}^{-\mathrm{j}\frac{2\pi}{N}k\mu}\end{align}for $0\leq k \leq N-1$.Furthermore, assume that $x[k]$ results from continuous-time signal $x(t)$ using sampling frequency of $f_s=10$ Hz.1. Plot the discrete-time signals that correspond to the DFT and the DTFT spectrum.2. Calculate the DFT spectrum $X[\mu]$ of $x[k]$ and visualise the real and imaginary part as well as the magnitude and the phase of $X[\mu]$ over $0\leq\mu\leq N-1$.3. Implement the above mentioned interpolation towards the DTFT and visualise the resulting magnitude spectra $|X[\mu]|$, $|X(\Omega)|$ over frequency axes $\mu$, $\Omega$ as well as the physical frequency $f$.  Task 1: Generate Signal and Plot<jupyter_code>N = 8

Om = 4 * 2*np.pi/N  # play with the factor 4

k = np.arange(N)
x = np.exp(1j*(Om*k+1))  # pi/pi=1

koffs = [-2, -1, 1, 2]
plt.subplot(2, 1, 1)
plt.plot(k, x.real, 'C0o-', label='real', ms=5)
plt.plot(k, x.imag, 'C1o-', label='imag', ms=5)
# the DFT spectrum corresponds to a periodic signal of period N
# we exemplarily indicate this
for koffsi in koffs:
    plt.plot(k+koffsi*N, x.real, 'C0o-', ms=3)
    plt.plot(k+koffsi*N, x.imag, 'C1o-', ms=3)
plt.xlim(koffs[0]*N, koffs[-1]*N+N-1)
plt.xlabel(r'$k$')
plt.ylabel(r'$x_\mathrm{DFT}[k]$')
plt.legend()
plt.grid(True)
plt.subplot(2, 1, 2)
plt.plot(k, x.real, 'C0o-', ms=5)
plt.plot(k, x.imag, 'C1o-', ms=5)
# the DTFT spectrum corresponds to singular event of length N for k=0...N-1
# we exemplarily indicate this
for koffsi in koffs:
    plt.plot(k+koffsi*N, 0*x.real, 'C0o-', ms=3)
    plt.plot(k+koffsi*N, 0*x.imag, 'C1o-', ms=3)
plt.xlim(koffs[0]*N, koffs[-1]*N+N-1)
plt.xlabel(r'$k$')
plt.ylabel(r'$x_\mathrm{DTFT}[k]$')
plt.grid(True)<jupyter_output><empty_output><jupyter_text>Task 2: DFT and Plot<jupyter_code>mu = np.arange(N)
w = np.ones(N)

# optional playground for windowing
# see dft_windowing_tutorial.pdf for the equations of these windows
# note that the scipy.signal.windows with same name behave different !
# thus standardised names does NOT imply standardised characteristics
if False:
    # Hann, two unused zeros
    w = (1 - np.cos(2*np.pi/N*(k+1/2))) / 2 
    # Hamming, put notch into first sidelobe, thus all zeros used
    w = 0.54 - 0.46 * np.cos(2*np.pi/N*(k+1/2))

X = fft(x*w)
plt.figure(figsize=(8, 4))

plt.subplot(2, 2, 1)
plt.plot(mu, np.real(X), 'o')
plt.xlim(0, N-1)
plt.xlabel(r'$\mu$')
plt.ylabel(r'$\Re(X[\mu])$')
plt.grid(True)

plt.subplot(2, 2, 2)
plt.plot(mu, np.imag(X), 'o')
plt.xlim(0, N-1)
plt.xlabel(r'$\mu$')
plt.ylabel(r'$\Im(X[\mu])$')
plt.grid(True)

plt.subplot(2, 2, 3)
plt.plot(mu, np.abs(X), 'o')
plt.xlim(0, N-1)
plt.xlabel(r'$\mu$')
plt.ylabel(r'$|X[\mu]|$')
plt.grid(True)

# note that only the phase at mu=4 is a meaningful result (=1)
# all other values are due to numerical noise / trash
# since Re and Im are approx zero at these mu
plt.subplot(2, 2, 4)
plt.plot(mu, np.angle(X), 'o')
plt.xlim(0, N-1)
plt.xlabel(r'$\mu$')
plt.ylabel(r'$\angle X[\mu]$')
plt.grid(True)<jupyter_output><empty_output><jupyter_text>Task 3: DFT -> DTFT Interpolation, Plot<jupyter_code>fs = 10
Ni = 2**8
W = np.arange(Ni) * 2*np.pi/Ni  # Omega for which to interpolate to DTFT
f = W / (2*np.pi) * fs  # physical frequency based on W using fs
mui = W / (2*np.pi) * N  # frequency vector normalized to integer frequencies
Xi = dft2dtft(X, W)  # DTFT interpolation from DFT

plt.figure(figsize=(10, 10))

plt.subplot(3, 2, 1)
plt.plot(mu, np.abs(X), 'C0o', label='DFT')
plt.plot(mui, np.abs(Xi), 'C1-', label='DTFT')
plt.xlabel(r'$\mu$')
plt.ylabel('|X|')
plt.legend()
plt.grid(True)

plt.subplot(3, 2, 2)
plt.plot(mu, 20*np.log10(np.abs(X)), 'C0o', label='DFT')
plt.plot(mui, 20*np.log10(np.abs(Xi)), 'C1-', label='DTFT')
plt.ylim(-100, 20)
plt.xlabel(r'$\mu$')
plt.ylabel('20lg|X| / dB')
plt.grid(True)

plt.subplot(3, 2, 3)
plt.plot(2*np.pi/N*mu, np.abs(X), 'C0o', label='DFT')
plt.plot(W, np.abs(Xi), 'C1-', label='DTFT')
plt.xlabel(r'$\Omega$')
plt.ylabel('|X|')
plt.grid(True)

plt.subplot(3, 2, 4)
plt.plot(2*np.pi/N*mu, 20*np.log10(np.abs(X)), 'C0o', label='DFT')
plt.plot(W, 20*np.log10(np.abs(Xi)), 'C1-', label='DTFT')
plt.ylim(-100, 20)
plt.xlabel(r'$\Omega$')
plt.ylabel('20lg|X| / dB')
plt.grid(True)

plt.subplot(3, 2, 5)
plt.plot(mu*fs/N, np.abs(X), 'C0o', label='DFT')
plt.plot(f, np.abs(Xi), 'C1-', label='DTFT')
plt.xlabel(r'f / Hz')
plt.ylabel('|X|')
plt.grid(True)

plt.subplot(3, 2, 6)
plt.plot(mu*fs/N, 20*np.log10(np.abs(X)), 'C0o', label='DFT')
plt.plot(f, 20*np.log10(np.abs(Xi)), 'C1-', label='DTFT')
plt.ylim(-100, 20)
plt.xlabel(r'f / Hz')
plt.ylabel('20lg|X| / dB')
plt.grid(True)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.19787533597849738, 'fraction_numerical': 0.030718034045821067, 'mean_word_length': 4.703649635036497, 'pattern_counts': {'":': 0, '<': 23, '<?xml version=': 0, '>': 24, 'https://': 2, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '5', 'dir': 'jupyter-structured-clean-dedup', 'id': '436192', 'license': 'CC-BY-4.0', 'n_tokens_mistral': 3509, 'n_tokens_neox': 3227, 'n_words': 772, 'path': 'dft/dft_to_dtft_interpolation.ipynb', 'repo_name': 'spatialaudio/digital-signal-processing-exercises', 'size': '7840'}