# Exercise 2: Sinusoids and the DFT

Doing this exercise you will get a better understanding of the basics elements and operations that take place in the Discrete Fourier Transform (DFT). There are five parts: 1) Generate a sinusoid, 2) Generate a complex sinusoid, 3) Implement the DFT, 4) Implement the IDFT, and 5) Compute the magnitude spectrum of an input sequence.

### Relevant Concepts

A real sinusoid in discrete time domain can be expressed by:

\begin{equation}
x[n] = A\cos(2 \pi fnT + \varphi)
\end{equation}

where, $x$ is the array of real values of the sinusoid, $n$ is an integer value expressing the time index, $A$ is the amplitude value of the sinusoid, $f$ is the frequency value of the sinusoid in Hz, $T$ is the sampling period equal to $1/fs$, fs is the sampling frequency in Hz, and $\varphi$ is the initial phase of the sinusoid in radians.

A complex sinusoid in discrete time domain can be expressed by:

\begin{equation}
\bar{x}[n] = Ae^{j(\omega nT + \varphi)} = A\cos(\omega nT + \varphi)+ j A\sin(\omega nT + \varphi)
\end{equation}

where, $\bar{x}$ is the array of complex values of the sinusoid, $n$ is an integer value expressing the time index, $A$ is the amplitude value of the sinusoid, $e$ is the complex exponential number, $\omega$ is the frequency of the sinusoid in radians per second (equal to $2 \pi f$), $T$ is the sampling period equal $1/fs$, fs is the sampling frequency in Hz and $\varphi$ is the initial phase of the sinusoid in radians.

The $N$ point DFT of a sequence of real values $x$ (a sound) can be expressed by:

\begin{equation}
X[k] = \sum_{n=0}^{N-1} x[n]e^{-j2 \pi kn/N}  \hspace{1cm} k=0,...,N-1
\end{equation}

where $n$ is an integer value expressing the discrete time index, $k$ is an integer value expressing the discrete frequency index, and $N$ is the length of the DFT.

The IDFT of a spectrum $X$ of length $N$ can be expressed by:

\begin{equation}
x[n] = \frac{1}{N} \sum_{k=0}^{N-1} X[k]e^{j2 \pi kn/N}  \hspace{1cm} n=0,...,N-1
\end{equation}

where, $n$ is an integer value expressing the discrete time index, $k$ is an integer value expressing the discrete frequency index, and $N$ is the length of the spectrum $X$.

The magnitude of a complex spectrum $X$ is obtained by taking its absolute value: $|X[k]| $


## Part 1 - Generate a sinusoid

Complete the function `gen_sine` to generate a real sinusoid (use `np.cos()`) given its amplitude `A`, frequency `f` (Hz), initial phase `phi` (radians), sampling rate `fs` (Hz) and duration `t` (seconds). 

All the input arguments to this function `(A, f, phi, fs and t)` are real numbers such that `A`, `t` and `fs` are positive, and `fs > 2*f` to avoid aliasing. The function should return a numpy array `x` of the generated sinusoid. 

Use the function `cos` of the numpy package to compute the sinusoidal values.


```python
import numpy as np
```


```python
def gen_sine(A, f, phi, fs, t):
    """
    Inputs:
        A (float) =  amplitude of the sinusoid
        f (float) = frequency of the sinusoid in Hz
        phi (float) = initial phase of the sinusoid in radians
        fs (float) = sampling frequency of the sinusoid in Hz
        t (float) =  duration of the sinusoid (is second)
    Output:
        x (numpy array) = generated sinusoid
    """
    ### your code here

    
```

Now call and test the `gen_sine()` function. If you use `A=1.0, f = 10.0, phi = 1.0, fs = 50` and `t = 0.1`, the output numpy array should be:

```
array([ 0.54030231, -0.63332387, -0.93171798,  0.05749049,  0.96724906])
```

To generate a sinewave that you can hear, it should be longer and with a higher sampling rate. For example you can use `A=1.0, f = 440.0, phi = 1.0, fs = 5000` and `t = 0.5`. To play it import the `Ipython.display` package and use `ipd.display(ipd.Audio(data=x, rate=fs))`.


```python
import IPython.display as ipd

### your code here

```

## Part 2 - Generate a complex sinusoid 

Complete the `gen_complex_sine()` function to generate the complex sinusoid that is used in DFT computation of length `N` (samples), corresponding to the frequency index `k`. Note that the complex sinusoid used in DFT computation has a negative sign in the exponential function.

The amplitude of such a complex sinusoid is `1`, the length is `N`, and the frequency in radians is `2*pi*k/N`.

The input arguments to the function are two positive integers, `k` and `N`, such that `k < N-1`. The function should return `c_sine`, a numpy array of the complex sinusoid. Use the function `exp()` of the numpy package to compute the complex sinusoidal values.


```python
def gen_complex_sine(k, N):
    """
    Inputs:
        k (integer) = frequency index of the complex sinusoid of the DFT
        N (integer) = length of complex sinusoid in samples
    Output:
        c_sine (numpy array) = generated complex sinusoid (length N)
    """
    ### your code here

```

