<jupyter_start><jupyter_text>Exercise 3: Fourier propertiesWith this exercise you will get a better understanding of some of the Fourier theorems and of some useful properties of the DFT. You will write code to implement and verify several properties of the DFT that are discussed in the lectures. You will also learn to use the `dftModel.py` module of sms-tools, which contains the basic python functions implementing the DFT. There are five parts in the exercise: 1) Minimize energy spread in DFT of sinusoids, 2) Optimal zero-padding, 3) Symmetry properties of the DFT, 4) Suppressing frequency components using DFT model, and 5) FFT size and zero-padding.  Relevant Concepts__DFT of sinusoids:__ When a real sinusoid has an integer number of cycles in $N$ samples, the frequency of the sinusoid exactly matches one of the bin frequencies in an $N$ point DFT. Hence the DFT spectrum of the sinusoid has a value of zero at every DFT bin except at the two bins that match the frequency of the sinusoid. Otherwise, the energy of the sinusoid is spread over all the bins. When there are multiple sinusoids, the equations extend to each sinusoid. \begin{eqnarray} x[n]&=&A_{0}\cos\left(2\pi k_{0}n/N\right)=\frac{A_{0}}{2}{\textstyle e}^{j2\pi k_{0}n/N}+\frac{A_{0}}{2}{\textstyle e}^{-j2\pi k_{0}n/N}\\ X[k] &=& \frac{A_0}{2} \,\,\, \mathrm{for} \,\,\, k = k_0, -k_0; \,\,\,\, 0 \,\,\,\mathrm{otherwise}\end{eqnarray}__Zero-padding:__ Zero-padding a signal is done by adding zeros at the end of the signal. If we perform zero-padding to a signal before computing its DFT, the resulting spectrum will be an interpolated version of the spectrum of the original signal.  In most implementations of the DFT (including the FFT algorithms) when the DFT size is larger than the length of the signal, zero-padding is implicitly done. __Zero phase windowing:__ Zero phase windowing of a frame of a signal puts the centre of the signal at the zero time index for DFT computation. By moving the centre of the frame to zero index by a circular shift, the computed DFT will not have the phase offset which would have otherwise been introduced (recall that a shift of the signal causes the DFT to be multiplied by a complex exponential, which keeps the magnitude spectrum intact but changes the phase spectrum). When used in conjunction with zero-padding, zero phase windowing is also useful for the creation of a frame of length of power of 2 for FFT computation (`fftbuffer`). If the length of the signal $x$ is $M$ and the required DFT size is $N$, the zero phase windowed version of the signal, `dftbuffer`, for DFT computation can be obtained by (works for both even and odd $M$):     hM1 = floor((M+1)/2)    hM2 = floor(M/2)    dftbuffer = zeros(N)    dftbuffer[:hM1] = x[hM2:]                                  dftbuffer[-hM2:] = x[:hM2]__Real, even and odd signals:__ A signal is real when it does not have any imaginary component, and all sounds are real signals. A signal $x$ is even if $x[n] = x[-n]$, and odd if $x[n] = -x[-n]$. For a signal of length $M$ (and $M$ is odd), in the context of a zero phase windowed signal and its DFT, the signal is even if $x[n] = x[M-n]$ and odd if $x[n] = -x[M-n]$, $1 \leq n \leq M-1$. The DFT properties show that for real input signals, the magnitude spectrum is even and the phase spectrum is odd. Furthermore, when the input signal is both real and even, the DFT is real valued, with an even magnitude spectrum and imaginary component equal to zero. In summary, if $x$ is an input signal of length $M$ ($M$ is odd) and $X = \mathrm{DFT}(x,M)$, then for $1 \leq k \leq M-1$If $x$ is real, $\left|X[k]\right| = \left|X[M-k]\right|$ and $\boldsymbol{<}\!X[k] = -\boldsymbol{<}\!X[M-k]$If $x$ is real and even, $\left|X[k]\right| = \left|X[M-k]\right|$ and $\mathrm{imag}(X[k]) = 0$__Positive half of the DFT spectrum:__ Audio signals are real signals. Due to the symmetry properties of the DFT of a real signal, it is sufficient to store only one half of the magnitude and phase spectra. To save on both storage and computation, we will just store just the half spectrum when possible. From an $N$ point DFT ($N$ even), we can obtain the positive half of the spectrum by considering only the first $(N/2)+1$ samples of the DFT. We can compute the magnitude spectrum of the positive half (in dB) as $m_X = 20\log_{10}\left|X[0:(N/2)+1]\right|$, where $X$ is the DFT of the input. __Filtering:__ Filtering involves selectively suppressing certain frequencies present in the signal. Filtering is often performed in the time domain by the convolution of the input signal with the impulse response of a filter. The same operation can also be done in the DFT domain using the properties of DFT, by multiplying the DFT of the input signal by the DFT of the impulse response of the filter. In this assignment, we will consider a very simple illustrative filter that suppresses some frequency components by setting some DFT coefficients to zero. It is to be noted that the convolution operation here is circular convolution with a period $N$, the size of the DFT. If $x_1[n] \Leftrightarrow X_1[k]$ and $x_2[n] \Leftrightarrow X_2[k]$, $x_1[n] * x_2[n] \Longleftrightarrow X_1[k]\,X_2[k]$  Part 1 - Minimize energy spread in DFT of sinusoidsGiven an input signal consisting of two sinusoids, complete the function `minimize_energy_spread_dft()` to select the first `M` samples from the signal and return the positive half of the DFT magnitude spectrum (in dB), such that it has only two non-zero values. `M` is to be calculated as the smallest positive integer for which the positive half of the DFT magnitude spectrum has only two non-zero values. To get the positive half of the spectrum, first compute the `M` point DFT of the input signal (for this you can use the `fft()` function of `scipy.fftpack`). Consider only the first `(M/2)+1` samples of the DFT, computing the magnitude spectrum of the positive half (in dB) as `mX = 20*log10(abs(X[:M/2+1]))`, where `X` is the DFT of the input signal.The input arguments to this function are the input signal `x` (of length W >= M) consisting of two sinusoids of frequency `f1` and `f2`, the sampling frequency `fs` and the value of frequencies `f1` and `f2`. The function should return the positive half of the magnitude spectrum `mX`. For this question, you can assume the input frequencies `f1` and `f2` to be positive integers and factors of `fs`, and that `M` is even. Due to the precision of the FFT computation, the zero values of the DFT are not zero but very small values < 1e-12 (or -240 dB) in magnitude. For practical purposes, all values with absolute value less than 1e-6 (or -120 dB) can be considered to be zero. HINT: The DFT magnitude spectrum of a sinusoid has only one non-zero value (in the positive half of the DFT spectrum) when its frequency coincides with one of the DFT bin frequencies. This happens when the DFT size (`M` in this question) contains exactly an integer number of periods of the sinusoid. Since the signal in this question consists of two sinusoids, this condition should hold true for each of the sinusoids, so that the DFT magnitude spectrum has only two non-zero values, one per sinusoid. `M` can be computed as the Least Common Multiple (LCM) of the sinusoid periods (in samples). The LCM of two numbers `x`, `y` can be computed as: `x*y/gcd(x,y)`, where gcd denotes the greatest common divisor.<jupyter_code>from scipy.fftpack import fft, fftshift
import numpy as np
from math import gcd, ceil, floor
import sys
sys.path.append('../software/models/')
from dftModel import dftAnal, dftSynth
from scipy.signal import get_window
import matplotlib.pyplot as plt
def minimize_energy_spread_dft(x, fs, f1, f2):
    """
    Inputs:
        x (numpy array) = input signal 
        fs (float) = sampling frequency in Hz
        f1 (float) = frequency of first sinusoid component in Hz
        f2 (float) = frequency of second sinusoid component in Hz
    Output:
        positive half of DFT spectrum (in dB) (numpy array) 
    """
    ## Your code here<jupyter_output><empty_output><jupyter_text>Run and test the function `minimize_energy_spread_dft()`. Plot the output spectrum to validate the result. Here are two possible examples to use._Test case 1:_ For an input signal `x` sampled at `fs = 10000`Hz that consists of sinusoids of frequencies `f1 = 80`Hz and `f2 = 200`Hz, you need to select `M = 250` samples of the signal to meet the required condition. In this case, output `mX` is 126 samples in length and has non-zero values at bin indices 2 and 5 (corresponding to the frequency values of 80 and 200 Hz, respectively). You can create a test signal `x` by generating and adding two sinusoids of the given frequencies._Test case 2:_ For an input signal `x` sampled at `fs = 48000` Hz that consists of sinusoids of frequencies `f1 = 300`Hz and `f2 = 800`Hz, you need to select `M = 480` samples of the signal to meet the required condition. In this case, output `mX` is 241 samples in length and has non-zero values at bin indices 3 and 8 (corresponding to the frequency values of 300 and 800 Hz, respectively). You can create a test signal `x` by generating and adding two sinusoids of the given frequencies.<jupyter_code>## Your code here<jupyter_output><empty_output><jupyter_text>Part 2 - Optimal zero-paddingGiven a sinusoid, complete the function `optimal_zeropad()` to compute the DFT of the sinusoid after zero-padding and return the positive half of the magnitude spectrum (in dB). Zero-padding needs to be done such that one of the bin frequencies of the DFT coincides with the frequency of the sinusoid. Choose the minimum zero-padding length for which this condition is satisfied. The input arguments are the sinusoid `x` of length `M`, sampling frequency `fs` and the frequency of the sinusoid `f`. The output is the positive half of the magnitude spectrum `mX` computed using the `N` point DFT `(N >= M)` of `x` after zero-padding `x` to length `N` appropriately as required. To get the positive half of the spectrum, first compute the `N` point DFT of the zero-padded input signal (for this you can use the `fft()` function of `scipy.fftpack`, which is already imported in this script). Consider only the first `(N/2)+1` samples of the DFT and compute the magnitude spectrum of the positive half (in dB) as `mX = 20*log10(abs(X[:(N/2)+1]))`, where `X` is the `N` point DFT of the zero-padded input.For this exercise, you can assume that the frequency of the sinusoid `f` is a positive integer and a factor of the sampling rate `fs`. The input parameters will be given in such a way that `N` will be even.  Note that the numerical value of `f` is an integer but the data type is float, for example 1.0, 2.0, 55.0 etc. This is to avoid issues in python related with division by a integer.Due to the precision of the FFT computation, the zero values of the DFT are not zero but very small values < 1e-12 (or -240 dB) in magnitude. For practical purposes, all values with absolute value less than 1e-6 (or -120 dB) can be considered to be zero. HINT: One of the DFT bin frequencies coincides with the frequency `f` of a sinusoid when the DFT size (`N` in this question) contains exactly an integer number of periods of the sinusoid. For example, if `f = 100` Hz and `fs = 1000` Hz, one period of the sinusoid has 10 samples. Then given a signal of length `M = 25` samples, there are 2.5 periods in it. The minimum zero-padding length here would be 5 samples (0.5 period), so that the DFT size `N = 30` corresponds to 3 periods of a sinusoid of frequency `f=100` Hz.<jupyter_code>def optimal_zeropad(x, fs, f):
    """
    Inputs:
        x (numpy array) = input signal of length M
        fs (float) = sampling frequency in Hz
        f (float) = frequency of sinusoid in Hz
    Output:
        positive half of DFT spectrum of N point DFT after zero-padding (numpy array)
    """
    ## Your code here
    
    T = fs/f
    M = len(x)
    N = int(T*ceil(M/T))
    xw = np.concatenate((x,np.zeros(N-M)),axis = 0)
    X = fft(x)
    return 20*np.log10(abs(X[:int(N/2+1)]))<jupyter_output><empty_output><jupyter_text>Run and test the function `optimal_zeropad()`, validating the result by ploting the output array. Here are two possible examples._Test case 1:_ For a sinusoid `x` with `f = 100.0` Hz, `M = 25` samples and `fs = 1000.0` Hz, you will need to zero-pad by 5 samples and compute an `N = 30` point DFT. In the magnitude spectrum, you can see a maximum value at bin index 3 corresponding to the frequency of 100 Hz. The output array you return is 16 samples in length. _Test case 2:_ For a sinusoid `x` with `f = 250.0` Hz, `M = 210` samples and `fs = 10000.0` Hz, you will need to zero-pad by 30 samples and compute an `N = 240` point DFT. In the magnitude spectrum, you can see a maximum value at bin index 6 corresponding to the frequency of 250 Hz. The output array you return is 121 samples in length.<jupyter_code>## Your code here<jupyter_output><empty_output><jupyter_text>Part 3 - Symmetry properties of the DFTComplete the function `test_real_even()` a function to check if the input signal is real and even using the symmetry properties of its DFT. The function will return the result of this test, the zerophase windowed version of the input signal (`dftbuffer`), and its DFT. Given an input signal `x` of length `M`, do a zero phase windowing of `x` without any zero-padding. Then compute the `M` point DFT of the zero phase windowed signal and use the symmetry of the computed DFT to test if the input signal `x` is real and even. Return the result of the test, the `dftbuffer` computed, and the DFT of the `dftbuffer`. The input argument is a signal `x` of length `M`. The output is a tuple with three elements `(isRealEven, dftbuffer, X)`, where `isRealEven` is a boolean variable which is `True` if `x` is real and even, else `False`. `dftbuffer` is the `M` length zero phase windowed version of `x`. `X` is the `M` point DFT of the `dftbuffer`. To make the problem easier, we will use odd length input sequence in this question (`M` is odd). Due to the precision of the FFT computation, the zero values of the DFT are not zero but very small values < 1e-12 in magnitude. For practical purposes, all values with absolute value less than 1e-6 can be considered to be zero. Use an error tolerance of 1e-6 to compare if two floating point arrays are equal. Caveat: Use the imaginary part of the spectrum instead of the phase to check if the input signal is real and even.<jupyter_code>def test_real_even(x):
    """
    Inputs:
        x (numpy array)= input signal of length M (M is odd)
    Output:
        tuple including:
        isRealEven (boolean) = True if input x is real and even, and False otherwise
        dftbuffer (numpy array, possibly complex) = M point zero phase windowed version of x 
        X (numpy array, possibly complex) = M point DFT of dftbuffer 
    """
    ## Your code here<jupyter_output><empty_output><jupyter_text>Run and test the function `test_real_evenI ()`. Here are two possible examples._Test case 1:_ If `x = np.array([ 2, 3, 4, 3, 2 ])`, which is a real and even signal (after zero phase windowing), the function returns ```(True, array([ 4., 3., 2., 2., 3.]), array([14.0000+0.j, 2.6180+0.j, 0.3820+0.j, 0.3820+0.j, 2.6180+0.j])) (values are approximate)```_Test case 2:_ If `x = np.array([1, 2, 3, 4, 1, 2, 3])`, which is not an even signal (after zero phase windowing), the function returns ```(False,  array([ 4.,  1.,  2.,  3.,  1.,  2.,  3.]), array([ 16.+0.j, 2.+0.69j, 2.+3.51j, 2.-1.08j, 2.+1.08j, 2.-3.51j, 2.-0.69j])) (values are approximate)```To get a more realistic example use a longer input signal and plot the real and imaginary parts of the output spectrum `X`. For example, use `x = get_window('hanning', 51, fftbins=False)`, which is real an even, and plot `x`and the real and imaginary part of the spectrum `X`.<jupyter_code>## Your code here<jupyter_output><empty_output><jupyter_text>Part 4 - Suppressing frequency components using DFT modelGiven a signal as input, complete the function `supress_freq_dft_model()` to suppress the frequency components <= 70Hz using the DFT. Return the filtered signal in the time domain.Use the DFT to implement a very basic form of frequency domain filtering. Use the functions `dftAnal()` and `dftSynth()` provided in the `dftModel.py` module. Use `dftAnal()` to obtain the magnitude spectrum (in dB) and phase spectrum of the audio signal. Set the values of the magnitude spectrum that correspond to frequencies <= 70 Hz to -120dB (there may not be a bin corresponding exactly to 70Hz, choose the nearest bin of equal or higher frequency, e.g., using `np.ceil()`).Use `dftSynth()` to synthesize the filtered output signal. Then return the filtered signal. Use a hamming window to smooth the signal. Hence, do not forget to scale the output signals by the sum of the window values (as done in `software/models_interface/dftModel_function.py`).  Please note that this question is just for illustrative purposes and filtering is not usually done this way - such sharp cutoffs introduce artifacts in the output. The input is a `M` length signal `x`, sampling frequency is `fs` and the FFT size `N`. The output is the filtered signal.<jupyter_code>def suppress_freq_dft_model(x, fs, N):
    """
    Inputs:
        x (numpy array) = input signal of length N (odd)
        fs (float) = sampling frequency (Hz)
        N (positive integer) = FFT size
    Outputs:
       output signal with filtering (N samples long) (numpy array) 
    """
    N = len(x)
    w = get_window('hamming', N)
    outputScaleFactor = sum(w)
    
    ## Your code here<jupyter_output><empty_output><jupyter_text>Run and test the function `suppress_freq_dft_model()`. Plot the spectrum of the output over the spectrum of the input signals to validate the result. Here are two possible examples to try._Test case 1:_ For an input signal with 40Hz, 100Hz, 200Hz, 1000Hz components, the output should only contain 100Hz, 200Hz and 1000Hz components. _Test case 2:_ For an input signal with 23Hz, 36Hz, 230Hz, 900Hz, 2300Hz components, the output should only contain 230Hz, 900Hz and 2300Hz components. To understand the effect of filtering, you can plot the magnitude spectra of the input and output signals superposed.<jupyter_code>## Your code here<jupyter_output><empty_output><jupyter_text>Part 5 - Window-size, FFT-size and zero-paddingComplete the function `zp_fft_size_expt()`to take an input signal, compute three different magnitude spectra (with different parameters) and return them. This function should provide some insights into the effects window-size, FFT-size, and zero-padding on the spectrum of a signal. Te input signal should be of size 512 samples, the sampling rate should be 1000Hz, and the analysis window used should be hamming. The three set of analysis parameters should be:1. window-size = 256, FFT-size = 256 (no zero-padding)2. window-size = 512, FFT-size = 512 (no zero-padding)2. window-size = 256, FFT-size = 512 (zero-padding of 256 samples)Use `dftAnal()` to obtain the positive half of the magnitude spectrum (in dB). Return the 3 magnitude spectra in dB.<jupyter_code>def zp_fft_size_expt(x):
    """
    Inputs:
        x (numpy array) = input signal (512 samples long)
    Output:
        magnitude spectrum with window-size = 256, FFT-size = 256 (numpy array) 
        magnitude spectrum with window-size = 512, FFT-size = 512 (numpy array) 
        magnitude spectrum with window-size = 256, FFT-size = 512 (numpy array) 
    """
    
    N = len(x)    # should be 512 
    w1 = get_window('hamming', N)
    w2 = get_window('hamming', N//2)
    
    ## Your code here<jupyter_output><empty_output><jupyter_text>Run and test the function `zp_fft_size_expt()`. Plot the outputs to understand and validate the results. Here is one possible example._Test case 1:_ Use as input `x = .2*np.cos(2*np.pi*300*n)+.2*np.cos(2*np.pi*600*n)` where `n=np.arange(512)/fs` and the sampling rate `fs=1000`. Call the function with `mX1, mX2, mX3 = zp_fft_size_expt(x)`To understand better, plot the output of `dftAnal()` for each case on a common frequency axis with different colors. You will see that `mX3` is the interpolated version of `mX1` (zero-padding leads to interpolation of the DFT). You will also observe that the 'mainlobe' of the magnitude spectrum in `mX2` will be narrower than that in `mX1` and `mX3`. This shows that having a longer window leads to a narrower mainlobe with better frequency resolution and less spreading of the energy of the sinusoid. If we were to estimate the frequency of the sinusoid using its DFT, a first principles approach is to choose the frequency value of the bin corresponding to the maximum in the DFT magnitude spectrum. If you were to take this approach, which of the magnitude spectra will give you a better estimate of the frequency of the sinusoid? Comment and discuss.<jupyter_code>## Your code here<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.08065053022017214, 'fraction_numerical': 0.022778068381758523, 'mean_word_length': 4.114299610894942, 'pattern_counts': {'":': 0, '<': 48, '<?xml version=': 0, '>': 43, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '10', 'dir': 'jupyter-structured-clean-dedup', 'id': '93754', 'license': 'MIT', 'n_tokens_mistral': 6390, 'n_tokens_neox': 5788, 'n_words': 3215, 'path': 'stanford/sms-tools/notebooks/E3-Fourier-properties.ipynb', 'repo_name': 'phunc20/dsp', 'size': '21029'}