<jupyter_start><jupyter_text>Load all modules<jupyter_code>import numpy as np
import matplotlib.pyplot as plt
import scipy.linalg as la
from scipy.fft import fft as fft_sp
from scipy.fft import ifft as ifft_sp
from numpy.fft import fft as fft_np
from numpy.fft import ifft as ifft_np
# sparse complex matrices look nicer
np.set_printoptions(precision=3, suppress=True)<jupyter_output><empty_output><jupyter_text>DFT Matrix<jupyter_code>def dft_matrix(N):
    k = np.arange(N)
    return np.exp(-1j * 2 * np.pi / N * np.outer(k, k))<jupyter_output><empty_output><jupyter_text>Tasks  Some parameters setting<jupyter_code>N = 2 ** 2 # must be even, should be power of two
Nh = N // 2 # interger division
W = dft_matrix(N) # full DFT matrix for analysis of signal with length N
Wh = dft_matrix(Nh) # half size DFT matrix that actually needs to be computed
print(W)
print(Wh)<jupyter_output>[[ 1.+0.j  1.+0.j  1.+0.j  1.+0.j]
 [ 1.+0.j  0.-1.j -1.-0.j -0.+1.j]
 [ 1.+0.j -1.-0.j  1.+0.j -1.-0.j]
 [ 1.+0.j -0.+1.j -1.-0.j  0.-1.j]]
[[ 1.+0.j  1.+0.j]
 [ 1.+0.j -1.-0.j]]<jupyter_text>Task 1: Twiddle Factor Define,$${\rm\bf W}_N^k = e^{-j\frac{2\pi}{N} k}\qquad k=\{0,1,2,\cdots,N-1\}$$then,$${\rm\bf W}_{N/2}^k = e^{\left(-j\frac{2\pi}{N} k\right)\cdot 2} = {\rm\bf W}_N^{2k}$$N/2-DFT consequence when giving the input signal ${\rm\bf x}_k$$${\rm\bf X}_\mu = \sum_{k=0}^{N/2-1}{\rm\bf x}_k {\rm\bf W}_{N/2}^{k\cdot \mu} = \sum_{k=0}^{N/2-1}{\rm\bf x}_k {\rm\bf W}_N^{2k\cdot \mu}$$  Task 2: Helping Matrices Identity matrix $${\rm\bf I}_{N/2}$$The diagonal matrix $${\rm\bf D}_{N/2} = {\rm diag}(e^{-j\frac{2\pi}{N}{\rm\bf b}})$$ using $${\rm\bf b} = [0,1,2,\cdots\ \frac{N}{2}-1]^{\rm T}$$<jupyter_code>Ih = np.identity(Nh)
b = np.arange(Nh)
Dh = np.diag(np.exp(-1j * (2 * np.pi / N) * b))
print(Ih)
print(Dh)<jupyter_output>[[1. 0.]
 [0. 1.]]
[[1.+0.j 0.+0.j]
 [0.+0.j 0.-1.j]]<jupyter_text>Task 3: Output Combination Matrix $${\rm\bf A} = \left[\begin{matrix}{\rm\bf I}_{N/2} & +{\rm\bf D}_{N/2} \\{\rm\bf I}_{N/2} & -{\rm\bf D}_{N/2}\end{matrix}\right]$$<jupyter_code>A = np.vstack((np.hstack((Ih, Dh)), np.hstack((Ih, -Dh))))<jupyter_output><empty_output><jupyter_text>Task 4: Validation of A<jupyter_code>print(A)<jupyter_output>[[ 1.+0.j  0.+0.j  1.+0.j  0.+0.j]
 [ 0.+0.j  1.+0.j  0.+0.j  0.-1.j]
 [ 1.+0.j  0.+0.j -1.-0.j -0.-0.j]
 [ 0.+0.j  1.+0.j -0.-0.j -0.+1.j]]<jupyter_text>Task 5: Sparse N/2-DFT Matrix $${\rm\bf B} = \left[\begin{matrix}{\rm\bf W}_{N/2} & 0 \\0 & {\rm\bf W}_{N/2}\end{matrix}\right]_{N\times N}$$<jupyter_code>Zeroh = np.zeros((Nh, Nh), dtype=complex)
B = np.vstack((np.hstack((Wh, Zeroh)), np.hstack((Zeroh, Wh))))<jupyter_output><empty_output><jupyter_text>Another way to contruct matrix ${\rm\bf B}$ using indices:<jupyter_code>B = np.zeros((N, N), dtype=complex)
B[:Nh, :Nh] = Wh
B[Nh:N, Nh:N] = Wh<jupyter_output><empty_output><jupyter_text>Task 6: Validation of B<jupyter_code>print(B)
print(Wh)<jupyter_output>[[ 1.+0.j  1.+0.j  0.+0.j  0.+0.j]
 [ 1.+0.j -1.-0.j  0.+0.j  0.+0.j]
 [ 0.+0.j  0.+0.j  1.+0.j  1.+0.j]
 [ 0.+0.j  0.+0.j  1.+0.j -1.-0.j]]