Now run an test the `gen_complex_sine()` function. If you run it using `k=1` and `N=5`, the function should return the following numpy array:

```
array([ 1.0 + 0.j,  0.30901699 - 0.95105652j, -0.80901699 - 0.58778525j, -0.80901699 + 0.58778525j, 0.30901699 + 0.95105652j])
```


```python
# call gen_complex_sine
### your code here

```

## Part 3 - Implement the discrete Fourier transform (DFT)

Complete the function `dft()` to implement the discrete Fourier transform (DFT) equation given above. Given a sequence `x` of length `N`, the function should return its spectrum of length `N` with the frequency indexes ranging from 0 to `N-1`.

The input argument to the function is a numpy array `x` and the function should return a numpy array `X`, the DFT of `x`.


```python
def dft(x):
    """
    Input:
        x (numpy array) = input sequence of length N
    Output:
        X (numpy array) = N point DFT of the input sequence x
    """
    ## Your code here

```

Now run and test the function `dft()`. If you run it using as input `x = np.array([1, 2, 3, 4])`, the function shoulds return the following numpy array:

```
array([10.0 + 0.0j,  -2. +2.0j,  -2.0 - 9.79717439e-16j, -2.0 - 2.0j])
```

Note that you might not get an exact 0 in the output because of the small numerical errors due to the limited precision of the data in your computer. Usually these errors are of the order 1e-15 depending on your machine.


```python
# call dft
### your code here

```

## Part 4 - Implement the inverse discrete Fourier transform (IDFT)

Complete the function `idft()` to implement the inverse discrete Fourier transform (IDFT) equation given above. Given a frequency spectrum `X` of length `N`, the function should return its IDFT `x`, also of length `N`. Assume that the frequency index of the input spectrum ranges from 0 to `N-1`.

The input argument to the function is a numpy array `X` of the frequency spectrum and the function should return a numpy array of the IDFT of `X`.

Remember to scale the output appropriately.


```python
def idft(X):
    """
    Input:
        X (numpy array) = frequency spectrum (length N)
    Output:
        x (numpy array) = N point IDFT of the frequency spectrum X
    """
    ### Your code here

```

Now run and test the `idft()` function. If you run it with the input `X = np.array([1 ,1 ,1 ,1])`, the function should return the following numpy array: 

```
array([  1.00000000e+00 +0.00000000e+00j,   -4.59242550e-17 +5.55111512e-17j,   0.00000000e+00 +6.12323400e-17j,   8.22616137e-17 +8.32667268e-17j])
```

Notice that the output numpy array is essentially `[1, 0, 0, 0]`. Instead of exact 0 we get very small numerical values of the order of 1e-15, which can be ignored. Also, these small numerical errors are machine dependent and might be different in your case.

In addition, an interesting test of the IDFT function can be done by providing the output of the DFT of a sequence as the input to the IDFT. See if you get back the original time domain sequence.


```python
# call idft
### your code here

```

## Part 5 - Compute the magnitude spectrum

Complete the function `gen_mag_spectrum()` to compute the magnitude spectrum of an input sequence `x` of length `N`. The function should return an `N` point magnitude spectrum with frequency index ranging from 0 to `N-1`.

The input argument to the function is a numpy array `x` and the function should return a numpy array of the magnitude spectrum of `x`.


```python
def gen_mag_spec(x):
    """
    Input:
        x (numpy array) = input sequence of length N
    Output:
        magX (numpy array) = magnitude spectrum of the input sequence x (length N)
    """
    ### Your code here

```

Now run and test the function `gen_mag_spec()`. If you run `gen_mag_spec()` using as input `x = np.array([1, 2, 3, 4])`, it should return the following numpy array:
```
array([10.0, 2.82842712, 2.0, 2.82842712])
```
For a more realistic use of `gen_mag_spec()` use as input a longer signal, such as `x = np.cos(2*np.pi*200.0*np.arange(512)/1000)`, and to get a visual representation of the input and output, import the `matplotlib.pyplot` package and use `plt.plot(x)` and `plt.plot(abs(X))`.


```python
import IPython.display as ipd
import matplotlib.pyplot as plt
```


```python
# call gen_mag_spec and plot result
### your code here


```

7a5bd2850ef299f1eb2c7ff4f878140aac8e4984
ipynb
Jupyter Notebook
stanford/sms-tools/notebooks/E2-Sinusoids-and-DFT.ipynb
phunc20/dsp
e7c496eb5fd4b8694eab0fc049cf98a5e3dfd886
MIT
2021-03-12T18:32:06.000Z
2021-03-12T18:32:06.000Z
stanford/sms-tools/notebooks/E2-Sinusoids-and-DFT.ipynb
phunc20/dsp
e7c496eb5fd4b8694eab0fc049cf98a5e3dfd886
MIT
stanford/sms-tools/notebooks/E2-Sinusoids-and-DFT.ipynb
phunc20/dsp
e7c496eb5fd4b8694eab0fc049cf98a5e3dfd886
MIT
Qwen/Qwen-72B
1. YES
2. YES