[Python] what's wrong with this Python assignment on signal processing - mostly Fourier series...

I'll try and detail as much as possible, please ask me if any info is missing in your opinions.

in this assigment i created the basic rect signal a[n] such that over the domain [-1000,1000] it's 1 only at |n|<100, meaning it's an array (complex one with zero in the imag part) that looks like this [0,0,...0,0,1,1,...,1,1,0,...0,0,0] where there are exactly 199 ones, 99 on positive and negative and one at 0, looks like that:



I've created also the following FT function, and a threshold function to clean Floating point errors from the results:

import numpy as np
import cmath
import matplotlib.pyplot as plt

D=1000
j = complex(0, 1)
pi = np.pi
N = 2 * D + 1

a=np.zeros(2*D+1)
for i in range(-99,100):
a[i+D] = 1
threshold = 1e-10
def clean_complex_array(arr, tol=threshold):
real = np.real(arr)
imag = np.imag(arr)

# Snap near-zero components
real[np.abs(real) < tol] = 0
imag[np.abs(imag) < tol] = 0
# Snap components whose fractional part is close to 0 or 1
real_frac = real - np.round(real)
imag_frac = imag - np.round(imag)

real[np.abs(real_frac) < tol] = np.round(real[np.abs(real_frac) < tol])
imag[np.abs(imag_frac) < tol] = np.round(imag[np.abs(imag_frac) < tol])

return real + 1j * imag

def fourier_series_transform(data, pos_range, inverse=False):
full_range = 2 * pos_range + 1
# Allocate result array
result = np.zeros(full_range, dtype=complex)

if inverse:
# Inverse transform: reconstruct time-domain signal from bk
for n in range(-pos_range, pos_range+ 1):
for k in range(-pos_range, pos_range+ 1):
result[n + pos_range] += data[k + pos_range] * cmath.exp(j * 2 * pi * k * n / full_range)
else:
# Forward transform: compute bk from b[n]
for k in range(-pos_range, pos_range+ 1):
for n in range(-pos_range, pos_range+ 1):
result[k + pos_range] += (1 / full_range) * data[n + pos_range] * cmath.exp(-j * 2 * pi * k * n / full_range)

return result

ak = fourier_series_transform(a, D)
ak = clean_complex_array(ak)


a_k looks like that: (a real sinc signal, which is to be expected)



i've checked that the threshold value is good, FPE starts at around e-14 and there's no significant contributions to the signal below e-8.

now for the part i had a problem with: we're asked to create the freq signal f_k such that f_k will be a_k padded with 4 zeros after each value and multiplied by 0.2, meaning it will look like this 0.2*[a_0,0,0,0,0,a_1,0,0,0,0,a_2,0,0,0,0,a_3,...], we want to show that doing so equals a streching of the signal in the time domain.

now when i did the math it checks out, you get 5 copies of the original signal over a range of [-5002,5002] (to create 10005 samples which is exactly 5*2001 which was the original number of samples of the signals), the following is the code for this section, to set f_k and f[n]:

stretch_factor = 5
f_k = np.zeros(stretch_factor * N, dtype=complex)
f_k[::stretch_factor] = 0.2 * ak # scale to keep energy in check
# New domain size after stretching
D_new = (len(f_k) - 1) // 2
# Inverse transform to get f[n]
f_n = fourier_series_transform(f_k, D_new, inverse=True)
f_n = clean_complex_array(f_n)


plt.figure()
plt.plot(np.arange(-D_new, D_new + 1), np.real(f_n), label='Real part')
plt.plot(np.arange(-D_new, D_new + 1), np.imag(f_n), label='Imaginary part', color='red')
plt.grid(True)
plt.title("Compressed signal $f[n]$ after frequency stretching")
plt.xlabel("n")
plt.ylabel("Amplitude")
plt.legend()


and this is what i get:



which is wrong, i should be getting a completly real signal, and as i said it should be 5 identical copies at distance of 2000 from each other, something like that:



I dont know why it does that, and i even tried to use AI to explain why it happens and how to fix it and it couldn't help with either, i would appriciate help here.

Continue reading...