[[ 1.+0.j  1.+0.j]
 [ 1.+0.j -1.-0.j]]<jupyter_text>Task 7: Matrix Product A B $${\rm\bf C} = {\rm\bf A} {\rm\bf B}$$<jupyter_code>C = A @ B
print(C)<jupyter_output>[[ 1.+0.j  1.+0.j  1.+0.j  1.+0.j]
 [ 1.+0.j -1.-0.j  0.-1.j -0.+1.j]
 [ 1.+0.j  1.+0.j -1.+0.j -1.+0.j]
 [ 1.+0.j -1.-0.j -0.+1.j  0.-1.j]]<jupyter_text>${\rm\bf C}$ is NOT a sparse matrix.  Task 8: Solving for the Permutation Matrix $${\rm\bf W} = {\rm\bf C} {\rm\bf D}$$${\rm\bf D}$ is the unknown permutation matrix.Because of the nonsingularity of the matrix ${\rm\bf C}$,$${\rm\bf D} = {\rm\bf C}^{-1} {\rm\bf W}$$<jupyter_code>D = np.linalg.inv(C) @ W # using Numpy's function
D = la.inv(C) @ W # using Scipy's function
np.allclose(W, C @ D) # Validate that W = CD
print(D)<jupyter_output>[[ 1.+0.j -0.-0.j -0.+0.j  0.-0.j]
 [ 0.+0.j  0.-0.j  1.+0.j -0.-0.j]
 [ 0.-0.j  1.+0.j -0.-0.j  0.+0.j]
 [ 0.+0.j  0.-0.j  0.-0.j  1.+0.j]]<jupyter_text>Task 9: Vector Permutation with Matrix Operation $${\rm\bf b} = {\rm\bf P}{\rm\bf a}$$$${\rm\bf P} = \left[\begin{matrix}0 & 1 & 0 \\0 & 0 & 1 \\1 & 0 & 0\end{matrix}\right]\quad {\rm\bf a} = \left[\begin{matrix}1 \\2 \\3\end{matrix}\right]\quad {\rm\bf b} = ?$$<jupyter_code>P = np.array([[0, 1, 0], [0, 0, 1], [1, 0, 0]])
