Fourier Transforms

# Create array with zeros
y = np.zeros(len(x))

# Change zeros to 1 based on a given period
# in radians (e.g. 2pi, 3pi)
for i in range(len(x)):
    if (x[i]/(period)) % 1 < 0.50:
        y[i] = 1.0
return y

Now let's try to express the square wave as a Fourier series. First, we calculate the coefficients using the expressions above. Carrying out the calculations, we find

$$a_0 = \frac{1}{\pi} \int_{0}^{\pi} dx = 1,\\
a_n = \frac{1}{\pi} \int_{0}^{\pi} \cos(nx)dx = 0, \,\,\,\,n\geq 1,$$

$$
b_n = \frac{1}{\pi} \int_{0}^{\pi} \sin(nx)\,dx = \left\{
    \begin{array}{ll}
        \frac{2}{n\pi} & \text{if } n \text{ is odd}, \\
        0 & \text{if } n \text{ is even}. \\
    \end{array}
\right.
$$


Now let's plot the Fourier series for $ n=10, n=50,$ and $n=100$.

```{code-cell} python
:tags: [hide-input]
# Define the cos and sin terms of the Fourier series
def cosTerm(n):
    # Always zero except for n=0
    if n==0: return 1.0
    return 0.

def sinTerm(n):

    if n%2: # n modulo 2 = 1 (True) then Odd
        ret = 2. / (n* np.pi)
    else:
        ret = 0.
    return ret

def fourier(n,x):
    #a_0 term, remember 1/2
    sum = cosTerm(0)/2.0 * np.ones(len(x))
    
    #all other terms
    for i in range(1, n+1):
        sum += sinTerm(i)*np.sin(i*x) + cosTerm(i)*np.cos(i*x)
    return sum

fig, axes = plt.subplots(1, 3, figsize=(14,4), sharey=True)

# Loop over each subplot
for idx, i in enumerate([10, 50, 100]):
    ax = axes[idx]
    ax.plot(x, fourier(i, x), label="n = %g" % (i))
    ax.plot(x, square(x, 2*np.pi), color="black", label="square wave")
    ax.legend(loc="upper right", fontsize="small")
    ax.set_xlabel("x")
    ax.set_ylim(-0.2, 1.2)
    
plt.tight_layout()
plt.show()

- Notice that it is not a *perfect* reconstruction. We only put in two specific frequencies, so we would expect to see two sharp lines. But if the waves are not exactly periodic when we cut off the time signal, the peaks acquire a bit of width.

- The time-domain signal can be as complicated as we want. Here is an example with four frequencies. Note that the *phase*, whether indicated by a sine or a cosine function, does not matter: both are periodic.

#### Noisy periodic data

```{code-cell} python
:tags: [hide-input]
import numpy as np
import matplotlib.pyplot as plt
from scipy.fftpack import fft

# Generate sample data
t = np.linspace(0.0, 1.0, 500)
y = np.sin(30.0 * 2.0 * np.pi * t) + 0.5 * np.cos(50.0 * 2.0 * np.pi * t)
y = y + 1.5 * np.sin(165.0 * 2.0 * np.pi * t) + np.cos(104 * 2.0 * np.pi * t)

# Compute the Fourier Transform
N = len(t)
yf = fft(y)
xf = np.fft.fftfreq(N, t[1] - t[0])[:N//2]

# Plot both time domain and frequency domain side by side
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(12, 6))

# Time domain plot
ax1.plot(t, y)
ax1.set_xlabel("Time")
ax1.set_ylabel("Amplitude")
ax1.set_title("Time Domain Signal")

# Frequency domain plot
ax2.plot(xf, 2.0 / N * np.abs(yf[:N//2]))
ax2.set_xlabel("Frequency")
ax2.set_ylabel("Amplitude")
ax2.set_title("Frequency Domain Signal")

# Show the plots
plt.tight_layout()
plt.show()

- Notice that while the peak heights changed a bit, the relative intensities and peak areas remain unchanged. So if you are doing spectroscopy, you need a few wave repeats, but it need not be a "continuous wave" spectrum. This is useful if you are worried about applying a high-intensity beam (e.g. an x-ray).

### Fourier transform and the uncertainty relation

**Define the Square Function to represent $\psi(x)$**: Start with a square function (or rectangular pulse) that has a width $w$ and height 1, centered around zero. This function can be represented as:

$$
   \psi(x) = \begin{cases} 
   1, & |x| \leq \frac{w}{2} \\
   0, & |x| > \frac{w}{2}
   \end{cases}
$$
   
**Fourier Transform gives us $\psi(p)$**: The Fourier transform $\hat{f}(k)$ of a square function is a sinc function:

   $$
   \hat{\psi}(k) = w \, \text{sinc}(kw/2)
   $$

- where $\text{sinc}(x) = \frac{\sin(x)}{x}$. This shows how the Fourier transform spreads out as the width of the square function changes.

:::{important} **Uncertainty Relation**

   $$
   \Delta x \Delta k \geq \frac{1}{2}
   $$
   
:::

- The uncertainty principle states that the product of the **standard deviations** of a **function and its Fourier transform** is **bounded by a minimum value**!
- Here, $\Delta x$ is the width of the function (related to the spread in position), and $ \Delta k$ is the spread of the Fourier transform (momentum space). As the width $w$ of the square function decreases, the spread of the sinc function in the Fourier domain increases, demonstrating the uncertainty relation.
- In quantum mechanics, using the de Broglie relation, we can recover the uncertainty relation since $k=2\pi/\lambda = 2\pi p/h = p/\hbar$.

```{code-cell} python
:tags: [hide-input]
import matplotlib.pyplot as plt
from matplotlib.animation import FuncAnimation
from IPython.display import HTML
import numpy as np

# Define the square function
def square_function(x, w):
    return np.where(np.abs(x) <= w / 2, 1, 0)

# Fourier transform of the square function (sinc function)
def fourier_transform_square_function(k, w):
    return w * np.sinc(k * w / (2 * np.pi))

# Define x and k ranges
x = np.linspace(-5, 5, 1000)
k = np.linspace(-20, 20, 1000)

# Set up the figure and axes for animation
fig, axes = plt.subplots(2, 1, figsize=(8, 8))

# Initialize the plots
line1, = axes[0].plot([], [], lw=2)
axes[0].set_ylim([-0.1, 1.1])
axes[0].set_title("Square Function")

line2, = axes[1].plot([], [], lw=2)
axes[1].set_ylim([-0.5, 2.5])
axes[1].set_title("Fourier Transform")

# Set the x-axis limits
axes[0].set_xlim(x.min(), x.max())
axes[1].set_xlim(k.min(), k.max())

# Initialization function
def init():
    line1.set_data([], [])
    line2.set_data([], [])
    return line1, line2

# Animation update function
def update(w):
    # Compute the square function and its Fourier transform
    sq_func = square_function(x, w)
    ft_sq_func = fourier_transform_square_function(k, w)
    
    # Update the data for both plots
    line1.set_data(x, sq_func)
    line2.set_data(k, ft_sq_func)
    
    return line1, line2

# Create the animation
ani = FuncAnimation(fig, update, frames=np.linspace(0.1, 5.0, 100), init_func=init, blit=True)

plt.close(fig)  # Prevents static display of the last frame
HTML(ani.to_jshtml())

### Delta Function

The delta function, $\delta(t)$, is defined as

$$
\delta(t) = \left\{
    \begin{array}{ll}
        0 & \text{if } t \neq 0, \\
        \infty & \text{if } t = 0.
    \end{array}
\right.
$$


Carrying out the transform, we see that the Fourier transform of the delta function is actually a constant, such that

$$ \mathcal{F}[{\delta}(t)] = 1 $$

Let's plot the function and its transform.

```{code-cell} python
:tags: [hide-input]
from scipy import signal
imp = signal.unit_impulse(100, 'mid') # creates the delta function
t = np.linspace(-5, 5, 100)

plt.subplot(1,2,1)
plt.plot(t, imp)
plt.title('Delta Function')

FT_omega = np.fft.fftfreq(100, t[1] - t[0])
FT = np.fft.fft(imp)

plt.subplot(1,2,2)
plt.plot(FT_omega, abs(FT))
plt.title('Fourier Transform')
plt.ylim([0,1.2])
plt.tight_layout()
plt.show()