diff --git a/dft/dft_complex_exponential.ipynb b/dft/dft_complex_exponential.ipynb index a51d52f..3fe26b2 100644 --- a/dft/dft_complex_exponential.ipynb +++ b/dft/dft_complex_exponential.ipynb @@ -19,7 +19,7 @@ "- lecture: https://github.com/spatialaudio/digital-signal-processing-lecture\n", "- tutorial: https://github.com/spatialaudio/digital-signal-processing-exercises\n", "\n", - "Feel free to contact lecturer jacob.thoenes@uni-rostock.de" + "Feel free to contact the lecturer jacob.thoenes@uni-rostock.de" ] }, { @@ -28,8 +28,8 @@ "metadata": {}, "outputs": [], "source": [ - "import numpy as np\n", "import matplotlib.pyplot as plt\n", + "import numpy as np\n", "from numpy.fft import fft, ifft\n", "\n", "# from scipy.fft import fft, ifft\n", @@ -100,9 +100,7 @@ "# print(np.allclose(x, x_tmp))\n", "\n", "# calc DFT -> DTFT interpolation\n", - "N_dtft = (\n", - " 2**10\n", - ") # number of frequencies along unit circle at which DTFT values are calc\n", + "N_dtft = 2**10 # number of frequencies along unit circle at which DTFT values are calc\n", "Om_dtft = np.arange(N_dtft) * 2 * np.pi / N_dtft # set up frequency vector\n", "X_dtft = dft2dtft(X, Om_dtft) # perform interp" ] @@ -193,7 +191,7 @@ "- the signal $x[k]$ is zero for $k<0$ and $k>N-1$, thus **non-periodic**\n", "- the signal is discrete\n", "- the DTFT spectrum is periodic in $2\\pi$\n", - "- the DTFT spectrum is continous" + "- the DTFT spectrum is continuous" ] }, { @@ -304,7 +302,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.10.12" + "version": "3.12.2" } }, "nbformat": 4, diff --git a/dft/dft_intro.ipynb b/dft/dft_intro.ipynb index 911d760..74df857 100644 --- a/dft/dft_intro.ipynb +++ b/dft/dft_intro.ipynb @@ -19,7 +19,7 @@ "- lecture: https://github.com/spatialaudio/digital-signal-processing-lecture\n", "- tutorial: https://github.com/spatialaudio/digital-signal-processing-exercises\n", "\n", - "Feel free to contact lecturer jacob.thoenes@uni-rostock.de" + "Feel free to contact the lecturer jacob.thoenes@uni-rostock.de" ] }, { @@ -28,10 +28,10 @@ "metadata": {}, "outputs": [], "source": [ - "import numpy as np\n", "import matplotlib.pyplot as plt\n", - "from numpy.linalg import inv\n", + "import numpy as np\n", "from numpy.fft import fft, ifft\n", + "from numpy.linalg import inv\n", "\n", "# from scipy.fft import fft, ifft" ] @@ -52,7 +52,7 @@ "Let us first define a **complex-valued signal** $x[k]$ of a certain block length $N$ ranging from $0\\leq k\\leq N-1$.\n", "\n", "The variable `tmpmu` defines the frequency of the signal. We will see later how this is connected to the DFT.\n", - "For now on, leave it with `tmpmu=1`. This results in exactly one period of cosine and sine building the complex signal. If `tmpmu=2` we get exactly two periods of cos/sin. We'll get get an idea of `tmpmu`..." + "From now on, leave it with `tmpmu=1`. This results in exactly one period of cosine and sine building the complex signal. If `tmpmu=2`, we get exactly two periods of cos/sin. We'll get an idea of `tmpmu`..." ] }, { @@ -87,7 +87,7 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "We will now perform an DFT of $x[k]$ since we are interested in the frequency spectrum of it.\n", + "We will now perform a DFT of $x[k]$ since we are interested in the frequency spectrum of it.\n", "\n", "## DFT Definition\n", "\n", @@ -97,7 +97,7 @@ "\\text{IDFT}: x[k]=\\frac{1}{N}&\\sum_{\\mu=0}^{N-1}X[\\mu]\\cdot\\mathrm{e}^{+\\mathrm{j}\\frac{2\\pi}{N}k\\mu}\n", "\\end{align}\n", "\n", - "Note the sign reversal in the exp()-function and the $1/N$ normalization in the IDFT. This convention is used by the majority of DSP text books and also in Python's `numpy.fft.fft()`, `numpy.fft.ifft()` and Matlab's `fft()`, `ifft()` routines." + "Note the sign reversal in the exp()-function and the $1/N$ normalization in the IDFT. This convention is used by the majority of DSP textbooks and also in Python's `numpy.fft.fft()`, `numpy.fft.ifft()` and Matlab's `fft()`, `ifft()` routines." ] }, { @@ -106,11 +106,11 @@ "source": [ "## DFT and IDFT with For-Loops\n", "\n", - "We are now going to implement the DFT and IDFT with for-loop handling. While this might be helpful to validate algorithms in its initial development phase, this should be avoided for practical used code in the field: for-loops are typically slow and very often more complicated to read than appropriate set up matrices and vectors. Especially for very large $N$ the computation time is very long.\n", + "We are now going to implement the DFT and IDFT with a for-loop handling. While this might be helpful to validate algorithms in their initial development phase, this should be avoided for practical use code in the field: for-loops are typically slow and very often more complicated to read than an appropriate setup of matrices and vectors. Especially for very large $N$, the computation time is very long.\n", "\n", "Anyway, the for-loop concept is: the DFT can be implemented with an outer for loop iterating over $\\mu$ and an inner for loop summing over all $k$ for a specific $\\mu$.\n", "\n", - "We use variable with _ subscript here, in order to save nice variable names for the matrix based calculation." + "We use a variable with _ subscript here, in order to save nice variable names for the matrix based calculation." ] }, { @@ -160,9 +160,9 @@ "source": [ "## DFT and IDFT with Matrix Multiplication\n", "\n", - "Now we do a little better: We should think of the DFT/IDFT in terms of a matrix operation setting up a set of linear equations.\n", + "Now we do a little better: We should think of the DFT/IDFT in terms of a matrix operation, setting up a set of linear equations.\n", "\n", - "For that we define a column vector containing the samples of the discrete-time signal $x[k]$\n", + "For that, we define a column vector containing the samples of the discrete-time signal $x[k]$\n", "\\begin{equation}\n", "\\mathbf{x}_k = (x[k=0], x[k=1], x[k=2], \\dots , x[k=N-1])^\\mathrm{T}\n", "\\end{equation}\n", @@ -214,7 +214,7 @@ "\\end{equation}\n", "containing all possible products $k\\,\\mu$ in a suitable arrangement.\n", "\n", - "For the simple case $N=4$ these matrices are\n", + "For the simple case $N=4$, these matrices are\n", "\\begin{align}\n", "\\mathbf{K} = \\begin{bmatrix}\n", "0 & 0 & 0 & 0\\\\\n", @@ -296,7 +296,7 @@ "source": [ "## Fourier Matrix Properties\n", "\n", - "The DFT and IDFT basically solve two sets of linear equations, that are linked as forward and inverse problem.\n", + "The DFT and IDFT basically solve two sets of linear equations, which are linked as a forward and inverse problem.\n", "\n", "This is revealed with the important property of the Fourier matrix\n", "\n", @@ -323,16 +323,16 @@ "(\\frac{\\mathbf{W}}{\\sqrt{N}})^\\mathrm{H} \\, (\\frac{\\mathbf{W}}{\\sqrt{N}}) = \\mathbf{I} =\n", "(\\frac{\\mathbf{W}}{\\sqrt{N}})^{-1} \\, (\\frac{\\mathbf{W}}{\\sqrt{N}})\n", "\\end{equation}\n", - "holds, i.e. the complex-conjugate, transpose is equal to the inverse\n", + "holds, i.e., the complex-conjugate transpose is equal to the inverse\n", "$(\\frac{\\mathbf{W}}{\\sqrt{N}})^\\mathrm{H} = (\\frac{\\mathbf{W}}{\\sqrt{N}})^{-1}$\n", - "and due to the matrix symmetry also\n", + "and due to the matrix symmetry, also\n", "$(\\frac{\\mathbf{W}}{\\sqrt{N}})^* =\n", "(\\frac{\\mathbf{W}}{\\sqrt{N}})^{-1}$\n", "is valid.\n", "\n", - "This tells that the matrix $\\frac{\\mathbf{W}}{\\sqrt{N}}$ is **orthonormal**, i.e. the matrix spans a orthonormal vector basis (the best what we can get in linear algebra world to work with) of $N$ normalized DFT eigensignals.\n", + "This tells that the matrix $\\frac{\\mathbf{W}}{\\sqrt{N}}$ is **orthonormal**, i.e., the matrix spans an orthonormal vector basis (the best we can get in the linear algebra world to work with) of $N$ normalized DFT eigensignals.\n", "\n", - "So, DFT and IDFT is transforming vectors into other vectors using the vector basis of the Fourier matrix.\n" + "So, DFT and IDFT are transforming vectors into other vectors using the vector basis of the Fourier matrix.\n" ] }, { @@ -348,7 +348,7 @@ "since we have intentionally set up the matrix this way.\n", "\n", "The plot below shows the eigensignal for $\\mu=1$, which fits again one signal period in the block length $N$.\n", - "For $\\mu=2$ we obtain two periods in one block.\n", + "For $\\mu=2$, we obtain two periods in one block.\n", "\n", "The eigensignals for $0\\leq \\mu \\leq N-1$ therefore exhibit a certain digital frequency, the so called DFT eigenfrequencies.\n", "\n", @@ -383,11 +383,11 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "The nice thing about the chosen eigenfrequencies, is that the eigensignals are **orthogonal**.\n", + "The nice thing about the chosen eigenfrequencies is that the eigensignals are **orthogonal**.\n", "\n", "This choice of the vector basis is on purpose and one of the most important ones in linear algebra and signal processing.\n", "\n", - "We might for example check orthogonality with the **complex** inner product of some matrix columns." + "We might, for example, check orthogonality with the **complex** inner product of some matrix columns." ] }, { @@ -419,13 +419,13 @@ "\n", "Let us synthesize a discrete-time signal by using the IDFT in matrix notation for $N=8$.\n", "\n", - "The signal should contain a DC value, the first and second eigenfrequency with different amplitudes, such as\n", + "The signal should contain a DC value, the first and second eigenfrequencies with different amplitudes, such as\n", "\n", "\\begin{equation}\n", "\\mathbf{x}_\\mu = [8, 2, 4, 0, 0, 0, 0, 0]^\\text{T}\n", "\\end{equation}\n", "\n", - "using large `X_test` in code." + "using a large `X_test` in code." ] }, { @@ -480,7 +480,7 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "We don't need summing the other columns, since their DFT coefficients in `X_test` are zero.\n", + "We don't need to sum up the other columns, since their DFT coefficients in `X_test` are zero.\n", "\n", "Finally, normalizing yields the IDFT." ] @@ -502,7 +502,7 @@ "source": [ "## Initial Example: DFT Spectrum Analysis for N=8\n", "\n", - "Now, let us calculate the DFT of the signal `x_test`. As result, we'd expect the DFT vector\n", + "Now, let us calculate the DFT of the signal `x_test`. As a result, we'd expect the DFT vector\n", "\n", "\\begin{equation}\n", "\\mathbf{x}_\\mu = [8, 2, 4, 0, 0, 0, 0, 0]^\\text{T}\n", @@ -552,17 +552,17 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "The analysis stage for the discrete-time signal domain, i.e. the DFT\n", + "The analysis stage for the discrete-time signal domain, i.e., the DFT\n", "can be reinvented by some intuition:\n", "How 'much' of the reference signal $\\mathbf{w}_{\\text{column i}}$\n", "(any column in $\\mathbf{W}$)\n", "is contained in the discrete-time signal $\\mathbf{x}_k$ that is to be analysed.\n", "\n", - "In signal processing / statistic terms we look for the amount of correlation\n", + "In signal processing/statistics terms, we look for the amount of correlation\n", "of the signals\n", "$\\mathbf{w}_{\\text{column i}}$ and $\\mathbf{x}_k$.\n", "\n", - "In linear algebra terms we are interested in the projection of $\\mathbf{x}_k$ onto\n", + "In linear algebra terms, we are interested in the projection of $\\mathbf{x}_k$ onto\n", "$\\mathbf{w}_{\\text{column i}}$, because the resulting length of this vector\n", "reveals the amount of correlation, which is precisely one DFT coefficient $X[\\cdot]$.\n", "\n", @@ -597,7 +597,7 @@ "X[\\mu=N-1] =& \\mathbf{w}_{\\text{column N}}^\\text{H} \\cdot \\mathbf{x}_k.\n", "\\end{align}\n", "\n", - "Naturally, all operations can be merged to one single\n", + "Naturally, all operations can be merged into one single\n", "matrix multiplication using the conjugate transpose of $\\mathbf{W}$.\n", "\n", "\\begin{equation}\n", @@ -635,7 +635,7 @@ "source": [ "Next, let us plot the magnitude of the spectrum over $\\mu$.\n", "\n", - "- We should play around with the variable `tmpmu` when defining the input signal at the very beginning of the notebook. For example we can check what happens for `tmpmu = 1`, `tmpmu = 2` and run the whole notebook to visualize the actual magnitude spectra.\n", + "- We should play around with the variable `tmpmu` when defining the input signal at the very beginning of the notebook. For example, we can check what happens for `tmpmu = 1`, `tmpmu = 2`, and run the whole notebook to visualize the actual magnitude spectra.\n", "\n", "We should recognize the link of the 'energy' at $\\mu$ in the magnitude spectrum with the chosen `tmpmu`.\n", "\n", @@ -685,7 +685,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.10.12" + "version": "3.12.2" } }, "nbformat": 4, diff --git a/dft/dft_to_dtft_interpolation.ipynb b/dft/dft_to_dtft_interpolation.ipynb index fb8c8e7..a3bb940 100644 --- a/dft/dft_to_dtft_interpolation.ipynb +++ b/dft/dft_to_dtft_interpolation.ipynb @@ -19,7 +19,7 @@ "- lecture: https://github.com/spatialaudio/digital-signal-processing-lecture\n", "- tutorial: https://github.com/spatialaudio/digital-signal-processing-exercises\n", "\n", - "Feel free to contact lecturer jacob.thoenes@uni-rostock.de\n" + "Feel free to contact the lecturer jacob.thoenes@uni-rostock.de" ] }, { @@ -61,7 +61,7 @@ "(-1)^{m(N-1)}&\\text{for }\\Omega=2\\pi m\\end{cases},\\,\\,m\\in\\mathbb{Z},\n", "\\end{align}\n", "\n", - "which is also known as aliased sinc and Dirichlet function.\n", + "which is also known as the aliased sinc and the Dirichlet function.\n", "\n", "Below, we give an example graph for $\\text{psinc}_N(\\Omega)$.\n", "Note that the orange dots indicate the DFT eigenfrequencies." @@ -132,13 +132,13 @@ "\\end{align}\n", "for $0\\leq k \\leq N-1$.\n", "\n", - "Furthermore, assume that $x[k]$ results from continuous-time signal $x(t)$ using sampling frequency of $f_s=10$ Hz.\n", + "Furthermore, assume that $x[k]$ results from a continuous-time signal $x(t)$ using a sampling frequency of $f_s=10$ Hz.\n", "\n", "1. Plot the discrete-time signals that correspond to the DFT and the DTFT spectrum.\n", "\n", - "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$.\n", + "2. Calculate the DFT spectrum $X[\\mu]$ of $x[k]$ and visualise the real and imaginary parts as well as the magnitude and the phase of $X[\\mu]$ over $0\\leq\\mu\\leq N-1$.\n", "\n", - "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$." + "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$." ] }, { @@ -331,19 +331,19 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "We observe that 7 of 8 DFT coeefficients are precisely zero, and only $|X[\\mu=4]|=8$. This is intentional and is expected for the chosen signal that oscillates with exactly a DFT eigenfrequency, here half of the sampling frequency.\n", + "We observe that 7 of 8 DFT coefficients are precisely zero, and only $|X[\\mu=4]|=8$. This is intentional and is expected for the chosen signal that oscillates with exactly a DFT eigenfrequency, here, half of the sampling frequency.\n", "\n", "The code above allows for two more playgrounds:\n", "\n", "- We might change $\\Omega=1\\cdot\\frac{2\\pi}{N}$, $\\Omega=2\\cdot\\frac{2\\pi}{N}$, $\\Omega=3\\cdot\\frac{2\\pi}{N}$\n", - "and observe how the spectrum moves around the frequency axis, however the shape is preserved.\n", + "and observe how the spectrum moves around the frequency axis, however, the shape is preserved.\n", "The main lobe precisely corresponds to the chosen $\\mu$.\n", "\n", "- We might change $\\Omega=1.5\\cdot\\frac{2\\pi}{N}$, $\\Omega=2.5\\cdot\\frac{2\\pi}{N}$, $\\Omega=3.5\\cdot\\frac{2\\pi}{N}$\n", - "and observe that - since these frequencies are **no** DFT eigenfrequencies - now all DFT coefficients exhibit energy. This suggests frequencies in the signal which are actually not there, but originate from cutting the signal to the chosen block size $N$. The effect is known as **leakage effect**.\n", + "and observe that - since these frequencies are **no** DFT eigenfrequencies - now all DFT coefficients exhibit energy. This suggests frequencies in the signal that are actually not there, but originate from cutting the signal to the chosen block size $N$. The effect is known as the **leakage effect**.\n", "\n", - "- In task 2 we might optionally apply a window to the signal. The plots for task 3 then reveal how the leakage effect is reduced but the main lobe gets broader.\n", - "Since we have chosen a mono-frequent, complex oscillation as input signal to the DFT / DTFT, the actual window spectrum can be directly seen in the plots." + "- In task 2, we might optionally apply a window to the signal. The plots for task 3 then reveal how the leakage effect is reduced, but the main lobe gets broader.\n", + "Since we have chosen a mono-frequent, complex oscillation as an input signal to the DFT / DTFT, the actual window spectrum can be directly seen in the plots." ] }, { @@ -372,7 +372,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.10.12" + "version": "3.12.2" } }, "nbformat": 4, diff --git a/dft/dft_windowing_complex_exponential.ipynb b/dft/dft_windowing_complex_exponential.ipynb index 7126644..d18fdc0 100644 --- a/dft/dft_windowing_complex_exponential.ipynb +++ b/dft/dft_windowing_complex_exponential.ipynb @@ -19,7 +19,7 @@ "- lecture: https://github.com/spatialaudio/digital-signal-processing-lecture\n", "- tutorial: https://github.com/spatialaudio/digital-signal-processing-exercises\n", "\n", - "Feel free to contact lecturer jacob.thoenes@uni-rostock.de" + "Feel free to contact the lecturer jacob.thoenes@uni-rostock.de" ] }, { @@ -34,7 +34,7 @@ "x[k] \\cdot w[k] \\circ-\\bullet \\frac{1}{2\\pi} X(\\Omega) \\circledast_{2\\pi} W(\\Omega)\n", "\\end{equation}\n", "\n", - "Typically for the window function we choose non-zero values in $w[k]$ for $0 \\leq k \\leq N-1$ and $w[k]=0$ elsewhere.\n", + "Typically, for the window function, we choose non-zero values in $w[k]$ for $0 \\leq k \\leq N-1$ and $w[k]=0$ elsewhere.\n", "\n", "Thus, we deal with a length $N$ signal that is non-zero at $0 \\leq k \\leq N-1$.\n", "We can apply an $N$ DFT onto this signal $x_w[k] = x[0 \\leq k \\leq N-1] \\cdot w[0 \\leq k \\leq N-1]$ to obtain its DFT spectrum $X_w[\\mu]$ with the baseband $0\\leq\\mu\\leq N-1$.\n", @@ -48,7 +48,7 @@ "source": [ "## Special Case: Complex Exponential\n", "\n", - "For a certain frequency $\\Omega_1$ we can invent a simple complex exponential signal\n", + "For a certain frequency $\\Omega_1$, we can invent a simple complex exponential signal\n", "\n", "\\begin{equation}\n", "x_1[k] = \\mathrm{e}^{+\\mathrm{j} \\Omega_1 k} \\circ-\\bullet \\bot\\bot\\bot (\\frac{\\Omega-\\Omega_1}{2\\pi})\n", @@ -56,10 +56,10 @@ "\n", "For certain window / DFT length $N$, the $\\Omega_1 = \\frac{2\\pi}{N} \\cdot \\mu$ for $0\\leq\\mu\\leq N-1$ correspond to the DFT eigenfrequencies in the spectrum baseband. On the other hand, if $\\mu$ is rather non-integer, the frequency is between two DFT eigenfrequencies.\n", "\n", - "Since the complex exponential signal exhibits a Dirac comb in DTFT domain, the above circular convolution for the windowing process is very convenient to discuss:\n", + "Since the complex exponential signal exhibits a Dirac comb in the DTFT domain, the above circular convolution for the windowing process is very convenient to discuss:\n", "\n", "The Dirac is the neutral element of the convolution.\n", - "Thus, to analyze the spectrum of $x_1[k]$ (actually its windowed part) we just need to interpret the DTFT spectrum $W(\\Omega-\\Omega_1)$, i.e. the window spectrum that is circular shifted along frequency axis." + "Thus, to analyze the spectrum of $x_1[k]$ (actually its windowed part), we just need to interpret the DTFT spectrum $W(\\Omega-\\Omega_1)$, i.e., the window spectrum that is circularly shifted along the frequency axis." ] }, { @@ -102,9 +102,9 @@ "metadata": {}, "outputs": [], "source": [ - "import numpy as np\n", "import matplotlib.pyplot as plt\n", - "from numpy.fft import fft, ifft, fftshift\n", + "import numpy as np\n", + "from numpy.fft import fft, fftshift, ifft\n", "\n", "# from scipy.fft import fft, ifft, fftshift\n", "from scipy.signal.windows import kaiser\n", @@ -206,21 +206,21 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "Note that the scaling of DTFT is aligned to the DFT data for better visual cues.\n", - "We should clarify this even more in future.\n", + "Note that the scaling of the DTFT is aligned to the DFT data for better visual cues.\n", + "We should clarify this even more in the future.\n", "\n", "In the upcoming examples, we always have two plots\n", "- the upper graph over DFT frequency index $\\mu$ shows the two magnitude DTFT spectra of $W(\\Omega-\\Omega_1)$, $W(\\Omega-\\Omega_2)$ and there corresponding DFT spectra (dots)\n", - "- the lower plots over DTFT frequency $\\Omega$ shows the DTFT / DFT spectrum of the signal / spectrum superposition\n", + "- the lower plots over DTFT frequency $\\Omega$ show the DTFT / DFT spectrum of the signal/spectrum superposition\n", "\n", - "Thus, in the upper plot we see how two Dirac impulses are smeared by the window spectrum. We can study the characteristics of\n", + "Thus, in the upper plot, we see how two Dirac impulses are smeared by the window spectrum. We can study the characteristics of\n", "- the so called mainlobe (the shape around the spectrum's maximum) and\n", "- the so called sidelobes (the shapes that are followed by the mainlobe)\n", - "- number of zeros (for optimum windows we have $N-1$ zeros in the DTFT spectrum). Note that $\\Omega=0=2\\pi$ is the same zero, due to periodicity.\n", + "- number of zeros (for optimum windows, we have $N-1$ zeros in the DTFT spectrum). Note that $\\Omega=0=2\\pi$ is the same zero, due to periodicity.\n", "- number of non-zeros coefficients in DFT\n", "\n", - "In the lower plot we see how the superposition of the two signals affects the resulting spectrum.\n", - "- Do we obtain a separate mainlobes for the two frequencies?\n", + "In the lower plot, we see how the superposition of the two signals affects the resulting spectrum.\n", + "- Do we obtain separate main lobes for the two frequencies?\n", "- Or do we obtain a rather broad single mainlobe, not able to tell that the signal is built from two frequencies?\n", "- Do we have zeros in the spectrum?\n", "- How about the sidelobe magnitude and the decay of the sidelobes?" @@ -385,7 +385,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.10.12" + "version": "3.12.2" } }, "nbformat": 4, diff --git a/dft/dft_windowing_tutorial/dft_windowing_tutorial_exercises.ipynb b/dft/dft_windowing_tutorial/dft_windowing_tutorial_exercises.ipynb index 277ff58..6f863d8 100644 --- a/dft/dft_windowing_tutorial/dft_windowing_tutorial_exercises.ipynb +++ b/dft/dft_windowing_tutorial/dft_windowing_tutorial_exercises.ipynb @@ -19,7 +19,7 @@ "- lecture: https://github.com/spatialaudio/digital-signal-processing-lecture\n", "- tutorial: https://github.com/spatialaudio/digital-signal-processing-exercises\n", "\n", - "Feel free to contact lecturer jacob.thoenes@uni-rostock.de" + "Feel free to contact the lecturer jacob.thoenes@uni-rostock.de" ] }, { @@ -246,7 +246,7 @@ "and specify the missing value $X[\\mu=3]$ so that the IDFT results in\n", "$x[k]\\in\\mathbb{R}$. Certain symmetries in the spectrum hold.\n", "\n", - "Then calculate the IDFT analytically as\n", + "Then, calculate the IDFT analytically as\n", "\\begin{equation}\n", "x[k]=\\frac{1}{N}\\sum_{\\mu=0}^{N-1}X[\\mu]\\cdot\\mathrm{e}^{\\mathrm{j}\\frac{2\\pi}{N}k\\mu}\n", "\\end{equation}\n", @@ -297,7 +297,7 @@ "source": [ "**c)**\n", "\n", - "Plot the real and imaginary part as well as the magnitude and the phase of $X[\\mu]$ over $\\mu$." + "Plot the real and imaginary parts as well as the magnitude and the phase of $X[\\mu]$ over $\\mu$." ] }, { @@ -343,7 +343,7 @@ "\\hspace{5mm}\\text{for}\\,\\,0\\leq k\\leq3\n", "\\end{equation}\n", "\n", - "by help of a linear combination using the Fourier matrix." + "by the help of a linear combination using the Fourier matrix." ] }, { @@ -434,7 +434,7 @@ "(-1)^{m(N-1)}&\\text{for }\\Omega=2\\pi m\\end{cases},\\,\\,m\\in\\mathbb{Z},\n", "\\end{align}\n", "\n", - "which is also known as aliased sinc and Dirichlet function.\n", + "which is also known as the aliased sinc and the Dirichlet function.\n", "This interpolation implies:\n", "- 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$ \n", "- the DFT spectrum is discrete and $N$ periodic\n", @@ -473,7 +473,7 @@ "\n", "Implement the DTFT interpolation equation above and visualise the $|X[\\mu]|$ and $|X(\\Omega)|$ over the frequency axes $\\mu$, $\\Omega$ and $f$.\n", "\n", - "How and why the results are different?" + "How and why are the results different?" ] }, { @@ -482,7 +482,7 @@ "source": [ "**Solution**\n", "\n", - "It is meaningful to code two functions for repeated calculations, the first is the DFT to DTFT interpolation, the second function creates the signal and calculates the DFT, DTFT and plots the requested graphs." + "It is meaningful to code two functions for repeated calculations, the first is the DFT to DTFT interpolation, the second function creates the signal and calculates the DFT, DTFT, and plots the requested graphs." ] }, { @@ -646,10 +646,10 @@ "\n", "Obviously, the LTI system is non-recursive.\n", "\n", - "The magnitude and phase spectrum is to be evaluated from\n", + "The magnitude and phase spectrum are to be evaluated from\n", "- DFT\n", "- DTFT\n", - "- DFT of zerpadded impulse response\n", + "- DFT of zero-padded impulse response\n", "\n", "Discuss the different approaches." ] @@ -734,7 +734,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.10.12" + "version": "3.12.2" } }, "nbformat": 4, diff --git a/dft/dtft_to_dft_sampling.ipynb b/dft/dtft_to_dft_sampling.ipynb index f59cdab..19b55ff 100644 --- a/dft/dtft_to_dft_sampling.ipynb +++ b/dft/dtft_to_dft_sampling.ipynb @@ -19,7 +19,7 @@ "- lecture: https://github.com/spatialaudio/digital-signal-processing-lecture\n", "- tutorial: https://github.com/spatialaudio/digital-signal-processing-exercises\n", "\n", - "Feel free to contact lecturer jacob.thoenes@uni-rostock.de\n" + "Feel free to contact the lecturer jacob.thoenes@uni-rostock.de" ] }, { @@ -114,9 +114,9 @@ "source": [ "By plotting the spectrum of the DTFT and the DFT, we see that the DFT coefficients exactly correspond to DTFT values at specific frequencies. This is due to the inherent sampling process, where we **sample the DTFT spectrum** to **precisely obtain the DFT coefficients**. Detailed derivation can be found in the [DFT tutorial in Ch. 1.7](https://github.com/spatialaudio/digital-signal-processing-exercises/blob/master/dft/dft_windowing_tutorial/dft_windowing_tutorial.pdf) \n", "\n", - "**Sampling in one domain** (here the spectrum) corresponds to **'make it periodic' in the other domain** (here the discrete-time signal). The single window function that corresponds to the DTFT spectrum (blue) is thus repeated over and over again (orange), which is required for the DFT's periodicity characteristic. Its periodicity is determined by the DFT size and not necessarily equal to `N_DTFT`.\n", + "**Sampling in one domain** (here the spectrum) corresponds to **'make it periodic' in the other domain** (here the discrete-time signal). The single window function that corresponds to the DTFT spectrum (blue) is thus repeated over and over again (orange), which is required for the DFT's periodicity characteristic. Its periodicity is determined by the DFT size and is not necessarily equal to `N_DTFT`.\n", "\n", - "Very **large zeropadding** can be used to obtain many spectral samples and by that a very good **approximation of the DTFT** spectrum **without** using **explicit interpolation** between the DFT coefficients. This is often used in practice, very conveniently by calling `fft(x, N)` where `N` is larger than length of `x`. Note however, that by doing so, we do not increase the spectral resolution and we also do not enhance the capability for spectral discrimination. We gain **no new information by oversampling**.\n", + "Very **large zeropadding** can be used to obtain many spectral samples and by that a very good **approximation of the DTFT** spectrum **without** using **explicit interpolation** between the DFT coefficients. This is often used in practice, very conveniently by calling `fft(x, N)` where `N` is larger than the length of `x`. Note, however, that by doing so, we do not increase the spectral resolution, and we also do not enhance the capability for spectral discrimination. We gain **no new information by oversampling**.\n", "\n", "Note that `N_DTFT=N_DFT` (no zeropadding) is the case of **critical sampling**." ] @@ -219,7 +219,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.10.12" + "version": "3.12.2" } }, "nbformat": 4, diff --git a/dft/parametric_windows.ipynb b/dft/parametric_windows.ipynb index 3a2688b..1c36318 100644 --- a/dft/parametric_windows.ipynb +++ b/dft/parametric_windows.ipynb @@ -19,7 +19,7 @@ "- lecture: https://github.com/spatialaudio/digital-signal-processing-lecture\n", "- tutorial: https://github.com/spatialaudio/digital-signal-processing-exercises\n", "\n", - "Feel free to contact lecturer jacob.thoenes@uni-rostock.de" + "Feel free to contact the lecturer jacob.thoenes@uni-rostock.de" ] }, { @@ -29,20 +29,20 @@ "# Parametric Windows: Dolph-Chebyshev, Slepian, Kaiser-Bessel Window\n", "\n", "There are numerous window types, all developed for special requirements.\n", - "In the initial [DFT & Windowing](dft_windowing.ipynb) exercise and the detailed [DFT tutorial](dft_windowing_tutorial/dft_windowing_tutorial.pdf) we have learned about two very simple computable and very often used **non-parametric** windows, the Hann and the Hamming window.\n", - "Non-parametric means that by desired window length $M$, the window and its DTFT spectrum are fully determined.\n", - "In other words, if another spectral characteristics - very often we need higher spectral resolution - is asked for, then the only variable to change is $M$.\n", + "In the initial [DFT & Windowing](dft_windowing.ipynb) exercise and the detailed [DFT tutorial](dft_windowing_tutorial/dft_windowing_tutorial.pdf), we have learned about two very simple, computable, and very often used **non-parametric** windows, the Hann and the Hamming window.\n", + "Non-parametric means that by the desired window length $M$, the window and its DTFT spectrum are fully determined.\n", + "In other words, if another spectral characteristic - very often we need higher spectral resolution - is asked for, then the only variable to change is $M$.\n", "\n", "The **Hann** window is **not optimal**, since it does not use two of its potential zeros to shape the sidelobes of its DTFT spectrum.\n", "The Hamming window introduces two additional zeros to the Hann window to reduce the level of the first sidelobe.\n", - "Note that the Hann window is still often used nowadays, not due to its non-optimum spectrum, but rather due to its simple calculation of the window signal $w[k]$ requiring only a cosine and weight of $1/2$.\n", + "Note that the Hann window is still often used nowadays, not due to its non-optimum spectrum, but rather due to its simple calculation of the window signal $w[k]$ requiring only a cosine and a weight of $1/2$.\n", "\n", - "So called **parametric** windows have additional parameters, that allow to meet certain constraints for a given overall design criterion.\n", - "Two of the most prominent - in fact with these we probably can manage the majority of windowing applications - are the **Dolph-Chebyshev** and the **Kaiser-Bessel** window.\n", + "So called **parametric** windows have additional parameters that allow for meeting certain constraints for a given overall design criterion.\n", + "Two of the most prominent - in fact, with these we probably can manage the majority of windowing applications - are the **Dolph-Chebyshev** and the **Kaiser-Bessel** window.\n", "These are optimum window designs.\n", "The Kaiser-Bessel window itself can be considered as an approximation of the so called **discrete prolate spheroidal sequences** (DPSS, aka Slepian) window.\n", "\n", - "We will discuss these windows below in terms of their design criteria and the resulting additional parameter that can be set up to meet a desired constraint.\n", + "We will discuss these windows below in terms of their design criteria and the resulting additional parameters that can be set up to meet a desired constraint.\n", "TBD..." ] }, @@ -141,10 +141,10 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "We might figure out the design criterion of the Dolph-Chebyshev by ourselves, when inspecting the DTFT spectrum of it.\n", - "Hint: How is the additional parameter linked to the DTFT spectrum.\n", + "We might figure out the design criterion of the Dolph-Chebyshev by ourselves when inspecting the DTFT spectrum of it.\n", + "Hint: How is the additional parameter linked to the DTFT spectrum?\n", "Vary `M` and `SidelobeAttenuation` and check for changes.\n", - "For small `M` (to make analysis convenient) check how the zeros are placed in the spectrum to meet the design criterion." + "For small `M` (to make analysis convenient), check how the zeros are placed in the spectrum to meet the design criterion." ] }, { @@ -154,9 +154,9 @@ "## DPSS Window aka Slepian Window\n", "\n", "The design criterion of the **Slepian** (also named **discrete prolate spheroidal sequences** or **digital prolate spheroidal sequences** (DPSS)) window is **maximum energy** concentration in the **main lobe** for a given mainlobe **bandwidth**.\n", - "Actually, this is what we typically ask for signal analysis, if no other specific constraints about the sidelobes positions and there levels are requested.\n", + "Actually, this is what we typically ask for signal analysis, if no other specific constraints about the sidelobe positions and their levels are requested.\n", "\n", - "Recall, that the DTFT spectrum for the 'ideal world' window is the Dirac impulse (steming from the practically not feasible infinite rectangular window), so mainlobe energy concentration seems to be a very good approach to get close to it.\n", + "Recall that the DTFT spectrum for the 'ideal world' window is the Dirac impulse (stemming from the practically not feasible infinite rectangular window), so mainlobe energy concentration seems to be a very good approach to get close to it.\n", "\n", "See\n", "\n", @@ -168,7 +168,7 @@ "- Julius O. Smith (2011): Spectral Audio Signal Processing, online lecture of CCRMA, Stanford University, https://ccrma.stanford.edu/~jos/sasp/Slepian_DPSS_Window.html\n", "\n", "\n", - "for treatments how to derive the Slepian window.\n", + "for treatments, how to derive the Slepian window.\n", "\n", "Challenging question:\n", "What window type results if we ask the Slepian window to produce mainlobe bandwidth $\\rightarrow 0$?\n", @@ -192,7 +192,7 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "For the chosen example, this Slepian window has (i) the same window length and (ii) the first sidelobe is at about -50 dB like the Dolph-Chebyshev window above." + "For the chosen example, this Slepian window has (i) the same window length and (ii) the first sidelobe is at about -50 dB, like the Dolph-Chebyshev window above." ] }, { @@ -251,7 +251,7 @@ "## Kaiser-Bessel Window\n", "\n", "The Kaiser-Bessel window is an **approximation** of the Slepian window **for large window lengths** $M$, note however that they will be **never identical**.\n", - "In the days of its **invention by Kaiser** it was much easier to compute it than the discrete prolate spheroidal sequences discussed above. \n", + "In the days of its **invention by Kaiser**, it was much easier to compute it than the discrete prolate spheroidal sequences discussed above. \n", "This is due to the explicit given equation for the Kaiser-Bessel window, whereas for the Slepian window an eigenwert problem for a $M/2$ matrix has to be numerically solved.\n", "The Kaiser-Bessel window requires the [zeroth-order modified **Bessel** function of the first kind](https://dlmf.nist.gov/10.25) $I_0(\\cdot)$ to calculate the window signal $w[k]$.\n", "Thus, the given name in the DSP literature.\n", @@ -329,7 +329,7 @@ "source": [ "## Comparison of the Windows\n", "\n", - "What are the differences between the discussed **parametric** window types of same length and about the same level of the first sidelobe?" + "What are the differences between the discussed **parametric** window types of the same length and about the same level of the first sidelobe?" ] }, { @@ -468,7 +468,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.10.12" + "version": "3.12.2" } }, "nbformat": 4, diff --git a/dft/window_zplane_frequency_response.ipynb b/dft/window_zplane_frequency_response.ipynb index ed045a7..2f05b35 100644 --- a/dft/window_zplane_frequency_response.ipynb +++ b/dft/window_zplane_frequency_response.ipynb @@ -19,7 +19,7 @@ "- lecture: https://github.com/spatialaudio/digital-signal-processing-lecture\n", "- tutorial: https://github.com/spatialaudio/digital-signal-processing-exercises\n", "\n", - "Feel free to contact lecturer jacob.thoenes@uni-rostock.de" + "Feel free to contact the lecturer jacob.thoenes@uni-rostock.de" ] }, { @@ -162,9 +162,9 @@ "- This rational function has $N-1$ poles in the origin, i.e. $z=0$\n", "- Furthermore, it has $N-1$ zeros somewhere in the complex z-plane\n", "\n", - "For real valued $w[k]$ the zeros are either complex conjugate or real-valued (they lie on the real axis in the z-plane).\n", + "For real valued $w[k]$, the zeros are either complex conjugate or real-valued (they lie on the real axis in the z-plane).\n", "\n", - "For the windowing process we are interested in the DTFT spectrum, which can be obtained when evaluating $W(z)$ along the unit circle, thus as\n", + "For the windowing process, we are interested in the DTFT spectrum, which can be obtained when evaluating $W(z)$ along the unit circle, thus as\n", "\n", "\\begin{equation}\n", "W(\\Omega) = W(z)\\big|_{z=\\mathrm{e}^{\\mathrm{j}\\Omega}}\n", @@ -175,30 +175,30 @@ "Naturally, this design task has many, many solutions. Some of them are more meaningful than others, which is why they\n", "were published as special window designs over the last decades.\n", "\n", - "For windowing we actually always aim at very low sidelobes. Thus all potential zeros of $W(z)$ should used as best as possible: putting them exactly onto the unit circle the impact is largest in the DTFT spectrum, since notches are realized. Windows that follow this design rule, are commonly termed optimum windows. The design constraints for the trade-off then tell us where exactly to put these zeros on the unit circle." + "For windowing, we actually always aim at very low sidelobes. Thus, all potential zeros of $W(z)$ should be used as best as possible: putting them exactly onto the unit circle, the impact is largest in the DTFT spectrum, since notches are realized. Windows that follow this design rule are commonly termed optimum windows. The design constraints for the trade-off then tell us where exactly to put these zeros on the unit circle." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "In the examples below we discuss very simple cases for rather small $N$.\n", + "In the examples below, we discuss very simple cases for rather small $N$.\n", "\n", "**Example 1**\n", "\n", - "Although, is has worst leakage effect, the rectangular window is actually an optimum window. The design criterion is simply: provide the most narrow main lobe that is possible for $N$. The price you pay is the highest sidelobe magnitude and very slow sidelobe decay rate.\n", + "Although it has the worst leakage effect, the rectangular window is actually an optimum window. The design criterion is simply: provide the narrowest main lobe that is possible for $N$. The price you pay is the highest sidelobe magnitude and a very slow sidelobe decay rate.\n", "\n", "**Example 2**\n", "\n", - "This is a! symmetric Hann window (in older literature misnamed as Hanning window, don't do this anymore!). We see that two zeros are real valued and not on the unit circle. They don't shape the window spectrum in optimum manner.\n", + "This is a symmetric Hann window (in older literature misnamed as Hanning window, don't do this anymore!). We see that two zeros are real valued and not on the unit circle. They don't shape the window spectrum in an optimum manner.\n", "\n", - "Note that for didactical convenience we use a symmetric window here. For practical spectral analysis we rather should use the flag `sym=False` for scipy windows.\n", + "Note that for didactical convenience, we use a symmetric window here. For practical spectral analysis we rather should use the flag `sym=False` for scipy windows.\n", "\n", "**Example 3**\n", "\n", "Idea of Hamming window:\n", "- taking the above Hann window and moving the two zeros onto the unit circle (they must be conjugate complex).\n", - "- position of zeros such that they put a notch at the two sidelobes of the Hann window left and right of the main lobe\n", + "- position of zeros such that they put a notch at the two sidelobes of the Hann window, left and right of the main lobe\n", "\n", "Mainlobe shape remains, sidelobe magnitude is overall reduced compared to the Hann window." ] @@ -312,7 +312,7 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "We have two zeros that are **not** on the unit circle. If we put them **onto** the unit circle their influence with regard the DTFT (i.e. the frequency response) can be made stronger. By that we can optimize the window, for example to attenuate a certain sidelobe. This idea is pursued by the Hamming window, which is shown next." + "We have two zeros that are **not** on the unit circle. If we put them **onto** the unit circle, their influence with regard to the DTFT (i.e., the frequency response) can be made stronger. By that, we can optimize the window, for example, to attenuate a certain sidelobe. This idea is pursued by the Hamming window, which is shown next." ] }, { @@ -415,7 +415,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.10.12" + "version": "3.12.2" } }, "nbformat": 4, diff --git a/filter_design/filter_examples.ipynb b/filter_design/filter_examples.ipynb index 10e5639..4b9b51c 100644 --- a/filter_design/filter_examples.ipynb +++ b/filter_design/filter_examples.ipynb @@ -19,7 +19,7 @@ "- lecture: https://github.com/spatialaudio/digital-signal-processing-lecture\n", "- tutorial: https://github.com/spatialaudio/digital-signal-processing-exercises\n", "\n", - "Feel free to contact lecturer jacob.thoenes@uni-rostock.de" + "Feel free to contact the lecturer jacob.thoenes@uni-rostock.de" ] }, { @@ -280,7 +280,7 @@ "source": [ "## Linear Phase Finite Impulse Response (FIR) Filters\n", "\n", - "In the examples for type I, III, II and IV the coefficients in b are made up by random values. Thus, also the Level and Impulse Response are of arbitrary shape, which we are not really interested in.\n", + "In the examples for types I, III, II, and IV, the coefficients in b are made up of random values. Thus, the Level and Impulse Response are of arbitrary shape, which we are not really interested in.\n", "\n", "Rather, it is important to realize that the **group delay is constant over frequency** and the **impulse response** has **certain symmetry** over time for these linear phase filter types." ] @@ -417,7 +417,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.10.12" + "version": "3.12.2" } }, "nbformat": 4, diff --git a/filter_design/fir_filter.ipynb b/filter_design/fir_filter.ipynb index c8b19a8..ac89954 100644 --- a/filter_design/fir_filter.ipynb +++ b/filter_design/fir_filter.ipynb @@ -19,7 +19,7 @@ "- lecture: https://github.com/spatialaudio/digital-signal-processing-lecture\n", "- tutorial: https://github.com/spatialaudio/digital-signal-processing-exercises\n", "\n", - "Feel free to contact lecturer jacob.thoenes@uni-rostock.de" + "Feel free to contact the lecturer jacob.thoenes@uni-rostock.de" ] }, { @@ -219,7 +219,7 @@ "=\\frac{b_0z^0+b_1z^{-1}+b_2z^{-2}+...+b_Mz^{-M}}{a_0z^0+a_1z^{-1}+a_2z^{-2}+...+a_Nz^{-N}}\n", "\\end{equation}\n", "with input $X(z)$ and output $Y(z)$.\n", - "Real input signals $x[k]$ that should end up as real output signals $y[k]$ (in terms of signal processing fundamentals this is a special case, though most often needed in practice) require real coefficients $b,a\\in\\mathbb{R}$.\n", + "Real input signals $x[k]$ that should end up as real output signals $y[k]$ (in terms of signal processing fundamentals, this is a special case, though most often needed in practice) require real coefficients $b,a\\in\\mathbb{R}$.\n", "This is only achieved with\n", "\n", "- single or multiple **real** valued\n", @@ -227,15 +227,15 @@ "\n", "of zeros and poles.\n", "\n", - "Furthermore, in practice we most often aim at (i) causal and (ii) bounded input, bound output (BIBO) stable LTI systems, which requires (i) $M \\leq N$ and (ii) poles inside the unit circle.\n", - "If all poles **and** zeros are **inside** the unit circle then the system is **minimum-phase** and thus $H(z)$ is straightforwardly **invertible**.\n", + "Furthermore, in practice, we most often aim at (i) causal and (ii) bounded input, bound output (BIBO) stable LTI systems, which requires (i) $M \\leq N$ and (ii) poles inside the unit circle.\n", + "If all poles **and** zeros are **inside** the unit circle, then the system is **minimum-phase** and thus $H(z)$ is straightforwardly **invertible**.\n", "\n", "Further concepts related to the transfer function are:\n", "\n", "- Analysis of the transfer characteristics is done by the DTFT\n", "$H(z=\\mathrm{e}^{\\mathrm{j}\\Omega})$, i.e. evaluation on the unit circle.\n", "\n", - "- We use $a_0=1$ according to convention in many textbooks.\n", + "- We use $a_0=1$ according to the convention in many textbooks.\n", "\n", "- The convention for arraying filter coefficients is straightforward with Python index starting at zero:\n", "$b_0=b[0]$, $b_1=b[1]$, $b_2=b[2]$, ..., $a_0=a[0]=1$, $a_1=a[1]$, $a_2=a[2]$.\n", @@ -262,8 +262,8 @@ "\\end{equation}\n", "needs to be implemented.\n", "\n", - "- A **pure non-recursive** system is obtained by ignoring the feedback paths, i.e. setting $a_1,a_2,...,a_N=0$.\n", - "- A **pure recursive** system is obtained by ignoring the forward paths, i.e. setting $b_0,b_1,...,b_M=0$. Then, the values of the state variables $z^{-1}, z^{-2}, ..., z^{-M}$ alone determine how the system starts to perform at $k=0$, since the system has no input actually. This system type can be used to generate (damped) oscillations.\n", + "- A **pure non-recursive** system is obtained by ignoring the feedback paths, i.e., setting $a_1,a_2,...,a_N=0$.\n", + "- A **pure recursive** system is obtained by ignoring the forward paths, i.e., setting $b_0,b_1,...,b_M=0$. Then, the values of the state variables $z^{-1}, z^{-2}, ..., z^{-M}$ alone determine how the system starts to perform at $k=0$, since the system has no input actually. This system type can be used to generate (damped) oscillations.\n", "\n", "Please note: A recursive system can have a finite impulse response, but this is very rarely the case.\n", "Therefore, literature usually refers to\n", @@ -276,11 +276,11 @@ "- a non-recursive part (feedforward paths, left $z^{-}$-path)\n", "- a recursive part (feedback paths, left $z^{-}$-path)\n", "\n", - "is depicted below (graph taken from Wikimedia Commons) as straightforward **direct form I**, i.e. directly following the difference equation.\n", + "is depicted below (graph taken from Wikimedia Commons) as straightforward **direct form I**, i.e., directly following the difference equation.\n", "\n", "\n", "\n", - "Such as second order section is usually termed a biquad." + "Such as a second-order section is usually termed a biquad." ] }, { @@ -289,8 +289,8 @@ "source": [ "# FIR Filter\n", "\n", - "If all coefficients $a_{1,...,N}=0$, the feedback paths are not existent in the signal flow chart above.\n", - "This yields a non-recursive system and has transfer function\n", + "If all coefficients $a_{1,..., N}=0$, the feedback paths do not exist in the signal flow chart above.\n", + "This yields a non-recursive system and has a transfer function\n", "\n", "\\begin{equation}\n", "H(z) = \\frac{Y(z)}{X(z)} = \\sum\\limits_{m=0}^M b_mz^{-m}\n", @@ -303,7 +303,7 @@ "y[k] = b_0 x[k] + b_1 x[k-1] + b_2 x[k-2] + ... + b_M x[k-M],\n", "\\end{equation}\n", "\n", - "from which we can directly observe that the impulse response (i.e. for $x[k] = \\delta[k]$) is\n", + "from which we can directly observe that the impulse response (i.e., for $x[k] = \\delta[k]$) is\n", "\n", "\\begin{equation}\n", "h[k] = b_0 \\delta[k] + b_1 \\delta[k-1] + b_2 \\delta[k-2] + ... + b_M \\delta[k-M].\n", @@ -311,18 +311,17 @@ "\n", "This constitutes $h[k]$ as the coefficients $b_k$ at sample instances $k$.\n", "\n", - "The impulse response for this non-recursive system has always finite length of $M+1$ samples.\n", + "The impulse response for this non-recursive system always has a finite length of $M+1$ samples.\n", "\n", - "Usually this filter type is referred to as finite impulse response (FIR) filter in literature.\n", + "Usually, this filter type is referred to as a finite impulse response (FIR) filter in the literature.\n", "\n", - "Very special recursive systems/filters can produce FIRs as well. This is however so rare, that\n", - "the common link FIR filter == non-recursive system is predominantly made.\n", + "Very special recursive systems/filters can produce FIRs as well. This is, however, so rare that the common link FIR filter == non-recursive system is predominantly made.\n", "\n", - "The filter **order** is $M$, the **number of coefficients** $b$ is $M+1$. Be cautious here and consistent with the naming, it sometimes gets confusing. Especially for linear phase filters (see below) it is really important if $M$ or $M+1$ is either even or odd.\n", + "The filter **order** is $M$, the **number of coefficients** $b$ is $M+1$. Be cautious here and consistent with the naming, as it sometimes gets confusing. Especially for linear phase filters (see below), it is really important if $M$ or $M+1$ is either even or odd.\n", "\n", - "Sometimes the **number of taps** $M+1$ is stated (rather than calling this number of coefficients). This refers to tapping $M+1$ delayed instances of the signal input signal $x$ to calculate the filtered output signal $y$. Note however, that tapping the signal for the first coefficient $b_0$ involves non-delayed $x[k]$.\n", + "Sometimes the **number of taps** $M+1$ is stated (rather than calling this number of coefficients). This refers to tapping $M+1$ delayed instances of the input signal $x$ to calculate the filtered output signal $y$. Note, however, that tapping the signal for the first coefficient $b_0$ involves non-delayed $x[k]$.\n", "\n", - "For FIR filters the magnitude at $\\Omega=0$ and $\\Omega=\\pi$ can be straightforwardly evaluated with the following equations. \n", + "For FIR filters, the magnitude at $\\Omega=0$ and $\\Omega=\\pi$ can be straightforwardly evaluated with the following equations. \n", "\n", "The DC magnitude is obtained by \n", "\n", @@ -330,7 +329,7 @@ "g_0 = \\sum_{k=0}^{M} h[k].\n", "\\end{align}\n", "\n", - "The magnitude at $\\Omega=\\pi$ (i.e. at half the sampling frequency) is obtained by \n", + "The magnitude at $\\Omega=\\pi$ (i.e., at half the sampling frequency) is obtained by \n", "\n", "\\begin{align}\n", "g_\\pi = \\sum_{k=0}^{M} (-1)^k h[k].\n", @@ -343,7 +342,7 @@ "source": [ "## Poles / Zeros of FIR Filter\n", "\n", - "To calculate poles and zeros of the transfer function $H(z)$ it is meaningful to rewrite\n", + "To calculate poles and zeros of the transfer function $H(z)$, it is meaningful to rewrite\n", "\n", "\\begin{equation}\n", "H(z) = \\sum\\limits_{m=0}^M b_mz^{-m} \\frac{z^M}{z^M}\n", @@ -358,7 +357,7 @@ "\\sum\\limits_{m=0}^M b_m z^{-m} z^M = b_M z^M+b_1 z^{M-1}+b_2 z^{M-2}+...+b_M = 0\n", "\\end{equation}\n", "\n", - "needs to be solved. The $M$-th order polynomial has $M$ zeros. Recall from above, that for $b\\in\\mathbb{R}$ only real or complex conjugate zeros can occur, but never single complex zeros." + "needs to be solved. The $M$-th order polynomial has $M$ zeros. Recall from above that for $b\\in\\mathbb{R}$ only real or complex conjugate zeros can occur, but never single complex zeros." ] }, { @@ -367,14 +366,14 @@ "source": [ "## Essence of FIR Filter Design\n", "\n", - "The fundamental concept (just as was it with window design for the DFT) is to place the $M$ available zeros in the $z$-plane such that a target magnitude **and** phase response results, which suits the desired filter characteristics.\n", + "The fundamental concept (just as it was with window design for the DFT) is to place the $M$ available zeros in the $z$-plane such that a target magnitude **and** phase response results, which suits the desired filter characteristics.\n", "\n", "It is important to note that (contrary to IIR filters) the **magnitude** and **phase** response can be **separately controlled** with FIR filters.\n", "\n", - "In the referenced [english monographs](../index.ipynb) we can find information on specific FIR design techniques such as\n", + "In the referenced [english monographs](../index.ipynb), we can find information on specific FIR design techniques such as\n", "\n", - "- FIR design with windowing method\n", - "- FIR design as minimax least-squares optimization problem\n", + "- FIR design with the windowing method\n", + "- FIR design as a minimax least-squares optimization problem\n", "- FIR design of frequency sampling a DTFT spectrum\n", "\n", "These are well covered in Python's Scipy and Matlab. We will later discuss the windowing method. " @@ -388,7 +387,7 @@ "\n", "Consider an input signal $x[k]$.\n", "\n", - "An FIR filter $h[k]$ is used for convolution (i.e. filtering process)\n", + "An FIR filter $h[k]$ is used for convolution (i.e., filtering process)\n", "\n", "\\begin{align}\n", "x[k] \\ast h[k] \\circ-\\bullet X(\\Omega) \\cdot H(\\Omega),\n", @@ -400,9 +399,9 @@ "x[k] \\cdot w[k] \\circ-\\bullet \\frac{1}{2\\pi} X(\\Omega) \\circledast_{2\\pi} W(\\Omega).\n", "\\end{align}\n", "\n", - "In the DTFT domain this results in multiplication and circular convolution, respectively.\n", + "In the DTFT domain, this results in multiplication and circular convolution, respectively.\n", "\n", - "So, for the finite-length sequences $h[k]$ and $w[k]$ the same design fundamentals hold: we must put zeros at suitable locations in the z-plane to realize a certain desired DTFT spectrum, either $H(\\Omega)$ or $W(\\Omega)$. Depending on the application, filtering or windowing, the DTFT design criteria might be very different, since the DTFT spectrum acts as multiplication or convolution onto the input signal's DTFT spectrum. However, the design concepts and algorithms itself are basically the same, this is sometimes not so obvious in textbooks. " + "So, for the finite-length sequences $h[k]$ and $w[k]$, the same design fundamentals hold: we must put zeros at suitable locations in the z-plane to realize a certain desired DTFT spectrum, either $H(\\Omega)$ or $W(\\Omega)$. Depending on the application, filtering or windowing, the DTFT design criteria might be very different, since the DTFT spectrum acts as a multiplication or a convolution of the input signal's DTFT spectrum. However, the design concepts and algorithms themselves are basically the same, which is sometimes not so obvious in textbooks. " ] }, { @@ -412,14 +411,14 @@ "# FIR Examples with M=1 and M=2\n", "\n", "It is meaningful to discuss a very simple case of FIR in detail in the first place.\n", - "Once getting familiar with the underlying principles and concepts it is comparably easy to increase the FIR filter order and see how complexity evolves.\n", + "Once you get familiar with the underlying principles and concepts, it is comparatively easy to increase the FIR filter order and see how complexity evolves.\n", "\n", "Almost all important issues on FIRs can be explained with\n", "\n", "- the filter order $M=1$, thus number of coefficients is $M+1=2$ and with\n", "- the filter order $M=2$, thus number of coefficients is $M+1=3$.\n", "\n", - "The calculus of zeros is not too tedious and can be still performed manually. Furthermore, the impact of $M$ zeros in the $z$-plain is comparably easy linked to the magnitude and phase response.\n", + "The calculus of zeros is not too tedious and can still be performed manually. Furthermore, the impact of $M$ zeros in the $ z$-plane is comparably easily linked to the magnitude and phase response.\n", "\n", "So, let's play around with some simple coefficient settings." ] @@ -430,7 +429,7 @@ "source": [ "## Example FIR M=1, b0=1, b1=1\n", "\n", - "This is the most simple case of an FIR (actually the most simple FIR would be using only $b_0$, which is a pure gain/attenuation for input signal).\n", + "This is the most simple case of an FIR (actually the most simple FIR would be using only $b_0$, which is a pure gain/attenuation for the input signal).\n", "\n", "The squared magnitude response can be given analytically as (try yourself, make use of $\\mathrm{e}^{-\\mathrm{j}\\Omega}=\\cos(\\Omega)-\\mathrm{j}\\sin(\\Omega)$ ) \n", "\n", @@ -438,13 +437,13 @@ "|H(\\Omega)|^2 = |1 + \\mathrm{e}^{-\\mathrm{j}\\Omega}|^2 = 2 \\cos(\\Omega) + 2 = 4 \\cos^2(\\frac{\\Omega}{2}).\n", "\\end{equation}\n", "\n", - "Thus the magnitude response is \n", + "Thus, the magnitude response is \n", "\n", "\\begin{equation}\n", "|H(\\Omega)| = 2 |\\cos(\\frac{\\Omega}{2})|,\n", "\\end{equation}\n", "\n", - "which is confirmed by the below plot (left, top). The magnitude response exhibits lowpass characteristics.\n", + "which is confirmed by the plot below (left, top). The magnitude response exhibits lowpass characteristics.\n", "\n", "The impulse response is simply\n", "\n", @@ -454,7 +453,7 @@ "\n", "confirmed by the plot (right, bottom).\n", "\n", - "For the $M$-th order polynomial the $M=1$ zero is easy to evaluate\n", + "For the $M$-th order polynomial, the $M=1$ zero is easy to evaluate\n", "\n", "\\begin{equation}\n", "b_0 z^1+b_1 z^{0} = z+1 = 0\\to z_{0,1} = -1\n", @@ -462,7 +461,7 @@ "\n", "There is $M=1$ pole in the origin, i.e. $z_{\\infty,1} = 0$. Zero and pole are shown in the $z$-plane (right, top).\n", "\n", - "This FIR has special characteristics on the phase response, namely it is linear phase (type II), see discussion below.\n", + "This FIR has special characteristics on the phase response, namely, it is a linear phase (type II), see discussion below.\n", "\n", "In the present case DC magnitude is $g_0=2$ and $f_s/2$-magnitude is $g_\\pi=0$." ] @@ -634,7 +633,7 @@ "source": [ "## Example FIR M=2, b0=1, b1=-2, b2=1\n", "\n", - "By reversing the sign for $b_1$ compared to the just discussed lowpass we obtain a highpass.\n", + "By reversing the sign for $b_1$ compared to the just discussed lowpass, we obtain a highpass.\n", "\n", "Filter order $M=2$, number of coefficients $M+1=3$.\n", "\n", @@ -674,10 +673,10 @@ "source": [ "## Example FIR M=2, b0=1, b1=1, b2=1/2\n", "\n", - "So, far all zeros were aligned **on** the unit circle. We are not restricted to those locations, as long as we ensure positioning only real and complex-conjugate pairs in the $z$-plane.\n", + "So far, all zeros were aligned **on** the unit circle. We are not restricted to those locations, as long as we ensure positioning only real and complex-conjugate pairs in the $z$-plane.\n", "\n", - "However, recall that we discussed so called optimum window design for DFT-based spectral analysis.\n", - "There, zeros **on** the unit circle was a good idea, since the zero has most impact to shape the amplitude of the sidelobes. If a window has all zeros on the unit circle, we called this **optimum window**.\n", + "However, recall that we discussed the so called optimum window design for DFT-based spectral analysis.\n", + "There, zeros **on** the unit circle were a good idea, since the zero has the most impact to shape the amplitude of the sidelobes. If a window has all zeros on the unit circle, we call this **optimum window**.\n", "\n", "Back to our FIR design example:\n", "\n", @@ -690,9 +689,9 @@ "$g_0=\\frac{5}{2}$,\n", "$g_\\pi=\\frac{1}{2}$\n", "\n", - "This 2nd order FIR **lowpass** becomes a little more complicated and has a ripple in the stopband. Most important, since the zeros are **not** on the unit circle the magnitude response exhibits **no** exact zeros. Second morst important: the impulse response does not have a symmetry of linear-phase filter types I-IV, thus it is a **non-linear-phase** FIR. In fact, since poles and zeros are all within unit circle, the filter is minimum-phase.\n", + "This 2nd order FIR **lowpass** becomes a little more complicated and has a ripple in the stopband. Most importantly, since the zeros are **not** on the unit circle, the magnitude response exhibits **no** exact zeros. Second most important: the impulse response does not have a symmetry of linear-phase filter types I-IV, thus it is a **non-linear-phase** FIR. In fact, since poles and zeros are all within the unit circle, the filter is minimum-phase.\n", "\n", - "Try to find and plot the inverse transfer function $H(z)^{-1}$. Why this must be a minimum phase filter as well and why is this an IIR filter. The code in Jupyter notebook `iir_filter.ipynb` might be helpful." + "Try to find and plot the inverse transfer function $H(z)^{-1}$. Why must this be a minimum phase filter as well, and why is this an IIR filter? The code in Jupyter notebook `iir_filter.ipynb` might be helpful." ] }, { @@ -713,7 +712,7 @@ "source": [ "## Example FIR M=2, b0=1/2, b1=1, b2=1\n", "\n", - "We can even align zeros outside the unit circle, since contrary to poles this has no impact on stability.\n", + "We can even align zeros outside the unit circle, since, contrary to poles, this has no impact on stability.\n", "\n", "Filter order $M=2$, number of coefficients $M+1=3$.\n", "\n", @@ -724,9 +723,9 @@ "$g_0=\\frac{5}{2}$,\n", "$g_\\pi=\\frac{1}{2}$\n", "\n", - "This 2nd order FIR **lowpass** has the same magnitude response as the above lowpass. This is due to the fact that in both cases their zeros have same distance to the unit circle.\n", + "This 2nd order FIR **lowpass** has the same magnitude response as the above lowpass. This is due to the fact that in both cases their zeros have the same distance to the unit circle.\n", "\n", - "Obviously, this filter is also non-linear-phase, but also **not minimum-phase**. In fact, the filter is **maximum-phase**, i.e. the largest phase excess for the given magnitude response. Due to its maximum-phase, this FIR cannot simply be inverted, since then poles would lie outside the unit circle, yielding a non-stable system." + "Obviously, this filter is also non-linear-phase, but also **not minimum-phase**. In fact, the filter is **maximum-phase**, i.e., the largest phase excess for the given magnitude response. Due to its maximum-phase, this FIR cannot simply be inverted, since then the poles would lie outside the unit circle, yielding a non-stable system." ] }, { @@ -749,7 +748,7 @@ "source": [ "## Example Playground\n", "\n", - "Try some other coefficient settings by yourself, use larger $M$. Note that the plot of $h[k]$ is hard coded up to $k=7$ only.\n", + "Try some other coefficient settings by yourself, use a larger $M$. Note that the plot of $h[k]$ is hard coded up to $k=7$ only.\n", "\n", "For example, how about a lowpass filter where the coefficients come from the Kaiser-Bessel window?!?!\n", "\n", @@ -772,7 +771,7 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "Or another lowpass using the rectangular window. Yes, the rectangular window has lowpass characteristics, in literature this is known as the most simplest form of running average. In the example below, the last 7 samples are taken into account for averaging with equal weights." + "Or another lowpass using the rectangular window. Yes, the rectangular window has lowpass characteristics; in literature, this is known as the simplest form of running average. In the example below, the last 7 samples are taken into account for averaging with equal weights." ] }, { @@ -793,7 +792,7 @@ "source": [ "# FIR Filters with Linear Phase\n", "\n", - "In the above examples a very important concept of FIR filters was already included, which we did not discuss in detail so far. We do this here: the **linear phase**. In practice we can realize this only in digital signal processing, since analog circuits do not have the precision that would be required to shape the required impulse response characteristics (we talk about certain symmetries here as we will see below).\n", + "In the above examples, a very important concept of FIR filters was already included, which we did not discuss in detail so far. We do this here: the **linear phase**. In practice, we can realize this only in digital signal processing, since analog circuits do not have the precision that would be required to shape the required impulse response characteristics (we talk about certain symmetries here, as we will see below).\n", "\n", "A linear-phase FIR filter exhibits the DTFT spectrum\n", "\n", @@ -819,7 +818,7 @@ "- filter order $M$ even, odd number $M+1$ of filter coefficients $b$\n", "- even symmetry of the impulse response $h[k]=h[M-k]$\n", "- $\\beta=0,\\,\\pi$\n", - "- No fixed zeros, therefore all filter types are possible\n", + "- No fixed zeros, therefore, all filter types are possible\n", "\n", "## FIR Type II\n", "\n", @@ -857,10 +856,10 @@ "source": [ "# Windowed FIR Design Low-Pass\n", "\n", - "Now, let us discuss one potential FIR design method that is straightforward and comes typically first in teaching & learning DSP. It is still used in practice, especially when low computational complexity is aimed for.\n", - "The basic idea is to cut a suitable finite-length sequence out of an infinite impulse response of the ideal filter and to apply a window onto the finite-length sequence in order to reduce certain artifacts. The most simple case is using the rect window.\n", + "Now, let us discuss one potential FIR design method that is straightforward and typically comes first in teaching & learning DSP. It is still used in practice, especially when low computational complexity is aimed for.\n", + "The basic idea is to cut a suitable finite-length sequence out of an infinite impulse response of the ideal filter and to apply a window onto the finite-length sequence in order to reduce certain artifacts. The simplest case is using the rect window.\n", "\n", - "So, let's discuss the technique with the ideal lowpass filter: the ideal lowpass (i.e. zero-phase and infinite impulse response) with cutoff frequency $0 < \\Omega_c < \\pi$ is given by inverse DTFT as\n", + "So, let's discuss the technique with the ideal lowpass filter: the ideal lowpass (i.e., zero-phase and infinite impulse response) with cutoff frequency $0 < \\Omega_c < \\pi$ is given by inverse DTFT as\n", "\n", "\\begin{equation}\n", "h[k] = \\frac{1}{2 \\pi} \\int\\limits_{-\\Omega_c}^{+\\Omega_c} 1 \\cdot \\mathrm{e}^{+\\mathrm{j} \\Omega k} \\mathrm{d}\\Omega.\n", @@ -872,9 +871,9 @@ "h[k]=\\frac{\\sin\\left(\\Omega_c k \\right)}{\\pi k} = \\frac{\\Omega_c}{\\pi}\\frac{\\sin\\left(\\Omega_c k \\right)}{\\Omega_c k},\n", "\\end{equation}\n", "\n", - "i.e. a weighted sinc-function. We should have expected a sinc-like shape, since rect and sinc correspond in terms of the Fourier transform. \n", + "i.e., a weighted sinc-function. We should have expected a sinc-like shape, since rect and sinc correspond in terms of the Fourier transform. \n", "\n", - "In order to obtain a practical (finite order M) and causal (peak at M/2) we can shift this impulse response by $M/2$ (time delay)\n", + "In order to obtain a practical (finite order M) and causal (peak at M/2), we can shift this impulse response by $M/2$ (time delay)\n", "\n", "\\begin{equation}\n", "h[k]=\\frac{\\sin\\left(\\Omega_c\\left(k-\\frac{M}{2}\\right)\\right)}{\\pi\\left(k-\\frac{M}{2}\\right)}\n", @@ -889,12 +888,12 @@ "\\end{cases}\n", "\\end{equation}\n", "\n", - "The plain cut out (i.e. a rectangular window) towards a FIR yields\n", + "The plain cut out (i.e., a rectangular window) towards a FIR yields\n", "\n", "- on the one hand, the steepest roll-off for $\\Omega>\\Omega_c$ towards the first zero in the magnitude response\n", - "- but on the other hand also the worst stop band damping in the magnitude response (maximum stopband level is only about -21 dB)\n", + "- but on the other hand, also the worst stop band damping in the magnitude response (maximum stopband level is only about -21 dB)\n", "\n", - "Most often, larger maximum stopband level is desired under acceptance of less steep initial roll-off. This can be achieved with other window functions. A very well suited window for this filter design task is the Kaiser-Bessel window. Kaiser figured out that a certain maximum stopband level (usually this is here the first side lobe) can be controlled by the parameter $\\beta$. This of course holds only if $M$ is chosen large enough to obtain such a damping at all.\n", + "Most often, a larger maximum stopband level is desired under the acceptance of a less steep initial roll-off. This can be achieved with other window functions. A very well suited window for this filter design task is the Kaiser-Bessel window. Kaiser figured out that a certain maximum stopband level (usually this is the first side lobe) can be controlled by the parameter $\\beta$. This, of course, holds only if $M$ is chosen large enough to obtain such a damping at all.\n", "\n", "For a maximum stopband level $\\gamma<-50$ in dB, the approximation\n", "\n", @@ -912,7 +911,7 @@ "\n", "hence the design method's naming.\n", "\n", - "In the examples below we use a Kaiser-Bessel window to control that stopband level does not exceed -54 dB.\n", + "In the examples below, we use a Kaiser-Bessel window to ensure that the stopband level does not exceed -54 dB.\n", "\n", "We design linear-phase type I FIR filters with the windowing method. \n", "\n", @@ -962,9 +961,9 @@ "source": [ "# Windowed FIR Design High-Pass\n", "\n", - "We can do the same approach for the ideal highpass filter:\n", + "We can do the same approach for the ideal high-pass filter:\n", "\n", - "The ideal highpass (i.e. zero-phase and infinite impulse response) with cutoff frequency $0 < \\Omega_c < \\pi$ is given by inverse DTFT as\n", + "The ideal highpass (i.e,. zero-phase and infinite impulse response) with cutoff frequency $0 < \\Omega_c < \\pi$ is given by inverse DTFT as\n", "\n", "\\begin{equation}\n", "h[k] = \\frac{1}{2 \\pi} \\int\\limits_{\\Omega_c}^{2\\pi-\\Omega_c} 1 \\cdot \\mathrm{e}^{+\\mathrm{j} \\Omega k} \\mathrm{d}\\Omega.\n", @@ -976,7 +975,7 @@ "h[k]=\\frac{\\sin\\left(\\pi k \\right)}{\\pi k}-\\frac{\\sin\\left(\\Omega_c k\\right)}{\\pi k}.\n", "\\end{equation}\n", "\n", - "In order to obtain a practical (finite order M) and causal (peak at M/2) we can shift this impulse response by $M/2$ (time delay)\n", + "In order to obtain a practical (finite order M) and causal (peak at M/2), we can shift this impulse response by $M/2$ (time delay)\n", "\n", "\\begin{equation}\n", "h[k]=\\frac{\\sin\\left(\\pi\\left(k-\\frac{M}{2}\\right)\\right)}{\\pi\\left(k-\\frac{M}{2}\\right)}-\\frac{\\sin\\left(\\Omega_c\\left(k-\\frac{M}{2}\\right)\\right)}{\\pi\\left(k-\\frac{M}{2}\\right)}\n", @@ -997,7 +996,7 @@ "h_w[k] = w[k] \\cdot h[k].\n", "\\end{equation}\n", "\n", - "Kaiser-Bessel window is also perfectly suited for the highpass. So we use precisely the same window as for the lowpass." + "The Kaiser-Bessel window is also perfectly suited for the highpass. So we use precisely the same window as for the lowpass." ] }, { @@ -1017,7 +1016,7 @@ "source": [ "Play around with $M$, $\\Omega_c$, `StopBandMaxLevel` ($\\beta$) and how they are linked. Then, you're ready to study FIR design methods in books in detail.\n", "\n", - "It is also worth to check the pole/zero plot of these lowpass / highpass filters." + "It is also worth checking the pole/zero plot of these lowpass/highpass filters." ] }, { diff --git a/filter_design/iir_biquad.ipynb b/filter_design/iir_biquad.ipynb index f4020ca..cdbd288 100644 --- a/filter_design/iir_biquad.ipynb +++ b/filter_design/iir_biquad.ipynb @@ -19,7 +19,7 @@ "- lecture: https://github.com/spatialaudio/digital-signal-processing-lecture\n", "- tutorial: https://github.com/spatialaudio/digital-signal-processing-exercises\n", "\n", - "Feel free to contact lecturer jacob.thoenes@uni-rostock.de" + "Feel free to contact the lecturer jacob.thoenes@uni-rostock.de" ] }, { @@ -28,8 +28,8 @@ "metadata": {}, "outputs": [], "source": [ - "import numpy as np\n", "import matplotlib.pyplot as plt\n", + "import numpy as np\n", "from matplotlib.markers import MarkerStyle\n", "from matplotlib.patches import Circle\n", "from scipy import signal\n", @@ -100,8 +100,8 @@ "\\label{eq:biquad}\n", "\\end{equation}\n", "\n", - "- The notation in form of the **first fraction** in above equation is suitable for identifying the structure of the filter (the necessary delay elements for the input (numerator) and the output (denominator) can directly be seen) and for creating block diagrams (direct form I/II, transposed direct forms, ...) from the difference equation.\n", - "- The form of the **second fraction** is more convenient for calculating poles, zeros and gain of the filter." + "- The notation in the form of the **first fraction** in the above equation is suitable for identifying the structure of the filter (the necessary delay elements for the input (numerator) and the output (denominator) can directly be seen) and for creating block diagrams (direct form I/II, transposed direct forms, ...) from the difference equation.\n", + "- The form of the **second fraction** is more convenient for calculating poles, zeros, and gain of the filter." ] }, { @@ -127,7 +127,7 @@ "source": [ "## Real Poles\n", "\n", - "For **real poles** (i.e. poles on the real-axis in the z-plane), the expression under the square root has to be positive\n", + "For **real poles** (i.e., poles on the real axis in the z-plane), the expression under the square root has to be positive\n", "\n", "\\begin{equation}\n", "z_{\\infty,1,2}=\\frac{-a_1}{2}\\pm\\frac{1}{2}\\sqrt{\\underbrace{a_1^2-4a_2}_{\\geq0}}\n", @@ -140,7 +140,7 @@ "a_2>|a_1|-1\n", "\\end{align}\n", "\n", - "must hold (try to figure out these conditions by yourself) to achieve stable IIR system with real valued coefficients $a_{1,2}$." + "must hold (try to figure out these conditions by yourself) to achieve a stable IIR system with real-valued coefficients $a_{1,2}$." ] }, { @@ -260,7 +260,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.10.12" + "version": "3.12.2" } }, "nbformat": 4, diff --git a/filter_design/iir_filter.ipynb b/filter_design/iir_filter.ipynb index 97a7719..2617d9d 100644 --- a/filter_design/iir_filter.ipynb +++ b/filter_design/iir_filter.ipynb @@ -19,7 +19,7 @@ "- lecture: https://github.com/spatialaudio/digital-signal-processing-lecture\n", "- tutorial: https://github.com/spatialaudio/digital-signal-processing-exercises\n", "\n", - "Feel free to contact lecturer jacob.thoenes@uni-rostock.de" + "Feel free to contact the lecturer jacob.thoenes@uni-rostock.de" ] }, { @@ -224,7 +224,7 @@ "=\\frac{b_0z^0+b_1z^{-1}+b_2z^{-2}+...+b_Mz^{-M}}{a_0z^0+a_1z^{-1}+a_2z^{-2}+...+a_Nz^{-N}}\n", "\\end{equation}\n", "with input $X(z)$ and output $Y(z)$.\n", - "Real input signals $x[k]$ that should end up as real output signals $y[k]$ (in terms of signal processing fundamentals this is a special case, though most often needed in practice) require real coefficients $b,a\\in\\mathbb{R}$.\n", + "Real input signals $x[k]$ that should end up as real output signals $y[k]$ (in terms of signal processing fundamentals, this is a special case, though most often needed in practice) require real coefficients $b,a\\in\\mathbb{R}$.\n", "This is only achieved with\n", "\n", "- single or multiple **real** valued\n", @@ -232,15 +232,15 @@ "\n", "of zeros and poles.\n", "\n", - "Furthermore, in practice we most often aim at (i) causal and (ii) bounded input, bound output (BIBO) stable LTI systems, which requires (i) $M \\leq N$ and (ii) poles inside the unit circle.\n", - "If all poles **and** zeros are **inside** the unit circle then the system is **minimum-phase** and thus $H(z)$ is straightforwardly **invertible**.\n", + "Furthermore, in practice, we most often aim at (i) causal and (ii) bounded input, bound output (BIBO) stable LTI systems, which requires (i) $M \\leq N$ and (ii) poles inside the unit circle.\n", + "If all poles **and** zeros are **inside** the unit circle, then the system is **minimum-phase** and thus $H(z)$ is straightforwardly **invertible**.\n", "\n", "Further concepts related to the transfer function are:\n", "\n", "- Analysis of the transfer characteristics is done by the DTFT\n", "$H(z=\\mathrm{e}^{\\mathrm{j}\\Omega})$, i.e. evaluation on the unit circle.\n", "\n", - "- We use $a_0=1$ according to convention in many textbooks.\n", + "- We use $a_0=1$ according to the convention in many textbooks.\n", "\n", "- The convention for arraying filter coefficients is straightforward with Python index starting at zero:\n", "$b_0=b[0]$, $b_1=b[1]$, $b_2=b[2]$, ..., $a_0=a[0]=1$, $a_1=a[1]$, $a_2=a[2]$.\n", @@ -267,8 +267,8 @@ "\\end{equation}\n", "needs to be implemented.\n", "\n", - "- A **pure non-recursive** system is obtained by ignoring the feedback paths, i.e. setting $a_1,a_2,...,a_N=0$.\n", - "- A **pure recursive** system is obtained by ignoring the forward paths, i.e. setting $b_0,b_1,...,b_M=0$. Then, the values of the state variables $z^{-1}, z^{-2}, ..., z^{-M}$ alone determine how the system starts to perform at $k=0$, since the system has no input actually. This system type can be used to generate (damped) oscillations.\n", + "- A **pure non-recursive** system is obtained by ignoring the feedback paths, i.e., setting $a_1,a_2,...,a_N=0$.\n", + "- A **pure recursive** system is obtained by ignoring the forward paths, i.e,. setting $b_0,b_1,...,b_M=0$. Then, the values of the state variables $z^{-1}, z^{-2}, ..., z^{-M}$ alone determine how the system starts to perform at $k=0$, since the system has no input actually. This system type can be used to generate (damped) oscillations.\n", "\n", "Please note: A recursive system can have a finite impulse response, but this is very rarely the case.\n", "Therefore, literature usually refers to\n", @@ -281,7 +281,7 @@ "- a non-recursive part (feedforward paths, left $z^{-}$-path)\n", "- a recursive part (feedback paths, left $z^{-}$-path)\n", "\n", - "is depicted below (graph taken from Wikimedia Commons) as straightforward **direct form I**, i.e. directly following the difference equation.\n", + "is depicted below (graph taken from Wikimedia Commons) as straightforward **direct form I**, i.e., directly following the difference equation.\n", "\n", "\n", "\n", @@ -294,8 +294,8 @@ "source": [ "# IIR Filter\n", "\n", - "In the following 2nd order IIR filters shall be discussed.\n", - "The transfer function is with usual convention $a_0=1$\n", + "In the following, 2nd order IIR filters shall be discussed.\n", + "The transfer function is with the usual convention $a_0=1$\n", "\n", "\\begin{equation}\n", "H(z)=\\frac{\\sum\\limits_{m=0}^2 b_mz^{-m}}{\\sum\\limits_{n=0}^2 a_nz^{-n}}\n", @@ -411,13 +411,13 @@ "\n", "## Mapping\n", "\n", - "A Laplace transfer function can be transformed to $z$-domain transfer function by \n", + "A Laplace transfer function can be transformed to a $z$-domain transfer function by \n", "\n", "\\begin{equation}\n", "s=2f_s\\frac{z-1}{z+1}.\n", "\\end{equation}\n", "\n", - "This is known as bilinear transform.\n", + "This is known as the bilinear transform.\n", "There are several derivations for this mapping. One way is to linearize a function, namely applying the Taylor series and then taking only the linear term:\n", "So, we know that the exact link between $s$ and $z$ is given by the definition\n", "\\begin{equation}\n", @@ -445,13 +445,13 @@ "- $s=0$ is mapped to $z=1$\n", "- infinite, positive $\\mathrm{j}\\omega$ axis of $s$-domain is mapped to the upper unit circle\n", "- infinite, negative $\\mathrm{j}\\omega$ axis of $s$-domain is mapped to the lower unit circle\n", - "- $n$-th order filter in $s$-domain yields also $n$-th filter in $z$-domain, due to same powers of $s$ and $z$ in the mapping. Question: what happens if one uses more terms of the Taylor approximation.\n", + "- $n$-th order filter in $s$-domain yields also $n$-th filter in $z$-domain, due to same powers of $s$ and $z$ in the mapping. Question: What happens if one uses more terms of the Taylor approximation?\n", "\n", "When rearranging the mapping to \n", "\\begin{equation}\n", "z = \\frac{2 + s T}{2 - s T}\n", "\\end{equation}\n", - "we can see, what happens for $s = \\pm \\mathrm{j}\\omega$ with $\\omega\\rightarrow\\infty$. In both cases the mapping $z=-1$ is achieved. This make the upper/lower semi-circle obvious." + "we can see, what happens for $s = \\pm \\mathrm{j}\\omega$ with $\\omega\\rightarrow\\infty$. In both cases, the mapping $z=-1$ is achieved. This makes the upper/lower semi-circle obvious." ] }, { @@ -488,16 +488,16 @@ "For the example given below (using $f_s=1$ Hz), consider the **tan-mapping** in the left plot:\n", "\n", "A digital bandwidth from $0 < \\Omega < 1$ maps approximately to $0 < \\omega < 1$, the digital bandwidth\n", - "$1 < \\Omega < 2$ however maps to $1 < \\omega < 3$ in analog domain. For digital bandwidth $2 < \\Omega < 3$ this bandwidth mismatch gets even worse and is not depicted in the graph anymore.\n", + "$1 < \\Omega < 2$ however maps to $1 < \\omega < 3$ in analog domain. For digital bandwidth $2 < \\Omega < 3$, this bandwidth mismatch gets even worse and is not depicted in the graph anymore.\n", "\n", - "Digital to analog bandwidth mapping: the higher the digital frequency, the larger the bandwidth mismatch compared to analog domain becomes, analog bandwidth appears larger.\n", + "Digital to analog bandwidth mapping: the higher the digital frequency, the larger the bandwidth mismatch compared to the analog domain becomes, and analog bandwidth appears larger.\n", "\n", "Now, consider the **atan-mapping** in the right plot:\n", "\n", "An analog bandwidth from $0<\\omega<1$ maps approximately to $0<\\Omega<1$, we should expect this, since this is the region where tan() and atan() are approximately linear.\n", "However, an analog bandwidth from $1<\\omega<2$ maps approximately to $1<\\Omega<\\frac{3}{2}$, meaning that the digital bandwidth gets squeezed.\n", "\n", - "Thus, analog to digital bandwidth mapping: the higher the analog frequency, the larger the bandwidth mismatch compared to digital domain becomes, digital bandwidth appears smaller." + "Thus, analog to digital bandwidth mapping: the higher the analog frequency, the larger the bandwidth mismatch compared to the digital domain becomes, and digital bandwidth appears smaller." ] }, { @@ -557,13 +557,13 @@ "metadata": {}, "source": [ "There is no way to avoid these non-linear mapping effects, they come inherently with the application of the bilinear transform.\n", - "However, one can be cautious in filter design basically having three options:\n", + "However, one can be cautious in filter design, basically having three options:\n", "\n", - "1. Make sure that the bandwidth(s) under discussion are far away from half of the sampling frequency. Or in other words, try to work in the linear-like part of the tan/atan-mapping functions\n", - "2. Since 1. cannot be fulfilled always, one could select exactly one frequency where the non-linear mapping could be compensated for in the frequency design. A meaningful choice is the cut-off frequency of low- or highpass filters or the mid-frequency of bandpass or -stop filters. This technique is known as frequency pre-warping.\n", - "3. For bandpass/stop filters, which have a Q-factor / bandwidth parameter, one could do a bandwidth pre-warping as well. However, compared to option 2 (which is exact), this effect works only as an approximation. It can never compensate the mapping mismatch completely.\n", + "1. Make sure that the bandwidth(s) under discussion are far away from half of the sampling frequency. Or, in other words, try to work in the linear-like part of the tan/atan-mapping functions\n", + "2. Since 1. cannot be fulfilled always, one could select exactly one frequency where the non-linear mapping could be compensated for in the frequency design. A meaningful choice is the cut-off frequency of low- or high-pass filters or the mid-frequency of band-pass or -stop filters. This technique is known as frequency pre-warping.\n", + "3. For bandpass/stop filters, which have a Q-factor / bandwidth parameter, one could do a bandwidth pre-warping as well. However, compared to option 2 (which is exact), this effect works only as an approximation. It can never compensate for the mapping mismatch completely.\n", "\n", - "Let us first check one of the most simplest examples: the analog lowpass filter of first order and the design of discrete-time version with bilinear transform, in order to check the effects." + "Let us first check one of the simplest examples: the analog lowpass filter of first order and the design of a discrete-time version with bilinear transform, in order to check the effects." ] }, { @@ -589,15 +589,15 @@ "\\begin{equation}\n", "H(z) = \\frac{\\frac{1}{1+\\frac{2 f_s}{\\omega_c}}+\\frac{1}{1+\\frac{2 f_s}{\\omega_c}}z^{-1}}{1+\\frac{1-\\frac{2 f_s}{\\omega_c}}{1+\\frac{2 f_s}{\\omega_c}} z^{-1}} = \\frac{b_0 + b_1 z^{-1}}{1 + a_1 z^{-1}},\n", "\\end{equation}\n", - "and thus a first order filter.\n", + "and thus a first-order filter.\n", "\n", - "The single, real zero is always given at $z_0 = -\\frac{b_1}{b_0} = -1$ independently from the chosen sampling frequency and cut-off-frequency. This is not by accident: since, a lowpass (in other words a high-cut) is intended, the bilinear transform forces the amplitude to zero at highest usable frequency, which is $f_s/2$. The single, real pole is given at $z_\\infty = -a_1$. It is simple to show, that for $f_c = 0$, the pole is at $z_\\infty = 1$. For $f_c = f_s/2$, the pole is at $z_\\infty = - \\frac{1-2/\\pi}{1+2/\\pi} \\approx -0.22203094070331453$.\n", - "Thus, for $f_c \\ll f_s/2$ the pole is close to the unit circle on the $\\Re(z)$-axis and then moving to left along the axis when increasing $f_c$.\n", - "From this simple example (and although not proving it), we see that a stable Laplace transfer function yields a stable discrete-time version using bilinear transform. The pole is always within the unit circle (besides the special case of $f_c=0$ Hz).\n", + "The single, real zero is always given at $z_0 = -\\frac{b_1}{b_0} = -1$ independently of the chosen sampling frequency and cut-off frequency. This is not by accident: since a lowpass (in other words, a high-cut) is intended, the bilinear transform forces the amplitude to zero at the highest usable frequency, which is $f_s/2$. The single, real pole is given at $z_\\infty = -a_1$. It is simple to show that for $f_c = 0$, the pole is at $z_\\infty = 1$. For $f_c = f_s/2$, the pole is at $z_\\infty = - \\frac{1-2/\\pi}{1+2/\\pi} \\approx -0.22203094070331453$.\n", + "Thus, for $f_c \\ll f_s/2$ the pole is close to the unit circle on the $\\Re(z)$-axis and then moving to the left along the axis when increasing $f_c$.\n", + "From this simple example (and although not proving it), we see that a stable Laplace transfer function yields a stable discrete-time version using the bilinear transform. The pole is always within the unit circle (besides the special case of $f_c=0$ Hz).\n", "\n", "Let's implement this little example to compare the analog and discrete-time level responses over frequency.\n", "\n", - "Play around with the cut-off frequency $\\omega_c$. The closer it gets to $f_s/2$ the more the level responses deviate. Can you explain this? Pay special attention to the lowpass slope and where the cut-off frequency of the digital filter is." + "Play around with the cut-off frequency $\\omega_c$. The closer it gets to $f_s/2$, the more the level responses deviate. Can you explain this? Pay special attention to the lowpass slope and where the cut-off frequency of the digital filter is." ] }, { @@ -653,7 +653,7 @@ "\n", "Linear mapping of $\\omega_c \\rightarrow \\Omega_c$ gives $\\Omega_c = \\frac{\\omega_c}{f_s}$.\n", "\n", - "This $\\Omega_c$ is now inserted to the mapping \n", + "This $\\Omega_c$ is now inserted into the mapping \n", "\n", "\\begin{equation}\n", "\\omega_{c,\\mathrm{PRE}} = 2 f_s \\tan{\\frac{\\Omega_c}{2}} = 2 f_s \\tan{\\frac{2\\pi f_c}{2 f_s}}.\n", @@ -665,7 +665,7 @@ "\\Omega_c = 2 \\mathrm{atan}{\\frac{\\omega_{c,\\mathrm{PRE}}}{2 f_s}}\n", "\\end{equation}\n", "\n", - "holds we when remap and thus we obtain our desired $\\Omega_c$.\n", + "holds when we remap, and thus we obtain our desired $\\Omega_c$.\n", "\n", "Thus, instead of using $\\omega_c$ in the Laplace transfer function, we use the so called **pre-warped** cut-off frequency $\\omega_{c,\\mathrm{PRE}}$.\n", "\n", @@ -745,18 +745,18 @@ "source": [ "## General Second Order Filter Bilinear Transform\n", "\n", - "The second order filter is very important in practical filter design, since higher-order filters are usually split into a cascade of second order sections (SOS). Thus, it is meaningful to explicitly derive the bilinear transform of a general SOS given in Laplace domain. \n", + "The second order filter is very important in practical filter design, since higher-order filters are usually split into a cascade of second order sections (SOS). Thus, it is meaningful to explicitly derive the bilinear transform of a general SOS given in the Laplace domain. \n", "\n", - "Let us denote the Laplace transfer function with real coefficients $B,A$ as\n", + "Let us denote the Laplace transfer function with real coefficients $B, A$ as\n", "\\begin{equation}\n", "H(s)=\\frac{B_0s^2+B_1s+B_2}{A_0s^2+A_1s+A_2}\n", "\\end{equation}\n", - "for a second order filter. This transfer function is to be transformed to a discrete-time filter via bilinear transform. Let us denote the $z$-transfer function with real coefficients $b,a$ as\n", + "for a second order filter. This transfer function is to be transformed to a discrete-time filter via the bilinear transform. Let us denote the $z$-transfer function with real coefficients $b,a$ as\n", "\\begin{equation}\n", "H(z)=\\frac{b_0+b_1z^{-1}+b_2z^{-2}}{1+a_1z^{-1}+a_2z^{-2}}.\n", "\\end{equation}\n", "\n", - "The aboved discussed mapping $s=2f_s\\frac{z-1}{z+1}$ can be inserted to $H(s)$ to yield a representation for $H(z)$. The coefficients for this $H(z)$ can be derived to\n", + "The above discussed mapping $s=2f_s\\frac{z-1}{z+1}$ can be inserted into $H(s)$ to yield a representation for $H(z)$. The coefficients for this $H(z)$ can be derived to\n", "\\begin{equation}\n", "\\begin{split}\n", "b_0=\\frac{B_2+2B_1f_s+4B_0f_s^2}{A_2+2A_1f_s+4A_0f_s^2}\\\\\n", @@ -801,7 +801,7 @@ "\\end{equation}\n", "and plot its magnitude in dB over frequency\n", "\n", - "5. Design $H_3(z)$ with bilinear transform, frequency pre-warping and additional so-called Q-factor pre-warping (here with the tan()-function, there are other mappings as well, see [audio filters](https://nbviewer.jupyter.org/github/spatialaudio/digital-signal-processing-lecture/blob/master/filter_design/audiofilter.ipynb))\n", + "5. Design $H_3(z)$ with bilinear transform, frequency pre-warping, and additional so-called Q-factor pre-warping (here with the tan()-function, there are other mappings as well, see [audio filters](https://nbviewer.jupyter.org/github/spatialaudio/digital-signal-processing-lecture/blob/master/filter_design/audiofilter.ipynb))\n", "\\begin{equation}\n", "2f_s\\tan{\\left(\\pi\\frac{f_0}{f_s}\\right)} \\rightarrow \\omega_{0}\n", "\\end{equation}\n", @@ -810,7 +810,7 @@ "\\end{equation}\n", "Plot its magnitude in dB over frequency.\n", "\n", - "6. Discuss the differences of the level responses. Especially pay attention to what happens at $f_0$, $f_s/2$ and with the bandwidth of the filter.\n", + "6. Discuss the differences in the level responses. Especially pay attention to what happens at $f_0$, $f_s/2$, and with the bandwidth of the filter.\n", "\n", "Ife02...Ifeachor, E.C.; Jervis, B.W. (2002): Digital Signal Processing. 2nd ed. Prentice Hall." ] @@ -917,7 +917,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.10.12" + "version": "3.12.2" } }, "nbformat": 4, diff --git a/fundamentals/linear_algebra.ipynb b/fundamentals/linear_algebra.ipynb index 9f4856e..c3a7753 100644 --- a/fundamentals/linear_algebra.ipynb +++ b/fundamentals/linear_algebra.ipynb @@ -20,7 +20,7 @@ "- lecture: https://github.com/spatialaudio/digital-signal-processing-lecture\n", "- tutorial: https://github.com/spatialaudio/digital-signal-processing-exercises\n", "\n", - "Feel free to contact lecturer jacob.thoenes@uni-rostock.de" + "Feel free to contact the lecturer jacob.thoenes@uni-rostock.de" ] }, { @@ -60,7 +60,7 @@ "### Matrix and Vector\n", "\n", "A matrix is a two-dimensional arrangement of numbers, expressed in rows and columns. \n", - "Let's assume very simple matrix $\\textbf{A}$ with a row dimension $m$ and comulmn dimension $n$.\n", + "Let's assume a very simple matrix $\\textbf{A}$ with a row dimension $m$ and column dimension $n$.\n", "\n", "$$\n", "\\boldsymbol {A}=\\begin{pmatrix}a_{11}&a_{12}&\\cdots &a_{1n}\\\\a_{21}&a_{22}&\\cdots &a_{2n}\\\\\\vdots &\\vdots &&\\vdots \\\\a_{m1}&a_{m2}&\\cdots &a_{mn}\\\\\\end{pmatrix}\n", @@ -287,8 +287,8 @@ "\\end{matrix}\n", "$$\n", "### Using Gaussian Elimination\n", - "If we want to solve this system of equations we could choose the [gaussian elimination](https://en.wikipedia.org/wiki/Gaussian_elimination#Example_of_the_algorithm) method.\n", - "Therefore we rewrite it for example like:\n", + "If we want to solve this system of equations, we could choose the [Gaussian elimination](https://en.wikipedia.org/wiki/Gaussian_elimination#Example_of_the_algorithm) method.\n", + "Therefore, we rewrite it, for example, like:\n", "\n", "\\begin{bmatrix}\n", "\\begin{array}{ccc|c}\n", @@ -361,7 +361,7 @@ "id": "b5fba7f9-0dcb-4e8e-885c-b1b752fc1381", "metadata": {}, "source": [ - "Matrix $\\boldsymbol{A}$ is an $m \\times n$ matrix, with $m = n$. The determinant is not zero, $\\boldsymbol{A}^\\text{T}\\boldsymbol{A} = \\boldsymbol{I}$, and the matrix has full rank because $\\text{rank}(\\boldsymbol{A}) = m = n$. So finally the matrix $\\boldsymbol{A}$ satisfies the conditions for being invertible and we can define $\\boldsymbol{b}$ and solve for $\\boldsymbol{x}$." + "Matrix $\\boldsymbol{A}$ is an $m \\times n$ matrix, with $m = n$. The determinant is not zero, $\\boldsymbol{A}^\\text{T}\\boldsymbol{A} = \\boldsymbol{I}$, and the matrix has full rank because $\\text{rank}(\\boldsymbol{A}) = m = n$. So finally the matrix $\\boldsymbol{A}$ satisfies the conditions for being invertible, and we can define $\\boldsymbol{b}$ and solve for $\\boldsymbol{x}$." ] }, { @@ -417,7 +417,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.10.12" + "version": "3.12.2" } }, "nbformat": 4, diff --git a/fundamentals/motivation_application_concepts.ipynb b/fundamentals/motivation_application_concepts.ipynb index 163195f..b17a888 100755 --- a/fundamentals/motivation_application_concepts.ipynb +++ b/fundamentals/motivation_application_concepts.ipynb @@ -20,7 +20,7 @@ "- lecture: https://github.com/spatialaudio/digital-signal-processing-lecture\n", "- tutorial: https://github.com/spatialaudio/digital-signal-processing-exercises\n", "\n", - "Feel free to contact lecturer jacob.thoenes@uni-rostock.de" + "Feel free to contact the lecturer jacob.thoenes@uni-rostock.de" ] }, { @@ -30,7 +30,7 @@ "source": [ "## 1. Motivation\n", "\n", - "Digital Signal Processing (DSP) is a rapidly growing field with numerous applications in modern technology and communication systems. The ability to process, analyze and manipulate signals in the digital domain has revolutionized fields such as audio and speech processing, image and video processing, communications, control systems, and biomedical engineering. Learning DSP enables one to understand the theory and techniques used to, for example process signals in real-time, which has a wide range of applications in areas such as speech recognition, noise reduction, image compression, and many others. Hence, it is a valuable skill for individuals interested in careers in electronics, telecommunications, and computer engineering.\n", + "Digital Signal Processing (DSP) is a rapidly growing field with numerous applications in modern technology and communication systems. The ability to process, analyze, and manipulate signals in the digital domain has revolutionized fields such as audio and speech processing, image and video processing, communications, control systems, and biomedical engineering. Learning DSP enables one to understand the theory and techniques used to, for example, process signals in real-time, which has a wide range of applications in areas such as speech recognition, noise reduction, image compression, and many others. Hence, it is a valuable skill for individuals interested in careers in electronics, telecommunications, and computer engineering.\n", "\n", "### Advantages of DSP\n", "_based on: [Holton,2021]_\n", @@ -101,7 +101,7 @@ "- Signal modification\n", "- Realization of recursive and non-recursive filter\n", "- Familiarity with signal enhancement techniques such as noise reduction.\n", - "- Fundamentals of designing and implementation of DSP algorithms and systems.\n", + "- Fundamentals of designing and implementing DSP algorithms and systems.\n", "\n", "The ultimate goal is to be able to analyze, design, and implement DSP systems and algorithms for various applications, and to understand the trade-offs and limitations involved in real-world DSP systems." ] @@ -133,7 +133,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.10.12" + "version": "3.12.2" } }, "nbformat": 4, diff --git a/index.ipynb b/index.ipynb index c44fb11..4cde789 100644 --- a/index.ipynb +++ b/index.ipynb @@ -18,7 +18,7 @@ "- lecture: https://github.com/spatialaudio/digital-signal-processing-lecture\n", "- tutorial: https://github.com/spatialaudio/digital-signal-processing-exercises\n", "\n", - "Feel free to contact lecturer jacob.thoenes@uni-rostock.de" + "Feel free to contact the lecturer jacob.thoenes@uni-rostock.de" ] }, { @@ -59,7 +59,7 @@ "source": [ "#### Objectives for Exercises 1 & 2\n", "\n", - "What we should learn in the exercise 1 and 2 is\n", + "What we should learn in exercises 1 and 2 is\n", "- DFT / IDFT as transform\n", "- DFT / IDFT as matrix operation, Fourier matrix as beautiful vector base\n", "- Interpolation of DFT towards DTFT\n", @@ -73,16 +73,16 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "The detailed [DFT Tutorial](dft/dft_windowing_tutorial/dft_windowing_tutorial.pdf) is worth to comprehend. It contains useful information and very detailed calculus on the items above. Furthermore, five exercise tasks including their solutions are given in the tutorial, that might serve as first insight into DFT used for spectral analysis. These exercises are also covered in the dedicated [Exercises of DFT Tutorial](dft/dft_windowing_tutorial/dft_windowing_tutorial_exercises.ipynb) notebook." + "The detailed [DFT Tutorial](dft/dft_windowing_tutorial/dft_windowing_tutorial.pdf) is worth comprehending. It contains useful information and a very detailed calculus on the items above. Furthermore, five exercise tasks, including their solutions, are given in the tutorial, which might serve as a first insight into DFT used for spectral analysis. These exercises are also covered in the dedicated [Exercises of DFT Tutorial](dft/dft_windowing_tutorial/dft_windowing_tutorial_exercises.ipynb) notebook." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "DFT as matrix operation is again covered in [DFT Fundamentals](dft/dft_intro.ipynb), we should not miss this viewpoint! It is also programmed in the [DFT / DTFT of Complex Exponential](dft/dft_complex_exponential.ipynb) notebook.\n", + "DFT as a matrix operation is again covered in [DFT Fundamentals](dft/dft_intro.ipynb); we should not miss this viewpoint! It is also programmed in the [DFT / DTFT of Complex Exponential](dft/dft_complex_exponential.ipynb) notebook.\n", "\n", - "Then, there are more Jupyter notebooks that cover stuff from the [DFT Tutorial](dft/dft_windowing_tutorial/dft_windowing_tutorial.pdf) in a more playground matter with some reference Python code. So if we need a break from all the equations in the pdf file, we might play around with these:\n", + "Then, there are more Jupyter notebooks that cover stuff from the [DFT Tutorial](dft/dft_windowing_tutorial/dft_windowing_tutorial.pdf) in a more playground manner with some reference Python code. So if we need a break from all the equations in the PDF file, we might play around with these:\n", "* [DFT Fundamentals](dft/dft_intro.ipynb)\n", "* [DFT Interpolation to DTFT](dft/dft_to_dtft_interpolation.ipynb)\n", "* [DFT / DTFT of Complex Exponential](dft/dft_complex_exponential.ipynb)\n", @@ -94,16 +94,16 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "After studying we should be able to answer\n", + "After studying, we should be able to answer\n", "\n", "- DFT / IDFT equation and its meaning\n", - "-the matrix versions of it with an idea what the Fourier matrix is and how it is set up\n", + "-the matrix versions of it, with an idea of what the Fourier matrix is and how it is set up\n", "- the viewpoint of i) DFT -> interpolation of spectral data -> DTFT and ii) DTFT -> sampling of spectrum -> DFT\n", - "- the N-DFT of an system with finite impulse response (FIR) of length N contains all spectral data. The frequency response of it is obtained by DFT -> interpolation of spectral data -> DTFT. Note that the reverse process can be used to design an FIR system in the frequency domain, we will study this later in the course.\n", - "- concept of DFT eigenfrequency or in other words: whole signal periods fit into a DFT window\n", - "- what the DFT is actually doing when analyzing signal containing only i) DFT eigenfrequencies or ii) not only DFT eigenfrequencies\n", + "- the N-DFT of a system with a finite impulse response (FIR) of length N contains all spectral data. The frequency response of it is obtained by DFT -> interpolation of spectral data -> DTFT. Note that the reverse process can be used to design an FIR system in the frequency domain. We will study this later in the course.\n", + "- concept of DFT eigenfrequency, or in other words: whole signal periods fit into a DFT window\n", + "- what the DFT is actually doing when analyzing a signal containing only i) DFT eigenfrequencies or ii) not only DFT eigenfrequencies\n", "- the latter brings us to the leakage effect that might be reduced by windowing\n", - "- the compromise between main lobe width and side lobe level of an window and its implications for the leakage effect\n", + "- the compromise between the main lobe width and the side lobe level of a window and its implications for the leakage effect\n", "- windowing in time domain is spectral convolution. The best window would have a Dirac-shaped spectrum, since this behaves neutral in convolution, right? Why is a perfect Dirac spectrum not possible in practice?\n", "- if we want to have a Dirac-like window spectrum, what windows would be suitable, then in practice. Are they good or bad for the leakage effect?\n", "- what are parametric windows, what do we want to parametrize, we should give examples based on the Dolph-Chebychev and the Kaiser Bessel window\n", @@ -114,7 +114,7 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "At this stage it is also worth to think about two fundamental applications heavily used in DSP\n", + "At this stage, it is also worth thinking about two fundamental applications heavily used in DSP\n", "\n", "- FIR filtering: convolution of signal $x$ with a finite impulse response (FIR) $h[k]$\n", "\n", @@ -124,9 +124,9 @@ "\n", "$$y[k] = x[k] \\cdot h[k] \\circ - \\bullet Y(\\Omega) = \\frac{1}{2\\pi} X(\\Omega) \\circledast_{2\\pi} W(\\Omega)$$\n", "\n", - "The discussion on the characteristics of the frequency response of either $H(\\Omega)$ or $W(\\Omega)$ and how these act on the DTFT spectrum of the input signal $x[k]$ is a vital part on DSP. So, when considering finite length signals for $h[k]$ and $w[k]$, we can use the same concepts and principles for designing such signals only using different design goals. A FIR filter application probably needs a different spectrum than an optimized window application.\n", + "The discussion on the characteristics of the frequency response of either $H(\\Omega)$ or $W(\\Omega)$ and how these act on the DTFT spectrum of the input signal $x[k]$ is a vital part of DSP. So, when considering finite length signals for $h[k]$ and $w[k]$, we can use the same concepts and principles for designing such signals, only using different design goals. A FIR filter application probably needs a different spectrum than an optimized window application.\n", "\n", - "So, what role plays the DFT here? We should recall that for finite length signals all spectrum information is stored in the DFT coefficients. The DTFT can then be interpolated from the DFT coefficients, see [DFT Interpolation to DTFT](dft/dft_to_dtft_interpolation.ipynb). So, by knowing the DFT of an FIR filter, we basically know the whole frequency response of this FIR system and we can tell which frequencies in the input signal are gained/attenuated in the convolution process $y[k] = x[k] \\ast h[k]$. And, by knowing the DFT of a window, we know the frequency response of the window and are able to estimate the smoothing effect that happens in the convolution $Y(\\Omega) = \\frac{1}{2\\pi} X(\\Omega) \\circledast_{2\\pi} W(\\Omega)$. That's why we should learn about DFT in very detail. We could check \n", + "So, what role does the DFT play here? We should recall that for finite length signals all spectrum information is stored in the DFT coefficients. The DTFT can then be interpolated from the DFT coefficients, see [DFT Interpolation to DTFT](dft/dft_to_dtft_interpolation.ipynb). So, by knowing the DFT of an FIR filter, we basically know the whole frequency response of this FIR system, and we can tell which frequencies in the input signal are gained/attenuated in the convolution process $y[k] = x[k] \\ast h[k]$. And, by knowing the DFT of a window, we know the frequency response of the window and are able to estimate the smoothing effect that happens in the convolution $Y(\\Omega) = \\frac{1}{2\\pi} X(\\Omega) \\circledast_{2\\pi} W(\\Omega)$. That's why we should learn about DFT in very detail. We could check \n", "[Pole/Zero Plot and Frequency Response of Windows](dft/window_zplane_frequency_response.ipynb)." ] }, @@ -143,13 +143,13 @@ "\n", "#### Objectives for Exercise 3\n", "\n", - "We should get an idea on\n", + "We should get an idea of\n", "- the probability density function (PDF)\n", "- what is a sample function of a random process\n", - "- first / second order ensemble averages (moments)\n", + "- first/second order ensemble averages (moments)\n", "- the concept of stationarity and ergodicity\n", "- the concept of temporal average vs. ensemble average\n", - "- the auto correlation / cross correlation as as a higher order ensemble and temporal averages" + "- the auto correlation/cross correlation as a higher order ensemble and temporal averages" ] }, { @@ -171,10 +171,10 @@ "We should get familiar with the following aspects\n", "- the correlation function and its DTFT spectrum, which is the so called power spectral density (PSD) \n", "- the concept of cross- vs. auto-correlation function, we make sure that we correctly understood the convention for the shifting direction\n", - "- the four fundamental **equations I-IV** how the correlation functions and their power spectral densities are connected for in-/output signal of an LTI system (we should not miss these as they are really **important fundamentals**)\n", + "- the four fundamental **equations I-IV**, how the correlation functions and their power spectral densities are connected for in-/output signal of an LTI system (we should not miss these as they are really **important fundamentals**)\n", "- correlation of sine/cosine signals to get the idea of maximum positive/negative correlation and no correlation at all\n", "- the auto-correlation function (ACF) of a simple short sequence in very detail\n", - "- analytic calculus of simple lowpass filters FIR and IIR (practical examples are much more complicated, such that analytic treatment is no fun anymore...also we here can make sure that we got the time mirroring concept right)\n", + "- analytic calculus of simple lowpass filters FIR and IIR (practical examples are much more complicated, such that analytic treatment is no fun anymore...also we can make sure that we got the time mirroring concept right)\n", "- simple method of estimating an LTI system's impulse response via cross correlation" ] }, @@ -184,7 +184,7 @@ "source": [ "### 5. Spectral Estimation of Random Signals\n", "\n", - "#### Material is part of homework assignment\n", + "#### Material is part of a homework assignment\n", "\n", "* [STFT](spectral_estimation/stft.ipynb)\n", "* Welch method" @@ -202,18 +202,18 @@ "\n", "#### Objectives for Exercise 5\n", "\n", - "In this exercise we will have a detailed look to the concept of amplitude quantization.\n", + "In this exercise, we will have a detailed look at the concept of amplitude quantization.\n", "\n", "We should get familiar with\n", "- the quantization error model and its underlying assumptions\n", - "- why some of the assumptions are violated; for example when quantizing a sine signal\n", + "- why some of the assumptions are violated; for example, when quantizing a sine signal\n", "- the uniform (saturated) midtread quantization curve (characteristics, shape, quantization step size)\n", - "- the link between number of quantization steps Q and number of encoding bits\n", + "- the link between the number of quantization steps Q and the number of encoding bits\n", "- the signal-to-'quantization noise'-ratio for different input signals and the resulting rule of thumb \"6dB better SNR per bit\"\n", "- the concept of dithering (why, how)\n", - "- the difference of using dithering noise with either a uniform PDF or a triangular PDF and why is triangular better\n", + "- the difference of using dithering noise with either a uniform PDF or a triangular PDF, and why is triangular better\n", "\n", - "We should listen to the audio examples for the dithering exercise. Especially the case when the sine amplitude is below the quantization step is interesting, as one would assume that this sine signal cannot be perceived at all. **Please, DO NOT HARM YOUR EARS! Start with low playback level.**" + "We should listen to the audio examples for the dithering exercise. Especially the case when the sine amplitude is below the quantization step is interesting, as one would assume that this sine signal cannot be perceived at all. **Please, DO NOT HARM YOUR EARS! Start with a low playback level.**" ] }, { @@ -228,7 +228,7 @@ "\n", "#### Objectives for Exercise 6\n", "\n", - "In this exercise we deal with FIR filters, more precisely with non-recursive LTI systems.\n", + "In this exercise, we deal with FIR filters, more precisely with non-recursive LTI systems.\n", "\n", "We should get familiar with\n", "- filter fundamentals (recursive vs. non-recursive systems, signal flow diagram, transfer function)\n", @@ -251,12 +251,12 @@ "\n", "#### Objectives for Exercise 7\n", "\n", - "In this exercise we deal with recursive systems, more precise with recursive systems that exhibit an infinite impulse response, which we call IIR filter (most textbooks do this as well).\n", + "In this exercise, we deal with recursive systems, more precisely, with recursive systems that exhibit an infinite impulse response, which we call an IIR filter (most textbooks do this as well).\n", "\n", "We should get familiar with\n", "- filter fundamentals, especially the fundamental difference between recursive and non-recursive systems in terms of their transfer function, poles, stability, signal flow graph\n", - "- simple examples what a zero is doing, what a pole is doing (practical filters are just more complicated, but the key idea is the same) \n", - "- bilinear transform as an exemplary straightforward, simple IIR filter design method when Laplace transfer function is given\n", + "- simple examples of what a zero is doing, what a pole is doing (practical filters are just more complicated, but the key idea is the same) \n", + "- bilinear transform as an exemplary, straightforward, simple IIR filter design method when the Laplace transfer function is given\n", "- concept of frequency pre-warping and bandwidth / Q-factor pre-warping when applying the bilinear transform" ] }, @@ -324,7 +324,7 @@ "- Klaus Göldner\n", "- Udo Zölzer\n", "\n", - "If one looks for more practical approaches, the books in German language by\n", + "If one looks for more practical approaches, the books in the German language by\n", "\n", "- Daniel Ch. von Grünigen\n", "- Martin Werner\n", @@ -339,7 +339,7 @@ "\n", "- Alan V. Oppenheim, MIT, https://ocw.mit.edu/resources/res-6-008-digital-signal-processing-spring-2011/\n", "- edX: https://www.edx.org/search?q=digital%20signal%20processing \n", - "- Julius O. Smith III, Standford University, Online Books and Courses at https://ccrma.stanford.edu/~jos/\n", + "- Julius O. Smith III, Stanford University, Online Books and Courses at https://ccrma.stanford.edu/~jos/\n", "- Gert Leus, TU Delft, https://ocw.tudelft.nl/courses/digital-signal-processing/\n", "\n", "Finally, although not traditionally related to DSP, it is recommended to have a look at\n", @@ -356,7 +356,7 @@ "## Contributions\n", "\n", "- Frank Schultz (concept, initial main author)\n", - "- Jacob Thönes (code)\n", + "- Jacob Thönes (code, proof read)\n", "- Robert Hauser (code for graphics generation on DFT material)\n", "- Vera Erbes (author, proof read, translation)\n", "- Sascha Spors (concept, proof read)" @@ -388,7 +388,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.10.12" + "version": "3.12.2" } }, "nbformat": 4, diff --git a/jupyter/jupyter_intro.ipynb b/jupyter/jupyter_intro.ipynb index 78864c5..fff4be4 100644 --- a/jupyter/jupyter_intro.ipynb +++ b/jupyter/jupyter_intro.ipynb @@ -19,7 +19,7 @@ "- lecture: https://github.com/spatialaudio/digital-signal-processing-lecture\n", "- tutorial: https://github.com/spatialaudio/digital-signal-processing-exercises\n", "\n", - "Feel free to contact lecturer jacob.thoenes@uni-rostock.de" + "Feel free to contact the lecturer jacob.thoenes@uni-rostock.de" ] }, { @@ -41,10 +41,10 @@ "source": [ "**Anaconda Environment**\n", "\n", - "The [Anaconda distribution](https://www.anaconda.com/distribution/) is a convenient solution to install a required environment, i.e. to have access to the Jupyter Notebook renderer with a Python interpreter on a personal computer. It is very likely that a very recent installation of Anaconda already delivers all required packages just using the `base` environment. It is however good practice to create a dedicated environment for each project. So, for this tutorial we might use a `mydsp` (or whatever name works for us) environment.\n", + "The [Anaconda distribution](https://www.anaconda.com/distribution/) is a convenient solution to install a required environment, i.e., to have access to the Jupyter Notebook renderer with a Python interpreter on a personal computer. It is very likely that a very recent installation of Anaconda already delivers all required packages just using the `base` environment. It is, however, good practice to create a dedicated environment for each project. So, for this tutorial,l we might use a `mydsp` (or whatever name works for us) environment.\n", "\n", - "- get into the folder where the exercises are located, e.g. `cd my_dsp_folder`\n", - "- in the subfolder `.binder` the `environment.yml` can be used to create a dedicated conda `mydsp` environment as\n", + "- get into the folder where the exercises are located, e.g., `cd my_dsp_folder`\n", + "- in the subfolder `.binder`, the `environment.yml` can be used to create a dedicated conda `mydsp` environment as\n", " - `conda env create -f environment.yml --force`\n", " - we can remove this environment with `conda env remove --name mydsp`\n", "- this should also have installed audio related libraries using pip\n", @@ -72,7 +72,6 @@ "# such as fft, signal\n", "import matplotlib as mpl\n", "import numpy as np\n", - "\n", "from matplotlib import pyplot as plt\n", "\n", "# from numpy import fft # use either numpy fft or\n", @@ -89,8 +88,8 @@ "The `np.matrix` package is (was) meant for linear algebra on matrices, which by\n", "definition are of dimension m rows $\\times$ n columns, i.e. 2D.\n", "So this might be what you're looking for when dealing with linear algebra.\n", - "**However**, the community does not recommend to use this package anymore and\n", - "it might be even removed in future." + "**However**, the community does not recommend using this package anymore, and\n", + "it might even be removed in the future." ] }, { @@ -154,11 +153,11 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "Note here, that the shape is not (4,1) as above but rather (4,), this array has one axis. Note that transpose is doing nothing for arrays of rank 1. So this representation differs from the concepts of a column vector $(m \\times 1)$ and a row vector $(1 \\times n)$, which by means of linear algebra and the matrix idea have two axis, since they are very special matrices.\n", + "Note here that the shape is not (4,1) as above, but rather (4,); this array has one axis. Note that transpose is doing nothing for arrays of rank 1. So this representation differs from the concepts of a column vector $(m \\times 1)$ and a row vector $(1 \\times n)$, which by means of linear algebra and the matrix idea have two axes, since they are very special matrices.\n", "\n", - "So, if we really need a *classic* column or row vector, we need to tell this numpy explicitly.\n", + "So, if we really need a *classic* column or row vector, we need to tell numpy explicitly.\n", "We have done this with the command `np.newaxis` in the above example.\n", - "Especially, having strong experience with Matlab and starting with the concept of rank-1 vs. rank-2 arrays might produce initial headache. However, the concept is very powerful, so we should get used to it.\n", + "Especially, having strong experience with Matlab and starting with the concept of rank-1 vs. rank-2 arrays might produce an initial headache. However, the concept is very powerful, so we should get used to it.\n", "\n", "Check https://numpy.org/devdocs/user/numpy-for-matlab-users.html for Matlab vs. numpy\n", "\n", @@ -304,7 +303,7 @@ "# Inner / Outer Product\n", "\n", "Once we get used to rank-1 arrays, we can do nice coding on (what we actually interpret as) vectors (and later\n", - "of course on matrices and tensors).\n", + "of course, on matrices and tensors).\n", "\n", "Consider the following matrix A" ] @@ -385,7 +384,7 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "The elegant part of the story is, that we directly see what the code shall do with the data. In other words, we know that the result's dimension is actually a scalar (`inner`) or a matrix (`outer`). We don't need to check the actual data to tell this. Just, compare this with typical Matlab code lines `a.' * b` or `a * b.'`. Here, we better should know the data dimensions to get an idea what the code is intended to do. This is only one advantage of numpy's array concept. " + "The elegant part of the story is that we directly see what the code shall do with the data. In other words, we know that the result's dimension is actually a scalar (`inner`) or a matrix (`outer`). We don't need to check the actual data to tell this. Just compare this with typical Matlab code lines `a.' * b` or `a * b.'`. Here, we better should know the data dimensions to get an idea of what the code is intended to do. This is only one advantage of numpy's array concept. " ] }, { @@ -594,9 +593,9 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "These are complex exponentials and due to the chosen parameters, periodic in $N$. In fact, these can be considered as DFT eigensignals for $N=32$, actually for $\\mu=1$ and $\\mu=2$.\n", + "These are complex exponentials and, due to the chosen parameters, periodic in $N$. In fact, these can be considered as DFT eigensignals for $N=32$, actually for $\\mu=1$ and $\\mu=2$.\n", "\n", - "Let's just plot these signals into one graph. The fastest way is to use `pyplot`, which is an API for `matplotlib` and similar to Matlab. `pyplot` is a good tool for quick'n dirty plots, whereas `matplotlib` gives us full access to any plotting objects, that's the professional way to plot. We leave it here with few simple calls of the `pyplot` API. " + "Let's just plot these signals into one graph. The fastest way is to use `pyplot`, which is an API for `matplotlib` and similar to Matlab. `pyplot` is a good tool for quick'n dirty plots, whereas `matplotlib` gives us full access to any plotting objects, that's the professional way to plot. We leave it here with a few simple calls of the `pyplot` API. " ] }, { @@ -688,7 +687,7 @@ "source": [ "We get expected results for orthonormal vectors, besides numerical precision errors.\n", "\n", - "If you don't like complex signals / complex vector space that much, check it\n", + "If you don't like complex signals / complex vector spaces that much, check it\n", "with plain unit amplitude cos() and sin() signals, periodic in $N$, where full\n", "periods fit into the signal length (= vector dimension)." ] @@ -763,7 +762,7 @@ "\n", "This is a simple example of a surface plot using\n", "- `pcolormesh` called with the `matplotlib` API\n", - "- discrete valued colorbar based on `viridis` colormap. You might also check the `plasma`, `inferno`,`magma`,`cividis` colormaps for perceptually uniform sequential colormaps. If you need a diverging colormap (such as for nicely indicating positive and negative amplitudes of a waveform) `RdBu`, `seismic`, `bwr` (for non-red/blue colorblind people) do a good job. Colormaps like `jet` or `hsv` are not recommended, they do not match very well with our visual perception of these colorspaces." + "- discrete valued colorbar based on `viridis` colormap. You might also check the `plasma`, `inferno`,`magma`,`cividis` colormaps for perceptually uniform sequential colormaps. If you need a diverging colormap (such as for nicely indicating positive and negative amplitudes of a waveform) `RdBu`, `seismic`, `bwr` (for non-red/blue colorblind people) do a good job. Colormaps like `jet` or `hsv` are not recommended, as they do not match very well with our visual perception of these colorspaces." ] }, { @@ -880,7 +879,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.10.12" + "version": "3.12.2" } }, "nbformat": 4, diff --git a/quantization/quantization.ipynb b/quantization/quantization.ipynb index f698dd2..b5db0bc 100644 --- a/quantization/quantization.ipynb +++ b/quantization/quantization.ipynb @@ -19,7 +19,7 @@ "- lecture: https://github.com/spatialaudio/digital-signal-processing-lecture\n", "- tutorial: https://github.com/spatialaudio/digital-signal-processing-exercises\n", "\n", - "Feel free to contact lecturer jacob.thoenes@uni-rostock.de" + "Feel free to contact the lecturer jacob.thoenes@uni-rostock.de" ] }, { @@ -50,13 +50,13 @@ "outputs": [], "source": [ "# most common used packages for DSP, have a look into other scipy submodules\n", - "import numpy as np\n", "import matplotlib as mpl\n", "import matplotlib.pyplot as plt\n", - "from scipy import signal\n", + "import numpy as np\n", "\n", "# audio write and play stuff\n", "import soundfile as sf # requires 'pip install soundfile'\n", + "from scipy import signal\n", "\n", "# last tested with soundfile-0.10.3" ] @@ -168,16 +168,16 @@ "\n", "Quantization generates signals that have discrete values $x_q[k]$, $x_q(t)$ from signals with continuous values $x[k]$, $x(t)$. \n", "\n", - "For quantization, the signals can be both, discrete and continuous in time.\n", + "For quantization, the signals can be both discrete and continuous in time.\n", "However, a signal that is discrete in time **and** discrete in value is termed a **digital** signal.\n", "Only digital signals can be processed by computers.\n", - "Here the quantization of discrete-time signals is treated due to practical importance.\n", + "Here, the quantization of discrete-time signals is treated due to practical importance.\n", "\n", "To describe quantization analytically, the model in the figure below is used.\n", "\n", "\"QuantizationModel.png\"\n", "\n", - "The input and output signal differ by the so called quantization error (quantization noise) $e[k]$, that is defined as\n", + "The input and output signal differ by the so called quantization error (quantization noise) $e[k]$, which is defined as\n", "\n", "\\begin{equation}\n", "e[k] = x_q[k] - x[k],\n", @@ -194,7 +194,7 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "To use this error model, some assumption have to be made.\n", + "To use this error model, some assumptions have to be made.\n", "The quantization noise shall be uniformly distributed, which then can be modeled with the probability density function (PDF) $p_e(\\theta) = \\frac{1}{\\Delta Q} \\mathrm{rect}(\\frac{\\theta_e}{\\Delta Q})$, where $\\Delta Q$ denotes the quantization step size and $\\theta_e$ the amplitudes of the quantization error signal.\n", "This PDF is shown below." ] @@ -289,11 +289,11 @@ "\n", "However, this is not observable in this example!\n", "\n", - "Instead, from the above plot, we can deduce that $e[k]$ is correlated to itself, i.e. it exhibits periodicity each 100 samples in phase, and each 50 sample out of phase.\n", - "The sine period is precisely 100 samples, thus the input signal and the quantization error are somewhat linked and not independent.\n", + "Instead, from the above plot, we can deduce that $e[k]$ is correlated to itself, i.e., it exhibits periodicity every 100 samples in phase, and every 50 samples out of phase.\n", + "The sine period is precisely 100 samples; thus, the input signal and the quantization error are somewhat linked and not independent.\n", "Thus, the error model assumption is violated. That is bad, since the sine signal allows for otherwise comparable simple analytical calculus.\n", "\n", - "The links between the signals can be furthermore confirmed with the help of the cross-correlation functions.\n", + "The links between the signals can be further confirmed with the help of the cross-correlation functions.\n", "Their oscillating characteristics reveal that quantization error is highly correlated. " ] }, @@ -328,14 +328,14 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "Therefore, the special case of sine signals is in fact not suited for the quantization model above.\n", - "Because of the simplicity of the involved calculation it is common practice to conduct this analysis for sine signals nevertheless, and signal-to-noise ratios in the data sheets of A/D converters are mostly stated for excitation with sine signals.\n", - "For random signals, the quantization model is only valid for high levels in the quantizer. For more information see\n", + "Therefore, the special case of sine signals is, in fact, not suited for the quantization model above.\n", + "Because of the simplicity of the involved calculation, it is common practice to conduct this analysis for sine signals; nevertheless, signal-to-noise ratios in the data sheets of A/D converters are mostly stated for excitation with sine signals.\n", + "For random signals, the quantization model is only valid for high levels in the quantizer. For more information, see\n", "[Udo Zölzer, Digital Audio Signal Processing, Wiley](https://onlinelibrary.wiley.com/doi/book/10.1002/9780470680018)\n", "(might be available as free access in your uni network)\n", "\n", "- Task:\n", - "Increase the (odd) number of quantization steps $Q$ and check what happens with the shape and amplitudes of the correlations functions. Hint: closer look to the amplitudes of the correlation signals." + "Increase the (odd) number of quantization steps $Q$ and check what happens with the shape and amplitudes of the correlation functions. Hint: a closer look at the amplitudes of the correlation signals." ] }, { @@ -375,10 +375,10 @@ "\n", "where for $\\lfloor \\cdot \\rfloor$ a rounding operation might used and $\\backslash$ denotes integer division.\n", "So the mapping functions $g$ and $f$ are simple multiplications.\n", - "At the beginning of this notebook, the function `my_quant` is implemented that realizes quantization based on this mapping.\n", + "At the beginning of this notebook, the function `my_quant` is implemented, which realizes quantization based on this mapping.\n", "The approach uses `numpy`'s `round` operation.\n", "When asking for rounding, care has to be taken, which [approach](https://en.wikipedia.org/wiki/Rounding) shall be used.\n", - "Numpy rounds to the nearest **even** integer in contrast to e.g. Matlab's rounding to nearest integer.\n", + "Numpy rounds to the nearest **even** integer in contrast to e.g., Matlab's rounding to nearest integer.\n", "\n", "Detailed analysis for `my_quant`:\n", "\n", @@ -396,7 +396,7 @@ "\n", "The case of **even** $Q$ is practically used for virtually all analog/digital (AD) and digital/analog (DA) converters.\n", "\n", - "When additionally to the above statements\n", + "When addition to the above statements\n", "\n", "\\begin{equation}\n", "\\log_2(Q)\\in\\mathbb{N}\n", @@ -468,17 +468,17 @@ "source": [ "### Plotting the Midtread Curve\n", "\n", - "We now can visualize the characteristic curve for a simple, made up input signal, i.e. a monotonic increasing signal between $x_{max} = -x_{min}$ using an equidistant increment $\\Delta x$ over sample index $k$.\n", + "We can now visualize the characteristic curve for a simple, made-up input signal, i.e., a monotonic increasing signal between $x_{max} = -x_{min}$ using an equidistant increment $\\Delta x$ over sample index $k$.\n", "\n", "Here, we use $x_{max} = 1.25$ and $\\Delta x=0.001$ and assume that we start with $x_{min} = -1.25$ at $k=0$.\n", - "If $\\Delta x$ is sufficiently small, the signal's amplitude can be interpreted as continuous straight line.\n", + "If $\\Delta x$ is sufficiently small, the signal's amplitude can be interpreted as a continuous straight line.\n", "This straight line is degraded in a quantization process.\n", - "Plotting the quantization result over the input, results in the characteristic curve, in our example in the curve of the uniform saturated midtread quantizer.\n", + "Plotting the quantization result over the input results in the characteristic curve, in our example, the curve of the uniform saturated midtread quantizer.\n", "\n", "Let us plot this.\n", "\n", "**Please note:**\n", - "The quantizer `uniform_midtread_quantizer` known from lecture and `my_quant` yield the same results besides a slight detail: `uniform_midtread_quantizer` always exhibits an **even** number of quantization steps $Q$.\n", + "The quantizer `uniform_midtread_quantizer` known from the lecture and `my_quant` yield the same results, besides a slight detail: `uniform_midtread_quantizer` always exhibits an **even** number of quantization steps $Q$.\n", "So, only for even $Q$ results are exactly identical.\n", "\n", "We might verify this in the next plots as well." @@ -532,14 +532,14 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "The following exercises used to be a homework assignment as exam's prerequisite.\n", + "The following exercises used to be a homework assignment as an exam's prerequisite.\n", "\n", "# Exercise 1: Uniform Saturated Midtread Characteristic Curve of Quantization\n", "\n", "## Task\n", "\n", "Check this quantizer curve for $Q=7$ and $Q=8$.\n", - "Make sure that you get the idea of the midtread concept (the zero is always quantized to zero) and saturation (for even $Q$) largest quantization step is saturated)." + "Make sure that you get the idea of the midtread concept (the zero is always quantized to zero) and saturation (for even $Q$) the largest quantization step is saturated)." ] }, { @@ -607,18 +607,18 @@ "source": [ "# Exercise 2: Quantization and Signal-to-Noise Ratio\n", "\n", - "From theory the **6dB / Bit rule of thumb** is well known for uniform quantization. It states that the signal-to-noise ratio increases by 6 dB for every additional bit that is spent to quantize the input data.\n", + "From theory, the **6dB / Bit rule of thumb** is well known for uniform quantization. It states that the signal-to-noise ratio increases by 6 dB for every additional bit that is spent to quantize the input data.\n", "Hence,\n", "\n", "\\begin{equation}\n", "\\text{SNR in dB} = 6 \\cdot B + \\gamma,\n", "\\end{equation}\n", "\n", - "where $\\gamma$ is a offset value in dB that depends on the PDF of the signal to be quantized.\n", - "Note that this rule of thumb assumes that the quantization error exhibits uniform PDF and is not correlated with the quantized signal.\n", - "We can see that this rule of thumb is inaccurate when quantizing a sine signal with small number of bits or an amplitude in the range of the quantization step. Then, the mentioned assumptions are not fulfilled. We will observe this in exercise 3.\n", + "where $\\gamma$ is an offset value in dB that depends on the PDF of the signal to be quantized.\n", + "Note that this rule of thumb assumes that the quantization error exhibits a uniform PDF and is not correlated with the quantized signal.\n", + "We can see that this rule of thumb is inaccurate when quantizing a sine signal with a small number of bits or an amplitude in the range of the quantization step. Then, the mentioned assumptions are not fulfilled. We will observe this in exercise 3.\n", "\n", - "We plot the function SNR of bits below for uniform, normal and Laplace PDF noises and a sine signal.\n", + "We plot the function SNR of bits below for uniform, normal, and Laplace PDF noises and a sine signal.\n", "We should observe the slope of always 6 dB per Bit.\n", "We should note the different absolute values of the SNR depending on the varying $\\gamma$. \n", "\n", @@ -785,7 +785,7 @@ "\n", "This dither signal with small amplitudes aims at de-correlating the quantization error $e[k]$ from the quantized signal $x_q[k]$, which is especially important for small amplitudes of $x[k]$.\n", "This technique is called **dithering**.\n", - "For $d[k]=0$ no dithering is applied.\n", + "For $d[k]=0$, no dithering is applied.\n", "\n", "Since the quantization error may be in the range $-\\frac{\\Delta Q}{2}\\leq e[k]\\leq \\frac{\\Delta Q}{2}$ (assuming uniform distribution), it appears reasonable to use a dither noise with a probability density function (PDF) of\n", "\n", @@ -795,7 +795,7 @@ "\n", "i.e. a **zero-mean, uniformly distributed noise** with maximum amplitude $|d[k]|=\\frac{\\Delta Q}{2}$.\n", "It can be shown that this dither noise improves the quality of the quantized signal.\n", - "However, there is still a noise modulation (i.e. a too high correlation between $x_q[k]$ and $e[k]$) that depends on the amplitude of the input signal.\n", + "However, there is still a noise modulation (i.e., a too high correlation between $x_q[k]$ and $e[k]$) that depends on the amplitude of the input signal.\n", "\n", "The noise modulation can be almost completely eliminated with a **zero-mean noise** signal exhibiting a **symmetric triangular PDF**\n", "\n", @@ -805,7 +805,7 @@ "\n", "with maximum amplitude $|d[k]|=Q$.\n", "By doing so, an almost ideal decorrelation between $x_q[k]$ and $e[k]$ is realized.\n", - "In audio, this technique is called TPDF-Dithering (Triangular Probability Density Function Dithering) and can be applied in the mastering process of audio material that is to be distributed e.g. on a CD or via streaming." + "In audio, this technique is called TPDF-Dithering (Triangular Probability Density Function Dithering) and can be applied in the mastering process of audio material that is to be distributed, e.g., on a CD or via streaming." ] }, { @@ -814,7 +814,7 @@ "source": [ "## Task \n", "\n", - "To get an impression on how dithering may be implemented and what quantized signals sound like, the following exercises shall be performed.\n", + "To get an impression of how dithering may be implemented and what quantized signals sound like, the following exercises shall be performed.\n", "\n", "- Generate the sine signal $x[k]$ defined above.\n", "\n", @@ -822,23 +822,23 @@ "\n", "- Generate the dither noise $d_\\text{TRI}[k]$ according to the PDF $p_\\text{TRI}(d)$. Check the resulting amplitude and distribution carefully. The length of $d_\\text{TRI}[k]$ and $x[k]$ must be equal.\n", "\n", - "- Add each dither noise $d_\\text{RECT}[k]$ and $d_\\text{TRI}[k]$ individually to $x[k]$. Together with the case of no dithering we now have three signals to be quantized.\n", + "- Add each dither noise $d_\\text{RECT}[k]$ and $d_\\text{TRI}[k]$ individually to $x[k]$. Together with the case of no dithering, we now have three signals to be quantized.\n", "\n", "- Quantize these signals individually using `my_quant(x,Q)` with $Q$ quantization steps.\n", "\n", "- Plot the midtread characteristic curve.\n", "\n", - "- Plot the histogram of the dither noises as estimate of its PDF.\n", + "- Plot the histogram of the dither noises as an estimate of its PDF.\n", "\n", - "- Plot the histogram of the error noises as estimate of its PDF.\n", + "- Plot the histogram of the error noises as an estimate of its PDF.\n", "\n", - "- Plot the sine signal, the dithered signal, the quantized signal and the quantization error signal in one diagram for all three cases.\n", + "- Plot the sine signal, the dithered signal, the quantized signal, and the quantization error signal in one diagram for all three cases.\n", "\n", "- Calculate and plot the CCF of the signals $x_q[k]$ and $e[k]$ for all three cases.\n", "\n", "- Interpret the obtained graphics.\n", "\n", - "- For each case, render WAV files from $x[k]$, $x[k]+d[k]$, $x_q[k]$ und $e[k]$ and listen to them. **Be careful! Do not harm your ears!** Pay special attention to the sound of the quantization error, how it is correlated with the quantized signal and how loud it appears.\n", + "- For each case, render WAV files from $x[k]$, $x[k]+d[k]$, $x_q[k]$ und $e[k]$ and listen to them. **Be careful! Do not harm your ears!** Pay special attention to the sound of the quantization error, how it is correlated with the quantized signal, and how loud it appears.\n", "\n", "- Consider the 5 cases\n", "\n", @@ -848,12 +848,12 @@ " 4. $B=3$ Bit, $A=\\Delta Q$\n", " 5. $B=3$ Bit, $A=\\frac{\\Delta Q}{2}$\n", "\n", - "In the last case the signal has amplitude even below the quantization step size $\\Delta Q$. You might verify by listening that the sine is still perceivable if dithering is applied, but not if no dithering is applied.\n", + "In the last case, the signal has amplitude even below the quantization step size $\\Delta Q$. You might verify by listening that the sine is still perceivable if dithering is applied, but not if no dithering is applied.\n", "\n", "**Again: Be careful! Do not harm your ears!**\n", - "The signal amplitude $A$ and chosen $B$ is directly related to the playback level!\n", + "The signal amplitude $A$ and the chosen $B$ are directly related to the playback level!\n", "\n", - "**Warning again: start with very very low playback level, find the loudest signal first and then increase volume to your convenience**" + "**Warning again: start with very, very low playback level, find the loudest signal first, and then increase volume to your convenience**" ] }, { @@ -1174,7 +1174,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.10.12" + "version": "3.12.2" } }, "nbformat": 4, diff --git a/random_signals/correlation.ipynb b/random_signals/correlation.ipynb index 0db4806..d3c1d02 100644 --- a/random_signals/correlation.ipynb +++ b/random_signals/correlation.ipynb @@ -18,7 +18,7 @@ "- lecture: https://github.com/spatialaudio/digital-signal-processing-lecture\n", "- tutorial: https://github.com/spatialaudio/digital-signal-processing-exercises\n", "\n", - "Feel free to contact lecturer jacob.thoenes@uni-rostock.de\n", + "Feel free to contact the lecturer jacob.thoenes@uni-rostock.de\n", "\n", "WIP..." ] @@ -88,9 +88,9 @@ "source": [ "## Normalization schemes for cross correlation of finite length signals\n", "\n", - "check cross correlation\n", + "Check cross correlation\n", "- of a cosine and a sine signal\n", - "- of a normal pdf process that exhibits some repetition" + "- of a normal PDF process that exhibits some repetition" ] }, { @@ -176,7 +176,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.10.12" + "version": "3.12.2" } }, "nbformat": 4, diff --git a/random_signals/ensemble_averages.ipynb b/random_signals/ensemble_averages.ipynb index d8001fd..3b3ada5 100644 --- a/random_signals/ensemble_averages.ipynb +++ b/random_signals/ensemble_averages.ipynb @@ -19,7 +19,7 @@ "- lecture: https://github.com/spatialaudio/digital-signal-processing-lecture\n", "- tutorial: https://github.com/spatialaudio/digital-signal-processing-exercises\n", "\n", - "Feel free to contact lecturer jacob.thoenes@uni-rostock.de\n", + "Feel free to contact the lecturer jacob.thoenes@uni-rostock.de\n", "\n", "WIP..." ] @@ -136,7 +136,7 @@ "which reads quadratic mean is linear mean plus variance.\n", "\n", "\n", - "For **stationary processes** these ensemble averages are not longer time-dependent, but rather $\\mu_x[k] = \\mu_x = \\mathrm{const}$, etc. holds.\n", + "For **stationary processes**, these ensemble averages are no longer time-dependent, but rather $\\mu_x[k] = \\mu_x = \\mathrm{const}$, etc. holds.\n", "This implies that the PDF describing the random process is **not changing over time**." ] }, @@ -163,7 +163,7 @@ "p_{xy}(\\theta_x, \\theta_y, k_x, k_y) = p_{xy}(\\theta_x, \\theta_y, \\kappa).\n", "\\end{equation}\n", "\n", - "For **stationary processes** two important cases lead to fundamental tools for random signal processing:\n", + "For **stationary processes**, two important cases lead to fundamental tools for random signal processing:\n", "\n", "- Case 1: $\\kappa = 0$, i.e. $k = k_x = k_y$\n", "- Case 2: $\\kappa \\neq 0$" @@ -238,7 +238,7 @@ "\n", "## Wide-Sense Ergodic\n", "\n", - "ergodicity holds for linear mapping\n", + "Ergodicity holds for a linear mapping\n", "\n", "\\begin{equation}\n", "\\overline{ x_n[k] \\cdot x_n[k-\\kappa] } = E\\{ x[k] \\cdot x[k-\\kappa] \\} \\;\\; \\forall n\n", @@ -278,7 +278,7 @@ "\\lim_{K \\to \\infty} \\frac{1}{2K + 1} \\sum_{k=-K}^{K} x[k] \\cdot x[k-\\kappa].\n", "\\end{equation}\n", "\n", - "These equations hold for power signals, i.e. the summation yields a finite value.\n" + "These equations hold for power signals, i.e., the summation yields a finite value.\n" ] }, { @@ -372,11 +372,11 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "# Example: Histogram of Gaussian Noise, Cosine and Rectangular Signal\n", + "# Example: Histogram of Gaussian Noise, Cosine, and Rectangular Signal\n", "\n", - "Here we use the numpy histogram with fixed number of bins and really the histogram mode rather than the density mode.\n", + "Here, we use the numpy histogram with a fixed number of bins and really the histogram mode rather than the density mode.\n", "\n", - "We here do not strictly deal with random sample functions (for the cosine and rect), but with amplitude values over time. We do this for practical purpose however, since it is nice to get an idea what a histogram looks like for known signals." + "We here do not strictly deal with random sample functions (for the cosine and rect), but with amplitude values over time. We do this for practical purposes, however, since it is nice to get an idea of what a histogram looks like for known signals." ] }, { @@ -725,7 +725,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.10.12" + "version": "3.12.2" } }, "nbformat": 4, diff --git a/random_signals_LTI_systems/random_LTI_time_domain.ipynb b/random_signals_LTI_systems/random_LTI_time_domain.ipynb index 1d0c278..0330ad7 100644 --- a/random_signals_LTI_systems/random_LTI_time_domain.ipynb +++ b/random_signals_LTI_systems/random_LTI_time_domain.ipynb @@ -19,7 +19,7 @@ "- lecture: https://github.com/spatialaudio/digital-signal-processing-lecture\n", "- tutorial: https://github.com/spatialaudio/digital-signal-processing-exercises\n", "\n", - "Feel free to contact lecturer jacob.thoenes@uni-rostock.de" + "Feel free to contact the lecturer jacob.thoenes@uni-rostock.de" ] }, { @@ -43,14 +43,14 @@ "\\varphi_{xy}[\\kappa]=\\frac{1}{N}\\sum\\limits_{k=0}^{N-1} x[k+\\kappa]\\cdot y^*[k] = \\frac{1}{N}\\sum\\limits_{k=0}^{N-1} x[k]\\cdot y^*[k-\\kappa],\n", "\\end{equation}\n", "\n", - "following the notation convention of the lecture and denoting $^*$ as complex-conjugate.\n", + "following the notation convention of the lecture and denoting $^*$ as the complex-conjugate.\n", "A very important property is\n", "\n", "\\begin{equation}\n", "\\varphi_{xy}[\\kappa]=\\varphi^*_{yx}[-\\kappa].\n", "\\end{equation}\n", "\n", - "Thus, only if the CCF is real-valued - which is the case if $x[k]$ and $y[k]$ are real valued - the equality $\\varphi_{xy}[\\kappa]=\\varphi_{yx}[-\\kappa]$ holds. Otherwise $\\varphi_{xy}[\\kappa]$ and $\\varphi_{yx}[\\kappa]$ are linked by time-mirroring **and** complex-conjugate." + "Thus, only if the CCF is real-valued, which is the case if $x[k]$ and $y[k]$ are real valued, the equality $\\varphi_{xy}[\\kappa]=\\varphi_{yx}[-\\kappa]$ holds. Otherwise $\\varphi_{xy}[\\kappa]$ and $\\varphi_{yx}[\\kappa]$ are linked by time-mirroring **and** complex-conjugate." ] }, { @@ -106,14 +106,14 @@ "\\varphi_{xx}[\\kappa]\\leq\\varphi_{xx}[\\kappa=0]=\\text{E}\\left\\{x[k]\\cdot x[k]\\right\\}=\\sigma_x^2+\\mu_x^2.\n", "\\end{equation}\n", "\n", - "This means, that evaluation of the ACF at $\\kappa=0$ results in the quadratic mean, i.e. the mean power. For all other values of $\\kappa$, the ACF can only be less or equal than the mean power.\n", - "Due to the IDTFT definition, also the relation\n", + "This means that evaluation of the ACF at $\\kappa=0$ results in the quadratic mean, i.e., the mean power. For all other values of $\\kappa$, the ACF can only be less than or equal to the mean power.\n", + "Due to the IDTFT definition, the relation\n", "\n", "\\begin{equation}\n", "\\varphi_{xx}[\\kappa=0]=\\frac{1}{2\\pi}\\int\\limits_{-\\pi}^\\pi\\Phi_{xx}(\\Omega)\\mathrm{d}\\Omega\n", "\\end{equation}\n", "\n", - "holds for the mean power, i.e. it is proportional to the area of the power spectral density." + "holds for the mean power, i.e., it is proportional to the area of the power spectral density." ] }, { @@ -134,7 +134,7 @@ "\\mathrm{I:}\\quad\\varphi_{yx}[\\kappa]=\\varphi_{xx}[\\kappa]\\ast h[\\kappa] \\circ - \\bullet \\Phi_{yx}(\\Omega)=\\Phi_{xx}(\\Omega)\\cdot H(\\Omega).\n", "\\end{equation}\n", "\n", - "Interchanging of indices allows to write\n", + "Interchanging indices allows us to write\n", "\n", "\\begin{equation}\n", "\\varphi_{xy}[\\kappa]=\\varphi_{xx}[\\kappa]\\ast h^*[-\\kappa] \\circ - \\bullet \\Phi_{xy}(\\Omega)=\\Phi_{xx}(\\Omega)\\cdot H^*(\\Omega).\n", @@ -150,7 +150,7 @@ "\\Phi_{aa}(\\Omega)= A(\\Omega) \\cdot A^*(\\Omega) = |A(\\Omega)|^2,\n", "\\end{equation}\n", "\n", - "and the further useful relations (actually from these two everything else can be deduced)\n", + "and the further useful relations (actually, from these tw,o everything else can be deduced)\n", "\n", "\\begin{equation}\n", "\\begin{split}\n", @@ -161,25 +161,25 @@ "\n", "are important to show the links between the functions.\n", "\n", - "The Wiener-Lee theorem links the input, the system response and output ACFs as\n", + "The Wiener-Lee theorem links the input, the system response, and the output ACFs as\n", "\n", "\\begin{equation}\n", "\\mathrm{IV:}\\quad\\varphi_{yy}[\\kappa]=\\varphi_{xx}[\\kappa]\\ast \\varphi_{hh}[\\kappa] \\circ - \\bullet \\Phi_{yy}(\\Omega)=\\Phi_{xx}(\\Omega)\\cdot |H(\\Omega)|^2,\n", "\\end{equation}\n", "\n", - "which can be also written as\n", + "which can also be written as\n", "\n", "\\begin{equation}\n", "\\varphi_{yy}[\\kappa] = \\varphi_{xx}[\\kappa] \\ast \\varphi_{hh}[\\kappa] \\circ - \\bullet |Y(\\Omega)|^2=|X(\\Omega)|^2\\cdot |H(\\Omega)|^2.\n", "\\end{equation}\n", "\n", - "The auto correlation function\n", + "The autocorrelation function\n", "\\begin{equation}\n", "\\varphi_{hh}[\\kappa]=h[\\kappa]\\ast h^*[-\\kappa]\n", "\\end{equation}\n", "is called *filter ACF*.\n", "\n", - "Equations $\\mathrm{I}$, $\\mathrm{II}$, $\\mathrm{III}$ and $\\mathrm{IV}$ describe the fundamentals of random signal processing with LTI systems." + "Equations $\\mathrm{I}$, $\\mathrm{II}$, $\\mathrm{III}$, and $\\mathrm{IV}$ describe the fundamentals of random signal processing with LTI systems." ] }, { @@ -239,13 +239,13 @@ "\n", "Let us shortly consider continuous-time signals, since analytical calculus is a little bit simpler to show the intended key facts.\n", "\n", - "Calculate the continuous-time auto-correlation function\n", + "Calculate the continuous-time autocorrelation function\n", "\n", "\\begin{equation}\n", "\\phi_{xx}(\\tau)=\\lim\\limits_{T\\rightarrow\\infty}\\frac{1}{2T}\\int\\limits_{-T}^T x(t+\\tau)\\,x^*(t)\\,\\mathrm{d}t\n", "\\end{equation}\n", "\n", - "for the power signal (i.e. a signal with a finite mean power)\n", + "for the power signal (i.e., a signal with a finite mean power)\n", "\n", "\\begin{equation}\n", "x(t)=A\\cdot\\sin(\\omega\\,t+\\phi)\n", @@ -309,7 +309,7 @@ "x(t)=A\\cdot\\sin(\\omega\\,t+\\phi)\n", "\\end{equation}\n", "\n", - "has the auto correlation function\n", + "has the autocorrelation function\n", "\n", "\\begin{equation}\n", "\\phi_{xx}(\\tau) = \\frac{A^2}{2}\\cos(\\omega\\tau).\n", @@ -367,13 +367,13 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "- We see that maximum **positive** auto correlation is obtained when sine signals exactly overlay each other. This is the case for $\\tau= m \\frac{2\\pi}{\\omega}$, $m\\in\\mathbb{Z}$, i.e. multiples of the sine signal's period.\n", + "- We see that maximum **positive** autocorrelation is obtained when sine signals exactly overlay each other. This is the case for $\\tau= m \\frac{2\\pi}{\\omega}$, $m\\in\\mathbb{Z}$, i.e. multiples of the sine signal's period.\n", "\n", - "- The maximum **negative** auto correlation is obtained when the two signals exhibit exact inverted polarity.\n", - "For a single sine signal this is the case for a phase shift of 180 deg or given in time by half of the period, thus at $\\tau = m \\frac{\\pi}{\\omega}$, $|m|\\in 1,3,5,7,\\dots$.\n", + "- The maximum **negative** autocorrelation is obtained when the two signals exhibit exact inverted polarity.\n", + "For a single sine signal, this is the case for a phase shift of 180 deg or given in time by half of the period, thus at $\\tau = m \\frac{\\pi}{\\omega}$, $|m|\\in 1,3,5,7,\\dots$.\n", "\n", - "- Correlation is zero when the phase shift exhibits $\\pm$90 deg, i.e. a quarter of a period.\n", - "Then the autocorrelation is evaluated for a sine and a cosine function (the sign of the cosine depends on the actual shift) of same frequency.\n", + "- Correlation is zero when the phase shift exhibits $\\pm$90 deg, i.e,. a quarter of a period.\n", + "Then the autocorrelation is evaluated for a sine and a cosine function (the sign of the cosine depends on the actual shift) of the same frequency.\n", "Since these two signals are known to be orthogonal, the zero correlation is not surprising, cf. the [Fourier series](http://mathworld.wolfram.com/FourierSeries.html)." ] }, @@ -488,12 +488,12 @@ "\n", "Then we can state:\n", "- Period is 1000 samples, so phase shift of $\\frac{\\pi}{2}$ corresponds to 250 samples\n", - "- For $\\kappa=0$ we have a positive sine signal and a negative cosine signal to be correlated. Since these two sequences are orthogonal, the CCF is zero here.\n", - "- ACF is also periodic, it is a negative sine signal actually\n", - "- If $x$ is shifted by $\\kappa=-250$, i.e. to the right: signals lie on top of each other $\\to$ maximum value in the ACF.\n", - "- If $y$ is shifted by $\\kappa=-250$, i.e. to the left: signals lie on top of each other $\\to$ maximum value in the ACF.\n", - "- If $x$ is shifted by $\\kappa=+250$, i.e. to the left: signals are exactly out of phase $\\to$ minimum value in the ACF.\n", - "- If $y$ is shifted by $\\kappa=+250$, i.e. to the right: signals are exactly out of phase $\\to$ minimum value in the ACF.\n" + "- For $\\kappa=0$, we have a positive sine signal and a negative cosine signal to be correlated. Since these two sequences are orthogonal, the CCF is zero here.\n", + "- ACF is also periodic; it is a negative sine signal, actually\n", + "- If $x$ is shifted by $\\kappa=-250$, i.e., to the right: signals lie on top of each other $\\to$ the maximum value in the ACF.\n", + "- If $y$ is shifted by $\\kappa=-250$, i.e., to the left: signals lie on top of each other $\\to$ the maximum value in the ACF.\n", + "- If $x$ is shifted by $\\kappa=+250$, i.e., to the left: signals are exactly out of phase $\\to$ the minimum value in the ACF.\n", + "- If $y$ is shifted by $\\kappa=+250$, i.e., to the right: signals are exactly out of phase $\\to$ the minimum value in the ACF.\n" ] }, { @@ -550,7 +550,7 @@ "source": [ "The CCF is zero for all $\\kappa$.\n", "The signals correspond to harmonics in a Fourier series (when considering time-continuous signals) or a DFT.\n", - "In our case we could think of dealing with a DFT of length $N$, where we created sine signals for the \n", + "In our case, we could think of dealing with a DFT of length $N$, where we created sine signals for the \n", "5th and 20th DFT eigenfrequency.\n", "\n", "Recall: It is a key feature of DFT that the $\\mathrm{e}^{\\mathrm{j}\\frac{2\\pi}{N}k\\mu}$ DFT eigensignals are orthogonal, i.e.\n", @@ -583,7 +583,7 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "Subtask h) Generate a white noise signal of length $N$, whose amplitude values are drawn from the standard normal distribution (i.e. ideally $\\mu_x=0, \\sigma^2=1$, $\\varphi_{xx}=\\delta[\\kappa]$, $\\Phi_{xx}(\\Omega)=1$ holds)." + "Subtask h) Generate a white noise signal of length $N$, whose amplitude values are drawn from the standard normal distribution (i.e., ideally $\\mu_x=0, \\sigma^2=1$, $\\varphi_{xx}=\\delta[\\kappa]$, $\\Phi_{xx}(\\Omega)=1$ holds)." ] }, { @@ -603,7 +603,7 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "We see that the ACF is a Dirac-like function, not ideally a Dirac since we are dealing with finite length sequences and thus only an estimator for the ACF based on the random sample function $x=y$.\n", + "We see that the ACF is a Dirac-like function, not ideally a Dirac, since we are dealing with finite-length sequences and thus only an estimator for the ACF based on the random sample function $x=y$.\n", "\n", "We will deal with this Dirac property a lot in the subsequent exercises." ] @@ -647,7 +647,7 @@ "We make sure that we get the equation right.\n", "\n", "We have $N=3$.\n", - "A general time shift / lag by $\\kappa$ can be written as\n", + "A general time shift/lag by $\\kappa$ can be written as\n", "\n", "\\begin{equation}\n", "x[k+\\kappa] = \\delta[k+\\kappa] + 2 \\delta[k+\\kappa-1] -3 \\delta[k+\\kappa-2].\n", @@ -673,9 +673,9 @@ "\\varphi_{xx}[\\kappa=-2] = \\frac{1}{3} (\\delta[k-2] + 2 \\delta[k-2-1] -3 \\delta[k-2-2]) \\cdot (\\delta[k] + 2 \\delta[k-1] -3 \\delta[k-2]) = \\frac{1}{3}(1\\cdot (-3)) = \\frac{-3}{3} = -1\n", "\\end{equation}\n", "\n", - "Of course, the results are also obtained by drawing the two signals on a paper, shift one by a certain $\\kappa$, multiply the values that overlap and sum them up (just as we learned for convolution but without the time reversal step).\n", + "Of course, the results are also obtained by drawing the two signals on paper, shifting one by a certain $\\kappa$, multiplying the values that overlap, and summing them up (just as we learned for convolution, but without the time reversal step).\n", "\n", - "From above we know that\n", + "From above, we know that\n", "\n", "\\begin{equation}\n", "\\varphi_{xx}[\\kappa]=\\varphi^*_{xx}[-\\kappa]\n", @@ -699,9 +699,9 @@ "\\varphi_{xx}[\\kappa=+1] = -\\frac{4}{3},\n", "\\end{equation}\n", "\n", - "which could be also be derived when using the manual calculus with the Dirac sifting that was deployed for $\\kappa=0,-1,-2$.\n", + "which could also be derived when using the manual calculus with the Dirac sifting that was deployed for $\\kappa=0,-1,-2$.\n", "\n", - "For completeness we should state that \n", + "For completeness, we should state that \n", "\n", "\\begin{equation}\n", "\\varphi_{xx}[\\kappa] = 0\\quad\\mathrm{for}\\quad |\\kappa|>2.\n", @@ -751,7 +751,7 @@ "\\begin{equation}\n", "\\varphi_{xx}[\\kappa=0]=\\text{E}\\left\\{x[k]\\cdot x[k]\\right\\}=\\sigma_x^2+\\mu_x^2 = \\frac{14}{3}.\n", "\\end{equation}\n", - "Trivial fact, but worth to recapitulate: this is always a positive number." + "Trivial fact, but worth recapitulating: this is always a positive number." ] }, { @@ -780,7 +780,7 @@ "We assume that $x[k]$ is drawn from a stationary, ergodic random process.\n", "\n", "- Make a sketch of the amplitude response of the transfer function and characterize the system characteristics.\n", - "- Calculate the ACF, the linear mean and the variance of the output signal $y[k]$." + "- Calculate the ACF, the linear mean, and the variance of the output signal $y[k]$." ] }, { @@ -821,7 +821,7 @@ "source": [ "- ACF of output: \n", "\n", - "Make use of the Wiener-Lee theorem in DTFT domain\n", + "Make use of the Wiener-Lee theorem in the DTFT domain\n", "\n", "\\begin{equation}\n", "\\begin{split}\n", @@ -830,7 +830,7 @@ "\\end{split}\n", "\\end{equation}\n", "\n", - "The squared transfer function of LTI system is simply\n", + "The squared transfer function of an LTI system is simply\n", "\n", "\\begin{equation}\n", "|H(\\Omega)|^2=\\begin{cases}\n", @@ -839,7 +839,7 @@ "\\end{cases}\n", "\\end{equation}\n", "\n", - "Make use of inverse DTFT\n", + "Make use of the inverse DTFT\n", "\n", "\\begin{equation}\n", "\\begin{split}\n", @@ -871,15 +871,15 @@ "\n", "- Linear mean of output:\n", "\n", - "For the output signal it is fair to assume that\n", + "For the output signal, it is fair to assume that\n", "\n", "\\begin{equation}\n", "E\\{y[k] y[k+\\kappa]\\} = E\\{y[k]\\}\\cdot E\\{y[k+\\kappa]\\}\\quad\\mathrm{for}\\quad \\kappa\\to \\infty,\n", "\\end{equation}\n", "\n", "meaning that the signal $y$ for $k$ and the time shift $\\kappa$ are not correlated.\n", - "Since we applied a white noise that should exhibit this characteristics, the output signal fulfills this assumption as well, since it was only shaped by the ideal lowpass filter.\n", - "Furthermore the process is assumed to be stationary, which means that the mean $\\mu_y[k]$ is not time-dependent.\n", + "Since we applied a white noise that should exhibit these characteristics, the output signal fulfills this assumption as well, since it was only shaped by the ideal lowpass filter.\n", + "Furthermore, the process is assumed to be stationary, which means that the mean $\\mu_y[k]$ is not time-dependent.\n", "Thus,\n", "\\begin{equation}\n", "E\\{y[k]\\} = E\\{y[k+\\kappa]\\} = \\mu_y,\n", @@ -888,7 +888,7 @@ "\\begin{equation}\n", "\\varphi_{yy}[\\kappa] = E\\{y[k]\\}\\cdot E\\{y[k+\\kappa]\\} = \\mu_y^2\\quad\\mathrm{for}\\quad \\kappa\\to \\infty.\n", "\\end{equation}\n", - "Applying results in the final result\n", + "Applying results to the final result\n", "\\begin{equation}\n", "\\mu_y^2 = \\varphi_{yy}[\\kappa\\to \\infty]=\\lim_{\\kappa\\to\\infty}\\Phi_0\\frac{\\sin(\\Omega_c\\kappa)}{\\kappa\\pi} = \\frac{\\mathrm{bounded}}{\\infty} = 0.\n", "\\end{equation}\n", @@ -918,17 +918,17 @@ "## Interpretation\n", "\n", "- The Dirac ACF of input $\\varphi_{xx}[\\kappa]=\\Phi_0 \\delta[\\kappa]$ is shaped to the output ACF $\\varphi_{yy}[\\kappa] = \\Phi_0\\frac{\\sin(\\Omega_c\\kappa)}{\\kappa\\pi}$, i.e. a sinc-like shape.\n", - "This is typical for lowpass filter characteristics.\n", + "This is typical for low-pass filter characteristics.\n", "\n", "- The squared mean of input $\\mu_x^2 = \\varphi_{xx}[\\kappa\\to\\infty]=\\lim_{\\kappa\\to\\infty}\\Phi_0 \\delta[\\kappa]=0$ tells us that the input signal is mean free (DC free).\n", - "Since, the (time-independent) means of input/output signal are linked by $\\mu_y = \\mu_x \\cdot H(\\Omega)\\big|_{\\Omega=0}$, we obtain $\\mu_y=0$ here. Note that in this example, this result is independent of the specific $H(\\Omega)\\big|_{\\Omega=0}$.\n", + "Since the (time-independent) means of input/output signal are linked by $\\mu_y = \\mu_x \\cdot H(\\Omega)\\big|_{\\Omega=0}$, we obtain $\\mu_y=0$ here. Note that in this example, this result is independent of the specific $H(\\Omega)\\big|_{\\Omega=0}$.\n", "\n", - "- If $\\Omega_c=\\pi$, i.e. no low-pass characteristic at all, the complete power is transmitted through the system.\n", + "- If $\\Omega_c=\\pi$, i.e., no low-pass characteristic at all, the complete power is transmitted through the system.\n", "For this case, verify that $h[k] = \\delta[k]$, which furthermore implies $y = x$, $\\varphi_{yy} = \\varphi_{xx}$.\n", "\n", "- For $\\Omega_c<\\pi$, i.e. certain low-pass characteristic, only the fraction $\\frac{\\Omega_c}{\\pi}$ of the whole power $\\Phi_0$ is transmitted.\n", "\n", - "- Plot $\\varphi_{yy}[\\kappa] = \\Phi_0\\frac{\\sin(\\Omega_c\\kappa)}{\\kappa\\pi}$ by yourself for different $0\\leq \\Omega_c\\leq \\pi$. What should be expected and is this reflected by the plots?" + "- Plot $\\varphi_{yy}[\\kappa] = \\Phi_0\\frac{\\sin(\\Omega_c\\kappa)}{\\kappa\\pi}$ by yourself for different $0\\leq \\Omega_c\\leq \\pi$. What should be expected, and is this reflected by the plots?" ] }, { @@ -942,7 +942,7 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "This exercise is very similar to the exercise 4, just that the LTI system is slightly more complicated and we discuss the spectral domain in more detail. The key principles remain the same, since we are again dealing with white noise input signal." + "This exercise is very similar to Exercise 4, just that the LTI system is slightly more complicated, and we discuss the spectral domain in more detail. The key principles remain the same, since we are again dealing with a white noise input signal." ] }, { @@ -953,7 +953,7 @@ "\n", "Calculate the PSD $\\Phi_{yy}(\\Omega)$ of the output signal $y[k]=x[k]\\ast h[k]$. Make a sketch of the PSD $\\Phi_{yy}(\\Omega)$ for $0\\leq\\Omega<2\\pi$.\n", "\n", - "- $x[k]$ is a stationary random input signal for which the (ideal) auto-correlation function is given as\n", + "- $x[k]$ is a stationary random input signal for which the (ideal) autocorrelation function is given as\n", "\n", "\\begin{equation}\n", "\\varphi_{xx}[\\kappa]=\\sigma_x^2\\cdot\\delta[\\kappa]\n", @@ -965,7 +965,7 @@ "h[k]=\\left(\\frac{3}{4}\\right)^k\\cdot\\epsilon[k]\n", "\\end{equation}\n", "\n", - "In signals & systems we learned how to discuss the system characteristics in terms of impulse/step response, frequency response, pole/zero plot, bode plot.\n", + "In signals & systems, we learned how to discuss the system characteristics in terms of impulse/step response, frequency response, pole/zero plot, and Bode plot.\n", "We will need some of this stuff here again." ] }, @@ -975,7 +975,7 @@ "source": [ "## Solution\n", "\n", - "We make use of the Wiener-Lee theorem, introduce the filter-ACF and find the links in spectral (DTFT) domain:\n", + "We make use of the Wiener-Lee theorem, introduce the filter-ACF, and find the links inthe spectral (DTFT) domain:\n", "\n", "\\begin{equation}\n", "\\varphi_{yy}[\\kappa]=\\varphi_{hh}[\\kappa]\\ast\\varphi_{xx}[\\kappa]=(h[\\kappa]\\ast h^*[-\\kappa])\\ast\\varphi_{xx}[\\kappa] \\circ - \\bullet \\Phi_{yy}(\\Omega)=|H(\\Omega)|^2\\cdot\\Phi_{xx}(\\Omega).\n", @@ -989,8 +989,8 @@ "\\Phi_{xx}(\\Omega)=\\text{DTFT}\\{\\varphi_{xx}[\\kappa]\\}=\\text{DTFT}\\{\\sigma_x^2\\cdot\\delta[\\kappa]\\}=\\sigma_x^2.\n", "\\end{equation*}\n", "\n", - "The squared magnitude of the transfer function can be derived by the DTFT of the impulse response $h[k]$.\n", - "The well known relation (have a look at signals and systems lecture / exercise) \n", + "The squared magnitude of the transfer function can be derived from the DTFT of the impulse response $h[k]$.\n", + "The well-known relation (have a look at the signals and systems lecture/exercise) \n", "\n", "\\begin{equation}\n", "\\text{DTFT}\\{a^k\\cdot\\epsilon[k]\\}=\\frac{1}{1-a\\cdot\\mathrm{e}^{-\\mathrm{j}\\Omega}} \\quad\\text{for}\\quad|a|<1\n", @@ -1029,15 +1029,15 @@ "\\end{split}\n", "\\end{equation}\n", "\n", - "We should plot this normalized by $\\sigma_x^2$, so in fact we plot $|H(\\Omega)|^2$ then for our example.\n", - "This also gives us a nice way to validate our analytic calculus. Instead of the DTFT, we can describe the LTI system response in the $z$-domain, i.e. with its $z$-transform\n", + "We should plot this normalized by $\\sigma_x^2$, so in fact we plot $|H(\\Omega)|^2$, then for our example.\n", + "This also gives us a nice way to validate our analytic calculus. Instead of the DTFT, we can describe the LTI system response in the $z$-domain, i.e., with its $z$-transform\n", "\n", "\\begin{equation}\n", "h[k] \\circ - \\bullet H(z) = \\frac{1}{1-\\frac{3}{4} z^{-1}},\n", "\\end{equation}\n", "\n", "where we recognize a recursive system with difference equation $a_0 y[k] = b_0 x[k] - a_1 y[k-1]$ using the coefficients $b_0=1$, $a_0=1$ and $a_1=-\\frac{3}{4}$.\n", - "The `freqz()` function conveniently allows to evaluate the magnitude response of (non)-recursive systems, i.e. it evaluates the DTFT numerically with certain frequency resolution. Here we use $\\Delta\\Omega = \\frac{2\\pi}{256}$." + "The `freqz()` function conveniently allows for evaluating the magnitude response of (non)-recursive systems, i.e., it evaluates the DTFT numerically with certain frequency resolution. Here we use $\\Delta\\Omega = \\frac{2\\pi}{256}$." ] }, { @@ -1074,13 +1074,13 @@ "source": [ "## Interpretation\n", "\n", - "- In the plot we see, that analytic and numeric solution of $|H(\\Omega)|^2$ are the same, so our calculus is very likely correct if we assume that freqz() is correctly implemented.\n", + "- In the plot, we see that the analytic and numeric solutions of $|H(\\Omega)|^2$ are the same, so our calculus is very likely correct if we assume that freqz() is correctly implemented.\n", "\n", - "- Furthermore in the plot, we sketched the $\\Phi_{xx}(\\Omega)$ normalized by $\\sigma_x^2$, which yields a line, i.e. magnitude one for all frequencies.\n", + "- Furthermore, in the plot, we sketched the $\\Phi_{xx}(\\Omega)$ normalized by $\\sigma_x^2$, which yields a line, i.e., magnitude one for all frequencies.\n", "\n", - "- We see that the input PSD $\\Phi_{xx}(\\Omega)$, which is constant over frequency, is spectrally shaped by the LTI system towards the output PSD $\\Phi_{yy}(\\Omega)$ .\n", + "- We see that the input PSD $\\Phi_{xx}(\\Omega)$, which is constant over frequency, is spectrally shaped by the LTI system towards the output PSD $\\Phi_{yy}(\\Omega)$.\n", "\n", - "- The output PSD $\\Phi_{yy}(\\Omega)$ is thus not longer constant over frequency, but amplified at lower frequencies and attenuated at higher frequencies.\n", + "- The output PSD $\\Phi_{yy}(\\Omega)$ is thus no longer constant over frequency, but amplified at lower frequencies and attenuated at higher frequencies.\n", "This filter exhibits certain lowpass characteristics with a DC gain of 12 dB, -3 dB cut frequency approximately at $\\Omega=\\frac{1}{10} \\pi$ and unity gain (0 dB) at about $\\frac{4}{10} \\pi$.\n" ] }, @@ -1104,11 +1104,11 @@ "source": [ "## Task\n", "\n", - "In this programming task we want to elaborate how to identify the impulse response of an LTI system with random signal in the time domain.\n", - "If white noise is used as input signal, the task becomes very convenient in terms of required signal processing steps. Thus, we\n", + "In this programming task, we want to elaborate on how to identify the impulse response of an LTI system with a random signal in the time domain.\n", + "If white noise is used as an input signal, the task becomes very convenient in terms of required signal processing steps. Thus, we\n", "\n", "- generate white noise signal $x[k]$ drawn from gaussian PDF\n", - "- create finite impulse response $h[k]$ of a simple LTI system, i.e. a lowpass (in practice this would be unknown)\n", + "- create a finite impulse response $h[k]$ of a simple LTI system, i.e., a lowpass (in practice, this would be unknown)\n", "- apply convolution $y[k] = x[k] \\ast h[k]$\n", "- estimate the impulse response $\\hat{h}[k]$ based on the concept of correlation functions." ] @@ -1126,12 +1126,12 @@ "source": [ "### Generate White Noise \n", "\n", - "We aim at an ACF\n", + "We aim for an ACF\n", "\n", "$$\\varphi_{xx}[\\kappa]\\approx \\delta[\\kappa]$$\n", "\n", "and thus constant (Auto)-PSD for an input signal block $x[k]$ that was drawn\n", - "- from Gaussian PDF with zero mean and unit variance (i.e. standard normal distribution)\n", + "- from a Gaussian PDF with zero mean and unit variance (i.e., standard normal distribution)\n", "- with N=$2^{16}$ samples length" ] }, @@ -1174,7 +1174,7 @@ "- cut frequency 4.8 kHz\n", "- sampling frequency 48 kHz\n", "- 45 taps\n", - "- one tap weight is manipulated manually to break axial-symmetry on purpose" + "- one tap weight is manipulated manually to break axial symmetry on purpose" ] }, { @@ -1246,15 +1246,15 @@ "source": [ "## Estimate the Impulse Response by Cross Correlation\n", "\n", - "The (most simple, and in practice rather not longer used) method here is based on the fact that the ACF of the input signal is ideally a Dirac impulse.\n", + "The (most simple, and in practice, rather not longer used) method here is based on the fact that the ACF of the input signal is ideally a Dirac impulse.\n", "\n", - "In essence we deal with:\n", + "In essence, we deal with:\n", "\n", "- CCF of $y[k]$ and $x[k]$ yields **either** $\\varphi_{yx}$ **or** $\\varphi_{xy}$ (note the index exchange) depending on the order of calling `my_xcorr2()`\n", "- check $\\varphi_{yx}[\\kappa]=h[+\\kappa]$\n", "- check $\\varphi_{xy}[\\kappa]=h[-\\kappa]$\n", "\n", - "More detailed, the relationships between CCF, ACF and IR are known from above\n", + "More detailed, the relationships between CCF, ACF, and IR are known from above\n", "\n", "\\begin{equation}\n", "\\begin{split}\n", @@ -1268,7 +1268,7 @@ "\\varphi_{xy}[\\kappa]=\\varphi_{xx}[\\kappa]\\ast h[-\\kappa].\n", "\\end{equation}\n", "\n", - "For above assumption that the random input signal fulfills\n", + "For the above assumption that the random input signal fulfills\n", "\n", "$$\\varphi_{xx}[\\kappa]\\approx \\delta[\\kappa]$$\n", "\n", @@ -1280,7 +1280,7 @@ "\\end{split}\n", "\\end{equation}\n", "\n", - "and since the Dirac impulse behaves neutral in convolution we have\n", + "and since the Dirac impulse behaves neutral in convolution, we have\n", "\n", "\\begin{equation}\n", "\\begin{split}\n", @@ -1343,12 +1343,12 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "By interchanging the order for cross correlation, i.e. `my_xcorr2(x,y)`, \n", + "By interchanging the order for cross correlation, i.e., `my_xcorr2(x,y)`, \n", "\n", "$$\\varphi_{xy}[\\kappa]=h[-\\kappa]$$\n", "\n", "we obtain the mirrored impulse response.\n", - "To conveniently see this, we implemented the `h[30] = 0` manipulation, to make slightly IR asymmetric." + "To conveniently see this, we implemented the `h[30] = 0` manipulation to make it slightly IR asymmetric." ] }, { @@ -1375,7 +1375,7 @@ "source": [ "## Interpretation\n", "\n", - "The estimated impulse can be compared against the reference impulse response (we only can do this here, since we have built one for demonstration. In practice this response is of course unknown, and one needs to be sure about the chosen method to obtain a valid estimation.). \n", + "The estimated impulse can be compared against the reference impulse response (we can only do this here, since we have built one for demonstration. In practice, this response is of course unknown, and one needs to be sure about the chosen method to obtain a valid estimation.). \n", "We can take the absolute error on the impulse response itself, here as\n", "\n", "\\begin{equation}\n", @@ -1420,14 +1420,14 @@ "source": [ "Without an application in mind, we cannot tell if these errors are sufficiently small or too large.\n", "\n", - "However, we can observe for that specific example, that the error $e_h$ somehow follows the shape of the impulse response, i.e. a sinc-like pattern.\n", - "Furthermore, we can state that for low frequencies there is a constant offset about 0.5 dB in the magnitude spectrum, slowly decreasing to zero dB in the region of the cut frequency. Above the cut frequency (except the strong peak) the error is somehow 'noisy' around 0 dB.\n", + "However, we can observe for that specific example that the error $e_h$ somehow follows the shape of the impulse response, i.e., a sinc-like pattern.\n", + "Furthermore, we can state that for low frequencies, there is a constant offset of about 0.5 dB in the magnitude spectrum, slowly decreasing to zero dB in the region of the cut frequency. Above the cut frequency (except for the strong peak), the error is somehow 'noisy' around 0 dB.\n", "\n", - "We might play around with the length of the input signal and the specific seed of the random number generator, to obtain different pictures.\n", - "For example instead of `N = 2**16` check what happens with `N = 2**8`.\n", - "It is very likely that the estimated impulse response is not precise enough for any application, since the deviations to the reference are too large.\n", + "We might play around with the length of the input signal and the specific seed of the random number generator to obtain different pictures.\n", + "For example, instead of `N = 2**16`, check what happens with `N = 2**8`.\n", + "It is very likely that the estimated impulse response is not precise enough for any application, since the deviations from the reference are too large.\n", "\n", - "Make yourself comfortable, why and how the parameter `N` changes the result and how the assumption of *ACF of $x$ being a Dirac* is affected. " + "Make yourself comfortable, why and how the parameter `N` changes the result, and how the assumption of *ACF of $x$ being a Dirac* is affected. " ] }, { @@ -1456,7 +1456,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.10.12" + "version": "3.12.2" } }, "nbformat": 4, diff --git a/spectral_estimation/stft.ipynb b/spectral_estimation/stft.ipynb index f88030e..010b509 100755 --- a/spectral_estimation/stft.ipynb +++ b/spectral_estimation/stft.ipynb @@ -20,7 +20,7 @@ "- lecture: https://github.com/spatialaudio/digital-signal-processing-lecture\n", "- tutorial: https://github.com/spatialaudio/digital-signal-processing-exercises\n", "\n", - "Feel free to contact lecturer jacob.thoenes@uni-rostock.de" + "Feel free to contact the lecturer jacob.thoenes@uni-rostock.de" ] }, { @@ -35,14 +35,14 @@ "With the STFT, the time-frequency varying characteristics of a signal can be analyzed.\n", "Combining the STFT as the analyzing step and a subsequent iSTFT as the synthesizing/reconstruction step, builds the foundation for so called **time-frequency processing**, which is heavily used in nowadays DSP algorithms.\n", "\n", - "For certain algorithms the computational load can be considerably decreased using time-frequency processing.\n", + "For certain algorithms, the computational load can be considerably decreased using time-frequency processing.\n", "Other approaches can only be realized due to explicit usage of time-frequency processing.\n", "A very prominent example is the invention of [Auto-Tune](https://en.wikipedia.org/wiki/Auto-Tune) in the late 1990s. Furthermore, all video and audio codecs work this way.\n", "IP-based streaming is only possible due to the transport of compressed time-frequency data rather than raw data. \n", "\n", - "If we consider only the analysis stage, i.e. the STFT itself, we can elaborate a nice link to the periodogram: we can interpret this as a special post-processing of STFT data.\n", + "If we consider only the analysis stage, i.e., the STFT itself, we can elaborate a nice link to the periodogram: we can interpret this as a special post-processing of STFT data.\n", "\n", - "It thus appears reasonable that we have a detailed look to STFT/iSTFT and their key concept within this homework assignment." + "It thus appears reasonable that we have a detailed look at STFT/iSTFT and their key concept within this homework assignment." ] }, { @@ -56,11 +56,11 @@ "- programming the Short-time Fourier Transform (STFT) by ourselves\n", "- programming the Inverse Short-time Fourier Transform (iSTFT) by ourselves\n", " - using standard overlap-add (OLA)\n", - " - using weigthed overlap-add (WOLA)\n", + " - using weighted overlap-add (WOLA)\n", "- normalize the STFT data either as a spectrum or as a power spectral density (PSD)\n", "- make a surface plot of the magnitude of the STFT spectrum\n", "- calculate PSDs from the STFT data\n", - "- apply the Welch method, i.e. averaging PSDs, to obtain a so called modified periodogram" + "- apply the Welch method, i.e., averaging PSDs, to obtain a so called modified periodogram" ] }, { @@ -198,11 +198,11 @@ "\n", "- $-\\infty \\leq n\\in\\mathbb{Z} \\leq +\\infty$ is the sample index for whole time range\n", "- $k\\in\\mathbb{Z}$ is --- as we usually define --- the sample index for a DFT block of length $N$, so $0\\leq k \\leq N-1$\n", - "- $N\\in\\mathbb{N}$ ist the DFT block size, here we use sizes corresponding to a power of two\n", + "- $N\\in\\mathbb{N}$ is the DFT block size, here, we use sizes corresponding to a power of two\n", "- $\\mu\\in\\mathbb{Z}$ is --- as we usually define --- the index for discrete DFT frequencies, $0\\leq \\mu \\leq N-1$\n", - "- $\\theta\\in\\mathbb{Z}$ is the index for discrete time blocks of the STFT, for ease of discussion we consider $n \\geq 0$ and indicate the very first time block with $\\theta=0$, thus $\\theta\\geq 0$; the mnemonic *theta* for *time* might help\n", - "- $H\\in\\mathbb{N}>0$ is the so called hop size, i.e. by how many samples we jump to get to the next block\n", - "- in software functions often the number of overlapped samples $N-H$ must be defined rather than the hop size" + "- $\\theta\\in\\mathbb{Z}$ is the index for discrete time blocks of the STFT, for ease of discussion, we consider $n \\geq 0$ and indicate the very first time block with $\\theta=0$, thus $\\theta\\geq 0$; the mnemonic *theta* for *time* might help\n", + "- $H\\in\\mathbb{N}>0$ is the so called hop size, i.e., by how many samples we jump to get to the next block\n", + "- in software functions, often the number of overlapped samples $N-H$ must be defined rather than the hop size" ] }, { @@ -213,7 +213,7 @@ "\n", "Fig. 1 depicts the fundamental signal processing chain for time-frequency processing. We will have a detailed look at the STFT and iSTFT itself. So we need to discuss the analysis and the synthesis stage.\n", "The theoretical foundations of the STFT analysis and synthesis were derived in the late 1970s and are referenced in the bibliography listing at the end of this document.\n", - "Keep this in mind when asking for detailed knowledge about STFT some time.\n", + "Keep this in mind when asking for detailed knowledge about STFT sometimes.\n", "Particularly helpful might be chapter 7 in [Zoe11] as a convenient summary of key concepts. \n", "\n", "\n", @@ -239,7 +239,7 @@ "\n", "The notation of using finite length DFT blocks inherently implies \n", "$w[k] = 0\\,\\,\\,\\text{for}\\,\\,\\,k<0\\,\\,\\,\\text{and}\\,\\,\\,k>N-1$.\n", - "We will use windows that have none zeropadding parts over the length $N$.\n", + "We will use windows that have no zero-padding parts over the length $N$.\n", "We should follow the convention that in a **row** the **time varies**, and in a **column** the **frequency changes** when aligning $X[\\mu,\\theta]$ as two-dimensional, complex-valued data matrix $\\mathbf{X}_{\\mu,\\theta}$.\n", "Thus, we have access to complex-valued data points aligned over frequency and time indices.\n", "This is the building block **STFT**, cf. the upper part of Fig. 1." @@ -264,9 +264,9 @@ "source": [ "- Recall that the obtained time-domain signal blocks $x_\\theta = x[\\theta\\cdot H+k] \\cdot w[k]$ hold for $0\\leq k\\leq N-1$. Note that we calculate $x_q$ blocks rather than $x[\\theta\\cdot H+k]$ in the iSTFT stage. Proper alignment and superposition of all $x_q$ blocks then returns the signal $x[n]$ under certain conditions, which we will invent below as so called perfect reconstruction. \n", "- Absolute time index is $n=\\theta\\cdot H+k$.\n", - "- Blocks might overlap, when the hop size $HN$.\n", - "In such a case certain time ranges are completely missing and no perfect reconstruction can be achieved at all.\n", + "In such a case, certain time ranges are completely missing, and no perfect reconstruction can be achieved at all.\n", "- The case $H=N$ is sometimes used for the Bartlett and Welch periodogram." ] }, @@ -280,16 +280,16 @@ "This is depicted in the lower part of Fig. 1 and is known as **overlap-add** (OLA).\n", "The OLA is already known from [segmented convolution](https://nbviewer.jupyter.org/github/spatialaudio/digital-signal-processing-lecture/blob/master/nonrecursive_filters/segmented_convolution.ipynb).\n", "Here it is used exactly in the same manner: we have blocks aligned at different time instances, these blocks might overlap, and their superposition gives us the desired output signal.\n", - "We should have an initial expectation that an analysis/synthesis-pair for perfect reconstruction needs carefully chosen parameters, such as\n", + "We should have an initial expectation that an analysis/synthesis pair for perfect reconstruction needs carefully chosen parameters, such as\n", "\n", "- window characteristics\n", - "- window length / block size\n", + "- window length/block size\n", "- hop size\n", "- signal pre-/post-processing steps\n", "\n", "The STFT literature reports **three different** concepts on how to properly approach OLA:\n", "\n", - "1. the **basic OLA** (commonly abbreviated with OLA), this is the approach discussed with formulas above so far, considered in [Smi01], [OS10], [Rei15], in Matlab this is an option for the `istft()`\n", + "1. the **basic OLA** (commonly abbreviated with OLA), this is the approach discussed with formulas above so far, considered in [Smi01], [OS10], [Rei15], in Matlab, this is an option for the `istft()`\n", "\n", "2. the **weighted OLA** (WOLA), this is what Fig. 1 actually shows, considered in [Smi01], [Zoe11], [Pul18], this is Matlab's default for `istft()`\n", "\n", @@ -298,7 +298,7 @@ "See [Gri84, eq. (6),(7),(8)] for a detailed elaboration.\n", "\n", "In order to reinvent the ideas of these three concepts (the proofs of them are more sophisticated, but the following simple example works quite well), consider the trivial input signal $x[n]=1$ for all $n$.\n", - "The STFT and subsequent iSTFT yields the rectangular signal shifted to different time instances and windowed\n", + "The STFT and subsequent iSTFT yield the rectangular signal shifted to different time instances and windowed\n", "\n", "\\begin{align}\n", "\\underbrace{\\mathrm{rect}_N[\\theta\\cdot H+k] \\, w[k]}_{x_\\theta} = \\frac{1}{N}\\sum_{\\mu=0}^{N-1}X[\\mu,\\theta]\\cdot\\mathrm{e}^{+\\mathrm{j}\\frac{2\\pi}{N}k\\mu}.\n", @@ -311,19 +311,19 @@ "source": [ "#### Basic OLA\n", "\n", - "For OLA we need to sum these blocks as \n", + "For OLA, we need to sum these blocks as \n", "\\begin{align}\n", "\\sum_\\theta \\mathrm{rect}_N[\\theta\\cdot H+k] \\, w[k] =\\sum_\\theta \\frac{1}{N}\\sum_{\\mu=0}^{N-1}X[\\mu,\\theta]\\cdot\\mathrm{e}^{+\\mathrm{j}\\frac{2\\pi}{N}k\\mu}.\n", "\\end{align}\n", - "The question for aiming at perfect reconstruction is: Under which circumstances this equation becomes equal to $x[n] = 1$ ?!?\n", + "The question for aiming at perfect reconstruction is: Under which circumstances does this equation become equal to $x[n] = 1$ ?!?\n", "\n", - "Recall, that the window is non-zero only for its arguments $[0...N-1]$.\n", + "Recall that the window is non-zero only for its arguments $[0...N-1]$.\n", "Then, we can get rid of the (unit amplitude) rectangular function and deploy shifted windows directly.\n", "This leads to the condition (absolute time index $n=\\theta\\cdot H+k$ rewritten as $k = n - \\theta\\cdot H$)\n", "\\begin{align}\n", "\\sum_\\theta w[n -\\theta\\cdot H] = 1\n", "\\end{align}\n", - "for perfect reconstruction of our trivial example signal $x[n]=1$.\n", + "for the perfect reconstruction of our trivial example signal $x[n]=1$.\n", "We can generalize this condition for any signal $x[n]$.\n", "It simply means that all overlapping windows must sum to unity. Or, at least they must sum to $\\sum_\\theta w[n -\\theta\\cdot H] = \\mathrm{const}$, by which the reconstructed signal can be normalized then.\n", "This is known as the constant-overlap-add (**[COLA](https://ccrma.stanford.edu/~jos/sasp/Mathematical_Definition_STFT.html#19930)**) constraint for windows." @@ -335,7 +335,7 @@ "source": [ "#### Weighted OLA\n", "\n", - "The WOLA adds **additional windowing after iSTFT**. For ease of discussion we assume this to be the same window as for the analysis stage. \n", + "The WOLA adds **additional windowing after iSTFT**. For ease of discussion, we assume this to be the same window as for the analysis stage. \n", "This yields\n", "\n", "\\begin{align}\n", @@ -365,7 +365,7 @@ "\n", "As an improvement of the WOLA for\n", "\n", - "- windows that not exactly fulfill the WOLA criterion and\n", + "- windows that do not exactly fulfill the WOLA criterion and\n", "- more generally, to find the least squares approximation of a desired STFT spectrum\n", "\n", "[Gri84] showed that\n", @@ -375,15 +375,15 @@ "\\end{align}\n", "\n", "yields the optimum solution.\n", - "So, instead of relying on the WOLA criterion (overlapping of squared windows must be constant), the reconstructed signal (numerator, WOLA approach) is divided by the window-overlap signal (denominator). However, this requires $\\sum_\\theta w^2[n -\\theta\\cdot H]\\neq 0$ for any sample, which is known as nonzero overlap-add (**[NOLA](https://docs.scipy.org/doc/scipy/reference/generated/scipy.signal.check_NOLA.html)**) constraint.\n", - "This approach is implemented in `scipy.signal.istft()`. This is convenient, since we don't have to care if WOLA criterion is met, but just need to check if NOLA holds." + "So, instead of relying on the WOLA criterion (overlapping of squared windows must be constant), the reconstructed signal (numerator, WOLA approach) is divided by the window-overlap signal (denominator). However, this requires $\\sum_\\theta w^2[n -\\theta\\cdot H]\\neq 0$ for any sample, which is known as the nonzero overlap-add (**[NOLA](https://docs.scipy.org/doc/scipy/reference/generated/scipy.signal.check_NOLA.html)**) constraint.\n", + "This approach is implemented in `scipy.signal.istft()`. This is convenient, since we don't have to care if the WOLA criterion is met, but just need to check if NOLA holds." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "In the first task section we want to implement our own STFT / iSTFT-pair for certain configurations.\n", + "In the first task section, we want to implement our own STFT / iSTFT-pair for certain configurations.\n", "These will then be investigated in terms of their COLA and WOLA constraints." ] }, @@ -404,7 +404,7 @@ "x[n] = \\cos(\\Omega_N \\, n) - \\frac{1}{2}\\sin(2\\, \\Omega_N \\, n)\n", "\\end{equation}\n", "\n", - "in the time range $0 \\leq n \\leq L-1$ using the fundamental digital frequency $\\Omega_N = \\frac{2\\pi}{N}$. Plot the signal over $n$ with proper axis labeling, x-ticks labeling with step size of 16 or 32 might be helpful." + "in the time range $0 \\leq n \\leq L-1$ using the fundamental digital frequency $\\Omega_N = \\frac{2\\pi}{N}$. Plot the signal over $n$ with proper axis labeling, x-ticks labeling with a step size of 16 or 32 might be helpful." ] }, { @@ -463,7 +463,7 @@ "It is explicitly required to **implement the STFT by our own**.\n", "\n", "For simplified programming, we might assume that $L/H\\in\\mathbb{N}$. Of course, the usage of the FFT (in Matlab: `fft()`, in Python: `numpy.fft.fft()` or `scipy.fft.fft()`) is allowed and meaningful.\n", - "We also might validate our STFT results against build-in functions.\n", + "We also might validate our STFT results against built-in functions.\n", "\n", "So, in this subtask 2.2, we create four STFT matrices in total:\n", "\n", @@ -475,7 +475,7 @@ "\n", "3. What are the dimensions of the STFT matrices?\n", "\n", - "4. We assume a sampling frequency $f_s=4$ Hz. For that state the physical time and frequency vectors, i.e. all values of $\\mu \\Delta f$ and $\\theta \\Delta t$. Recall that $\\mu$ and $\\theta$ start with zero by our convention.\n" + "4. We assume a sampling frequency $f_s=4$ Hz. For that state, the physical time and frequency vectors, i.e., all values of $\\mu \\Delta f$ and $\\theta \\Delta t$. Recall that $\\mu$ and $\\theta$ start with zero by our convention.\n" ] }, { @@ -516,12 +516,12 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "Note: We might consider to implement tasks 3 to 6 as a single block.\n", + "Note: We might consider implementing tasks 3 to 6 as a single block.\n", "Note that certain functionalities are repeatedly required, which could be conveniently encapsulated in functions.\n", "\n", "# Task 3: Inverse STFT\n", "\n", - "Implement own written source code that calculates the inverse STFT from the matrix $\\mathbf{X}_{\\mu,\\theta}$, i.e. all possible blocks $x_\\theta$ according to\n", + "Implement own written source code that calculates the inverse STFT from the matrix $\\mathbf{X}_{\\mu,\\theta}$, i.e., all possible blocks $x_\\theta$ according to\n", "\n", "\\begin{align}\n", "\\underbrace{x[\\theta\\cdot H+k] \\cdot w[k]}_{x_\\theta}= \\frac{1}{N}\\sum_{\\mu=0}^{N-1}X[\\mu,\\theta]\\cdot\\mathrm{e}^{+\\mathrm{j}\\frac{2\\pi}{N}k\\mu}\n", @@ -573,7 +573,7 @@ "source": [ "# Task 6: Check COLA Conditions\n", "\n", - "Implement own written source code to check the COLA / WOLA criteria\n", + "Implement your own written source code to check the COLA / WOLA criteria\n", "\n", "\\begin{align}\n", "c_\\mathrm{OLA}[n] = \\sum_\\theta w[n -\\theta\\cdot H] = \\mathrm{const},\n", @@ -669,7 +669,7 @@ "source": [ "# Task 7: Plot Signals\n", "\n", - "We created some signals which we now should plot to inspect and validate.\n", + "We created some signals that we now should plot to inspect and validate.\n", "Also, feel free to use plots of intermediate results as well.\n", "See the figure with the four subplots visualizing the four cases.\n", "\n", @@ -681,7 +681,7 @@ "- the cola signals $c_\\mathrm{OLA}[n]$ and $c_\\mathrm{WOLA}[n]$\n", "- the error signal $x_\\mathrm{reconstructed}[n] - x[n]$, e.g. `xa-x`\n", "\n", - "Pay attention to properly given legend, title, axis ticks and labeling, as well as meaningful plotting ranges." + "Pay attention to properly given legend, title, axis ticks, and labeling, as well as meaningful plotting ranges." ] }, { @@ -770,7 +770,7 @@ "\n", "Graphical visualization always requires interpretation.\n", "\n", - "One observation would be the fade in and fade out at the beginning and the end of the considered finite length.\n", + "One observation would be the fade-in and fade-out at the beginning and the end of the considered finite length.\n", "\n", "- Why is the case?\n", "\n", @@ -781,7 +781,7 @@ "- In which cases a,b,c,d perfect reconstruction is achieved **without** any further post-processing?\n", "- In which cases a,b,c,d perfect reconstruction is achieved **only with** further post-processing?\n", "- Name possible post-processing algorithms.\n", - "- Which other windows could have been used to obtain similar results. Why?" + "- Which other windows could have been used to obtain similar results? Why?" ] }, { @@ -790,9 +790,9 @@ "source": [ "# What we got so far...\n", "\n", - "By now, we have implemented the basic structure of time-frequency processing and hopefully got some detailed insights of the data structure of an STFT and the different windowing concepts.\n", + "By now, we have implemented the basic structure of time-frequency processing and hopefully gained some detailed insights into the data structure of an STFT and the different windowing concepts.\n", "\n", - "Note, that our above analysis and synthesis/reconstruction assumes **unmodified STFT data**.\n", + "Note that our above analysis and synthesis/reconstruction assume **unmodified STFT data**.\n", "\n", "For time-frequency processing, certain modifications in magnitude and phase (such as filtering) are vivid parts of the DSP job. Then the signal processing needs some **further processing steps**. These can be seen in Fig. 2 and comprise of\n", "\n", @@ -800,8 +800,8 @@ "- additional phase shift, such that the phase in each STFT frame refers to the original starting time of the input signal. This is especially required if we want to modify the phase\n", "\n", "before time-frequency processing, and afterwards the reversed operations.\n", - "Furthermore, to avoid time-aliasing often the blocks are zeropadded prior to the FFT. This is related to [linear convolution using cyclic convolution](https://nbviewer.jupyter.org/github/spatialaudio/digital-signal-processing-lecture/blob/master/nonrecursive_filters/fast_convolution.ipynb).\n", - "Please, be aware of this and elaborate on these facts once required sometime." + "Furthermore, to avoid time-aliasing, often the blocks are zeropadded prior to the FFT. This is related to [linear convolution using cyclic convolution](https://nbviewer.jupyter.org/github/spatialaudio/digital-signal-processing-lecture/blob/master/nonrecursive_filters/fast_convolution.ipynb).\n", + "Please, be aware of this and elaborate on these facts when required." ] }, { @@ -947,7 +947,7 @@ "source": [ "# STFT Spectrum\n", "\n", - "As a second task section, we should investigate the magnitude spectrum of an STFT as a surface plot, i.e. colored magnitude over time and frequency axis in one graph. We approach this by programing our own example using our own STFT. Once we've mastered that, we can have a closer look to other examples on our own." + "As a second task section, we should investigate the magnitude spectrum of an STFT as a surface plot, i.e., colored magnitude over time and frequency axis in one graph. We approach this by programming our own example using our own STFT. Once we've mastered that, we can have a closer look at other examples on our own." ] }, { @@ -971,9 +971,9 @@ "\\end{cases}\n", "\\end{align}\n", "\n", - "- Write own source code to generate the signal, store it in variable `x`.\n", + "- Write your own source code to generate the signal, store it in the variable `x`.\n", "\n", - "Before actual programming, we should ask what this signal actually consists of (this is a typical and of course intended example, where the analytic equation looks more complicated than a concise implementation in software, in Python this probably could be done with a one-liner). The following thoughts might help:\n", + "Before actual programming, we should ask what this signal actually consists of (this is a typical and, of course, intended example, where the analytic equation looks more complicated than a concise implementation in software; in Python, this probably could be done with a one-liner). The following thoughts might help:\n", "\n", "- What do $L$, $M$ and $N$ actually do?\n", "- Why do we use $\\mu$ and $n$ variables just as in a DFT?\n", @@ -1014,10 +1014,10 @@ "source": [ "# Task 10: Generate Gaussian-White Noise\n", "\n", - "- For the same signal length as $x[n]$ / `x` generate a random noise signal drawn from a normal distribution with zero mean and standard deviation $\\sigma = \\frac{3}{2}$. Store the data in signal vector `noise`.\n", - "- State the unbiased estimates of the mean, the standard deviation and the variance for `noise`. \n", + "- For the same signal length as $x[n]$ / `x`, generate a random noise signal drawn from a normal distribution with zero mean and standard deviation $\\sigma = \\frac{3}{2}$. Store the data in the signal vector `noise`.\n", + "- State the unbiased estimates of the mean, the standard deviation, and the variance for `noise`. \n", "\n", - "If we use `numpy` we might set up `np.random.seed(seed=1)` to obtain same results as referenced below. " + "If we use `numpy`, we might set up `np.random.seed(seed=1)` to obtain the same results as referenced below. " ] }, { @@ -1050,8 +1050,8 @@ "# Task 11: Mixing Input Signal with Noise Signal\n", "\n", "- Add `x` and `noise` and store the result as `xn`.\n", - "- What are the theoretical values for the resulting mean, standard deviation and variance of this signal mixture?\n", - "- State the unbiased estimates of the mean, the standard deviation and variance for `xn`. Do they fulfill our expectation?" + "- What are the theoretical values for the resulting mean, standard deviation, and variance of this signal mixture?\n", + "- State the unbiased estimates of the mean, the standard deviation, and the variance for `xn`. Do they fulfill our expectations?" ] }, { @@ -1075,7 +1075,7 @@ "source": [ "# Task 12: STFT Magnitude Spectrum\n", "\n", - "For the created signal mixture `xn` calculate the STFT matrix $\\mathbf{X}_{\\mu,\\theta}$ for\n", + "For the created signal mixture `xn`, calculate the STFT matrix $\\mathbf{X}_{\\mu,\\theta}$ for\n", "\n", "- DFT length $N=64$ samples\n", "- hop size $H=32$ samples\n", @@ -1086,13 +1086,13 @@ "\n", "- Store the STFT matrix in variable `XN`.\n", "\n", - "For **our** implementation, we normalize by the squared peak power gain of the utilized window $w$. Furthermore we calculate the magnitude of the complex data.\n", + "For **our** implementation, we normalize by the squared peak power gain of the utilized window $w$. Furthermore, we calculate the magnitude of the complex data.\n", "This reads\n", "\\begin{align}\n", "\\frac{|\\mathbf{X}_{\\mu,\\theta}|}{\\sqrt{\\left(\\sum_{k=0}^{N-1} w[k]\\right)^2}}\n", "\\end{align}\n", "\n", - "which should be stored into `XNMag`.\n", + "which should be stored in `XNMag`.\n", "Note that $\\left(\\sum_{k=0}^{N-1} w[k]\\right)^2 = W^2(\\mathrm{e}^{\\mathrm{j}\\Omega=0})$ is the so called peak power gain, cf. equation (10) in [Har78].\n", "The maximum value of 20 lg(`XNMag`) should be very close to 4.74 dB. The exact value depends on the random noise `noise`. \n", "Furthermore, we might validate that max(abs(`XNMag`))$\\approx\\sqrt{2}$ when the amplitude of the **noise signal** is set to **zero**. \n", @@ -1188,13 +1188,13 @@ "\n", "The surface plot might be realized with Matlab's `surf()` with `FaceColor='flat'` and `EdgeColor='none'` or \n", "Matplotlib's / PyPlot's `pcolormesh()` with `shading='flat'` and `edgecolors='face'`. \n", - "Note, that for plotting no further manual data post processing is required.\n", + "Note that for plotting, no further manual data post-processing is required.\n", "Only the correct time and frequency vectors must be set up for **flat shading**.\n", - "This requires that the plot starts at time and frequency values of zero (see the pyplot/Matlab help files how flat shading aligns the coordinates with the color patches). \n", + "This requires that the plot start at time and frequency values of zero (see the pyplot/Matlab help files on how flat shading aligns the coordinates with the color patches). \n", "\n", "Pay special attention that the colorbar has a discretized colormap.\n", - "Have a look to a matplotlib [surface plot example](https://github.com/spatialaudio/digital-signal-processing-exercises/blob/master/jupyter/jupyter_intro.ipynb) on how to do this.\n", - "Suitable colormaps are Matlab's new standard `parula` and matplotlib's `viridis`, `cividis`, `inferno`, `plasma` and `magma`. The figures below use the `inferno` colormap (weird naming?!). Please do not use `jet` (Matlab's earlier default) and `hsv`, since these colormaps do not consider that we have different perceptual sensitivity for different colors.\n", + "Have a look at a matplotlib [surface plot example](https://github.com/spatialaudio/digital-signal-processing-exercises/blob/master/jupyter/jupyter_intro.ipynb) on how to do this.\n", + "Suitable colormaps are Matlab's new standard `parula` and matplotlib's `viridis`, `cividis`, `inferno`, `plasma`, and `magma`. The figures below use the `inferno` colormap (weird naming?!). Please do not use `jet` (Matlab's earlier default) and `hsv`, since these colormaps do not consider that we have different perceptual sensitivity for different colors.\n", "\n", "- Interpret these plots." ] @@ -1274,9 +1274,9 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "By now, we got an idea, how to approach an analysis of time-frequency varying signal characteristics with the above surface plot.\n", + "By now, we have an idea of how to approach an analysis of time-frequency varying signal characteristics with the above surface plot.\n", "\n", - "In other analysis tasks we are interested in the power spectral density of a signal. We will approach this with the third task section in the following." + "In other analysis tasks, we are interested in the power spectral density of a signal. We will approach this with the third task section in the following." ] }, { @@ -1287,7 +1287,7 @@ "\n", "## Power Spectral Density Estimator\n", "\n", - "From the DSP lecture we know that the power spectral density (PSD) $\\hat{\\Phi}_{xx}(\\Omega)$ can be estimated by \n", + "From the DSP lecture, we know that the power spectral density (PSD) $\\hat{\\Phi}_{xx}(\\Omega)$ can be estimated by \n", "\n", "\\begin{equation}\n", "\\hat{\\Phi}_{xx}(\\Omega) = \\frac{1}{N \\cdot U} |X(\\Omega)|^2,\n", @@ -1314,8 +1314,8 @@ "\\hat{\\Phi}_{xx}(\\Omega) = \\frac{|X(\\Omega)|^2}{f_s \\cdot \\sum_{k=0}^{N-1} w^2[k]}.\n", "\\end{equation}\n", "\n", - "It has physical dimension of V$^2$ / Hz if $x(t)$ has dimension of V for physical sampling frequency in Hz.\n", - "If dealing in digital domain, the unit of the PSD and the averaged PSD is (signal amplitude)$^2$ / (rad/s)." + "It has a physical dimension of V$^2$ / Hz if $x(t)$ has a dimension of V for a physical sampling frequency in Hz.\n", + "If dealing in the digital domain, the unit of the PSD and the averaged PSD is (signal amplitude)$^2$ / (rad/s)." ] }, { @@ -1324,7 +1324,7 @@ "source": [ "## Periodogram\n", "\n", - "If we derive the PSD estimator from DFT data and for a rectangular window, the so called periodogram is obtained. For very large $N$ this estimator is unbiased. However, it is **not** a **[consistent](https://en.wikipedia.org/wiki/Consistent_estimator)** estimate (cf. the [lecture](https://github.com/spatialaudio/digital-signal-processing-lecture/blob/master/spectral_estimation_random_signals/introduction.ipynb))." + "If we derive the PSD estimator from DFT data and for a rectangular window, the so called periodogram is obtained. For very large $N$, this estimator is unbiased. However, it is **not** a **[consistent](https://en.wikipedia.org/wiki/Consistent_estimator)** estimate (cf. the [lecture](https://github.com/spatialaudio/digital-signal-processing-lecture/blob/master/spectral_estimation_random_signals/introduction.ipynb))." ] }, { @@ -1336,8 +1336,8 @@ "The periodogram exhibits more or less (but never zero) deviations from the true PSD.\n", "A naive approach would be taking an average of multiple periodograms in order to check what signal characteristics 'remains steady'.\n", "That's what we would intuitively do in practice, right?\n", - "Of course there is a strong theoretical foundation and statistical proof for this approach, see e.g. chapter 10.5 in [OS10].\n", - "By taking the mean over $K$ periodograms\n", + "Of course, there is a strong theoretical foundation and statistical proof for this approach, see e.g., chapter 10.5 in [OS10].\n", + "By taking the mean over the $K$ periodograms\n", "\n", "\\begin{equation}\n", "\\frac{1}{K}\\sum_K \\hat{\\Phi}_{xx}(\\Omega) = \\frac{1}{K}\\sum_K \\frac{|X(\\Omega)|^2}{f_s \\cdot \\sum_{k=0}^{N-1} w^2[k]}.\n", @@ -1364,7 +1364,7 @@ "source": [ "# Task 14: Welch Periodogram\n", "\n", - "Here the task is to calculate the **Welch** peridogram for the STFT data `XN` created in task 12. So step by step we proceed as follows: \n", + "Here, the task is to calculate the **Welch** periodogram for the STFT data `XN` created in task 12. So step by step, we proceed as follows: \n", "\n", "1. We use the same parameters as for the surface plot\n", "\n", @@ -1372,7 +1372,7 @@ " - hop size $H=32$ samples\n", " - normal, periodic (non-symmetric) Hamming window\n", "\n", - "to create STFT matrix `XN` of the input signal `xn`. \n", + "to create the STFT matrix `XN` of the input signal `xn`. \n", "\n", "2. The above discussed normalization to `XN` needs to be performed.\n", "\n", @@ -1380,7 +1380,7 @@ "\n", "We assume a sampling frequency of $2 \\pi$. This will give us the PSD estimate over **digital angular frequency** $\\Omega$.\n", "Of course, we have only data at discrete DFT frequencies, but this is ok for our purpose.\n", - "Though, we could interpolate DTFT data from DFT data, see [exercise](https://github.com/spatialaudio/digital-signal-processing-exercises/blob/master/dft/dft_to_dtft_interpolation.ipynb).\n", + "Though we could interpolate DTFT data from DFT data, see [exercise](https://github.com/spatialaudio/digital-signal-processing-exercises/blob/master/dft/dft_to_dtft_interpolation.ipynb).\n", "\n", "4. Store the Welch periodogram in variable `PxxWelch` and the angular frequencies of the DFT in `Om`." ] @@ -1432,7 +1432,7 @@ "# Task 15: Plot the Welch Periodogram\n", "\n", "See the figure below showing an exemplary Welch periodogram of `xn`.\n", - "This might slightly vary for different runtimes, when using a random noise signal.\n", + "This might slightly vary for different runtimes when using a random noise signal.\n", "As validation we might check that the maximum in this plot is about 2.5868, when **noise signal** amplitude is set to **zero**. By the way: is the plot for this case as we expect it?\n", "\n", "1. Generate the plot.\n", @@ -1506,15 +1506,15 @@ "\n", "To summarize:\n", "\n", - "- We have implemented the basic blocks for STFT-based analysis and synthesis using different overlap-add methods. It is comparably simple to add the further steps such that we can apply time-frequency processing for e.g. filtering.\n", + "- We have implemented the basic blocks for STFT-based analysis and synthesis using different overlap-add methods. It is comparably simple to add further steps so that we can apply time-frequency processing, for e.g., filtering.\n", "\n", - "- Then we checked in detail how a surface plot of STFT data is created. These insights and playing around with the parameters help to understand what is going on. We should get experience on our own regarding different time and different frequency resolution, zeropadding and impact of different windows. Once we set the basics here, these steps can be conveniently performed and better understood.\n", + "- Then we checked in detail how a surface plot of STFT data is created. These insights and playing around with the parameters help to understand what is going on. We should get experience on our own regarding different time and different frequency resolution, zero-padding, and the impact of different windows. Once we set the basics here, these steps can be conveniently performed and better understood.\n", "\n", - "- At last, we got into detail of the modified periodogram, more specifically the Welch periodogram. This is (besides parametric approaches, see the lecture) one of fundamental 'need to know' in DSP, when dealing with random data, which in practice is most often the case. Here, we also should play around with different averages, different DFT block sizes, different windows to get experience on the made up signal...and then trying other signals.\n", + "- At last, we got into the details of the modified periodogram, more specifically, the Welch periodogram. This is (besides parametric approaches, see the lecture) one of the fundamental 'need to know' in DSP, when dealing with random data, which in practice is most often the case. Here, we also should play around with different averages, different DFT block sizes, different windows to get experience on the made-up signal...and then try other signals.\n", "\n", - "Last reminder: After implementing stuff on our own, we might better understand the theoretical foundations and the equations of the discussed approaches. So, check back the lectures and text books, it will be enlightening.\n", + "Last reminder: After implementing stuff on our own, we might better understand the theoretical foundations and the equations of the discussed approaches. So, check back the lectures and textbooks, it will be enlightening.\n", "\n", - "Furthermore, we experienced that tiny details can make big differences. So, in-depth-checking things/tools, that we employ, is always a good idea when engineering." + "Furthermore, we experienced that tiny details can make big differences. So, in-depth checking of things/tools that we employ is always a good idea when engineering." ] }, { @@ -1583,7 +1583,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.10.12" + "version": "3.12.2" } }, "nbformat": 4,