a = np.array([1, 2, 3])
b = P @ a
print(b)<jupyter_output>[2 3 1]<jupyter_text>${\rm\bf P}$ is a special matrix, which exhibits only 3 ones, one in each row and one in each column. We call it a permutation matrix. If this kind of matrix is multiplied by any vector, which has the same number of column, it will just rearrange the vector without changing its value.Assume that the **ONES** in a permutation matrix locate in ${\rm\bf x}$th columns, ${\rm\bf y}$th rows. ${\rm\bf x}$ and ${\rm\bf y}$ are two vectors. $${\rm\bf x} = [x_1, x_2, \cdots, x_n]\quad {\rm\bf y} = [y_1, y_2, \cdots, y_n]$$${\rm\bf a}$ is the vector before rerrange, ${\rm\bf b}$ is the vector after rerrange. Then,$${\rm\bf b}[y_n] = {\rm\bf a}[x_n]$$  Task 10: Create Input Signal, Perform Permutation Input signal / column vector$${\rm\bf x}_k = [0,-1,+2,-3,+4,s,-(N-1)]^{\rm T}$$Permutation$${\rm\bf y}={\rm\bf D} {\rm\bf x}_k$$<jupyter_code>x = np.arange(N) * ((-1) ** np.arange(N))
y = D @ x
print(x)
print(y)<jupyter_output>[ 0 -1  2 -3]
[-0.+0.j  2.+0.j -1.-0.j -3.-0.j]<jupyter_text>The lengths of ${\rm\bf x}$ and ${\rm\bf y}$ are the same. Compared to ${\rm\bf x}$, the first N/2 columns of  ${\rm\bf y}$ is the same as the even columns of ${\rm\bf x}$ and the last N/2 columns of  ${\rm\bf y}$ is the same as the odd columns of ${\rm\bf x}$.<jupyter_code>y1 = x[::2] # even columns of x
y2 = x[1::2] # odd columns of x
print(y1)
print(y2)<jupyter_output>[0 2]
[-1 -3]<jupyter_text>We know that ${\rm\bf x}_k$ is a time serie, so we "decimate" a half of elements from this series and put them foreward, then leave the other half backward. In this way, we permutate the input serie.  Task 11: Computing N/2 DFTs When we divide the column matrix  ${\rm\bf y}$ into two parts:  ${\rm\bf y}_1$ and  ${\rm\bf y}_2$, it is much easier to multiply by this sparse matrix without calculating the whole part.$${\rm\bf z} = \left[\begin{matrix}{\rm\bf W}_{N/2} & 0 \\0 & {\rm\bf W}_{N/2}\end{matrix}\right]_{N\times N} {\rm\bf y}=\left[\begin{matrix}{\rm\bf W}_{N/2} & 0 \\0 & {\rm\bf W}_{N/2}\end{matrix}\right] \left[\begin{matrix}{\rm\bf y_1} \\{\rm\bf y_2}\end{matrix}\right] = \left[\begin{matrix}{\rm\bf W}_{N/2} {\rm\bf y_1} \\{\rm\bf W}_{N/2} {\rm\bf y_2}\end{matrix}\right]$$<jupyter_code>z = np.concatenate((Wh @ y1, Wh @ y2))<jupyter_output><empty_output><jupyter_text>The regular N-DFT needs $N^2$ times multiplications and $N(N-1)$ times additions, but the N/2-DFT just needs $N^2/2$ times multiplications and $N(N-2)/2$ times additions.  Task 12: Validating Step $${\rm\bf x}_\mu = {\rm\bf A}{\rm\bf z}$$<jupyter_code>x_mu = A @ z
x_mu_val = fft_np(x) # numpy's FFT
x_mu_val = fft_sp(x) # scipy's FFT
np.allclose(x_mu, x_mu_val)<jupyter_output><empty_output><jupyter_text>Task 13: Combination Matrix $${\rm\bf A} = \left[\begin{matrix}{\rm\bf I}_{N/2} & +{\rm\bf D}_{N/2} \\{\rm\bf I}_{N/2} & -{\rm\bf D}_{N/2}\end{matrix}\right]$$<jupyter_code>print(A)
print(z)
print(x_mu)<jupyter_output>[ 2.+0.j -2.-0.j -4.+0.j  2.+0.j]
[-2.+0.j -2.-2.j  6.+0.j -2.+2.j]<jupyter_text>$${\rm\bf x}_\mu = \left(\begin{matrix}z_1 + z_3 \\z_2 - jz_4 \\z_1 - z_3 \\z_2 + jz_4\end{matrix}\right)$$<jupyter_code>x_man = np.array([z[0]+z[2], z[1]-1j*z[3], z[0]-z[2], z[1]+1j*z[3]])
print(x_man)
np.allclose(x_man, x_mu)<jupyter_output><empty_output><jupyter_text>Task 14: Signal Flow Diagram       Task 15: Recursion If $N=2$,rerunning the cells above (copied the codes in one cell below for clear representation).<jupyter_code>N = 2 ** 1 # must be even, should be power of two
Nh = N // 2 # interger division
W = dft_matrix(N) # full DFT matrix for analysis of signal with length N
Wh = dft_matrix(Nh) # half size DFT matrix that actually needs to be computed

Ih = np.identity(Nh)
b = np.arange(Nh)
Dh = np.diag(np.exp(-1j * (2 * np.pi / N) * b))

A = np.vstack((np.hstack((Ih, Dh)), np.hstack((Ih, -Dh))))

B = np.zeros((N, N), dtype=complex)
B[:Nh, :Nh] = Wh
B[Nh:N, Nh:N] = Wh

C = A @ B

D = np.linalg.inv(C) @ W # using Numpy's function

x = np.arange(N) * ((-1) ** np.arange(N))
y = D @ x

y1 = x[::2] # even columns of x
y2 = x[1::2] # odd columns of x

z = np.concatenate((Wh @ y1, Wh @ y2))

x_mu = A @ z
print(D)
print(B)<jupyter_output>[[1.+0.j 0.-0.j]
 [0.+0.j 1.+0.j]]
[[1.+0.j 0.+0.j]
 [0.+0.j 1.+0.j]]<jupyter_text>${\rm\bf D}$ and ${\rm\bf B}$ are both identity matrices, so the input vector ${\rm\bf x}_k$ is unaltered.<jupyter_code>print(A)<jupyter_output>[[ 1.+0.j  1.+0.j]
 [ 1.+0.j -1.-0.j]]
{'fraction_non_alphanumeric': 0.2353542234332425, 'fraction_numerical': 0.04757039055404178, 'mean_word_length': 3.4557410217501263, 'pattern_counts': {'":': 0, '<': 69, '<?xml version=': 0, '>': 69, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '20', 'dir': 'jupyter-structured-clean-dedup', 'id': '130457', 'license': 'MIT', 'n_tokens_mistral': 4275, 'n_tokens_neox': 4000, 'n_words': 983, 'path': 'Task_FFT.ipynb', 'repo_name': 'wyzdwdz/digital-signal-processing-homework', 'size': '8808'}