import matplotlib.pyplot as plt
import numpy as np

Review wave basics in python

%matplotlib inline

x = np.linspace(0, 20, 201)
y = np.sin(x)

plt.figure(figsize = (8, 6))
plt.plot(x, y, 'b')
plt.ylabel('Amplitude')
plt.xlabel('Location (x)')
plt.show()

Look at the wave at different time

fig = plt.figure(figsize = (8,8))

times = np.arange(5)

n = len(times)

for t in times:
    plt.subplot(n, 1, t+1)
    y = np.sin(x + t)
    plt.plot(x, y, 'b')
    plt.plot(x[25], y [25], 'ro')
    plt.ylim(-1.1, 1.1)
    plt.ylabel('y')
    plt.title(f't = {t}')

plt.xlabel('location (x)')
plt.tight_layout()
plt.show()

The function for a sine wave

y(t)=Asin(ωt+ϕ)

ω=2π/T=2πf

Discrete sample's plot

# sampling rate in Hz
sr = 100.0
# sampling interval
ts = 1.0/sr
t = np.arange(0,1,ts)

# frequency of the signal
freq = 5   
y = np.sin(2*np.pi*freq*t)

plt.figure(figsize = (8, 8))
plt.subplot(211)
plt.plot(t, y, 'b')
plt.ylabel('Amplitude')

freq = 10   
y = np.sin(2*np.pi*freq*t)

plt.subplot(212)
plt.plot(t, y, 'b')
plt.ylabel('Amplitude')

plt.xlabel('Time (s)')
plt.show()

# frequency of the signal
freq = 5   
y = 5*np.sin(2*np.pi*freq*t)

plt.figure(figsize = (8, 8))
plt.subplot(211)
plt.plot(t, y, 'b')
plt.ylabel('Amplitude')

y = 10*np.sin(2*np.pi*freq*t + 10)

plt.subplot(212)
plt.plot(t, y, 'b')
plt.ylabel('Amplitude')

plt.xlabel('Time (s)')
plt.show()

The DFT

Sine waves with frequencies: 1 Hz, 4 Hz, and 7 Hz amplitudes: 3, 1 and 0.5 and phase all zeros. Sampling rate 100 Hz.

# sampling rate
sr = 100
# sampling interval
ts = 1.0/sr
t = np.arange(0,1,ts)

freq = 1.
x1 = 3*np.sin(2*np.pi*freq*t)

freq = 4
x2 = np.sin(2*np.pi*freq*t)

freq = 7   
x3 = 0.5* np.sin(2*np.pi*freq*t)

x = x1 + x2 + x3

plt.figure(figsize = (8, 6))
plt.subplot(411)
plt.plot(t, x1, 'g')
plt.subplot(412)
plt.plot(t, x2, 'g')
plt.subplot(413)
plt.plot(t, x3, 'g')
plt.subplot(414)
plt.plot(t, x, 'r')
plt.ylabel('Amplitude')

plt.show()

def DFT(x):
    """
    Function to calculate the 
    discrete Fourier Transform 
    of a 1D real-valued signal x
    """

    N = len(x)
    n = np.arange(N)
    k = n.reshape((N, 1))
    e = np.exp(-2j * np.pi * k * n / N)
    
    X = np.dot(e, x)
    
    return X

X = DFT(x)

# calculate the frequency
N = len(X)
n = np.arange(N)
T = N/sr
freq = n/T 

plt.figure(figsize = (8, 6))
plt.stem(freq, abs(X), 'b', \
         markerfmt=" ", basefmt="-b")
plt.xlabel('Freq (Hz)')
plt.ylabel('DFT Amplitude |X(freq)|')
plt.show()

Understanding the DFT

Analyse x1

# number of samples
N = len(x1)
N

100

# sample frequencies
k = np.arange(N)
# time component
t = k / N
# take the sample frequencies as column
k = k.reshape((N, 1))

# exponent gives an NxN matrix for each "t" and "k"
e = np.exp(-2j * np.pi * k * t)
# inner product of the vectors in e and x1;
# in other words, takes the sum of all samples "n" for a given frequency sample "k"
# this accounts for the positive and the negative values
X1 = np.dot(e, x1)
# total periods
T = N / sr
# frequencies
freq = np.arange(N) / T
plt.figure(figsize = (8, 6))
plt.stem(freq, abs(X1), 'b', \
         markerfmt=" ", basefmt="-b")
plt.xlabel('Freq (Hz)')
plt.ylabel('DFT Amplitude |X(freq)|')
plt.show()

[(i, item) for i, item in enumerate(X1) if abs(item) > 1]

[(1, (7.673197244665952e-16-150.00000000000003j)),
 (99, (9.013836182311416e-13+150.00000000000006j))]

Why and how does the values at index 1, 99 differs from the others? Doing it manually!

We want to calculate the Xk for each k independently of each other and look at the sine and cosine component of their values.

Let's take x2 this time and try to explain it piece by piece.

The N is same for all of the waves in this notebook so we don't calculate it again

# Multiplicative frequencies
k = np.arange(N)
# Time for corresponding sample as a column
t = k / N

For each sample of k (multiplicative frequency) at index (of each amplitude sample) calculate the sum of all dot(x2[index], position_vector_of_k (at index)).

Since all the x values are amplitudes of a sine wave in vertical axis we take only take the imaginary part of the position_vector_of_k from (corresponding to the vertical axis of) the complex plane, which corresponds to the "dot product". And will not be True for phase shifted waves, and in which case we will have to account for the phase shift and project the position_vector_of_k to the position_vector_of_the_original_sample. np.dot does it automatically given an NxN matrix of position_vector_of_k at time t.

nhfloor = N // 2
def apply_mfs(orig_sample, mfs, cosine=False):
    plt.figure(figsize=(8, N))
    # for each multiplicative frequency, mf; do (
    #     amplitude of mf * value of n-the sample(which already represents the amplitude)
    # )
    for index, mf in enumerate(mfs):
        amp_vectors = np.array([0.0 for _ in range(N)])
        for i, ti in enumerate(t):
            # position vector of the mf at time ti on the complex plane.
            pvki = np.exp(-2j * np.pi * mf * ti)
            # the projection of the position vector (in the direction of the position of the original sample)
            # times the original frequency amplitude at time ti
            if cosine:
                amp_vectors[i] = orig_sample[i] * pvki.real
            else:
                amp_vectors[i] = orig_sample[i] * pvki.imag
        if index + 1 <= nhfloor:
            marker = "bo"
            if abs(amp_vectors.sum()) > 1:
                print(mf, amp_vectors.sum())
                marker = "go"
            plt.subplot(nhfloor, 1, index+1)
            plt.stem(t, amp_vectors, marker)
    plt.show()
apply_mfs(x2, k)

4 -50.000000000000014

The above figure contains superimposition of all N // 2 frequencies to x2.

What is happening in the above plot?

We are basically superimposing N frequencies to the original frequency.

From the superimposition, if we take the sum of the area under the curve in the direction of the base line, for all frequency spectrum that does not match the frequency of the original frequency, become close to ZERO, given the period of the original frequency is less then the sampling interval. Also, We don't really take the area under the curve, since we are calculating them in a discrete manner, we only take the sum of magnitudes from the superimposition at each sample.

There should be some general pattern in the frequencies to superimpose, for the moment we consider that they must be N samples of n frequency. But the exception to this is the exploitation of symmetries in sinusoids, we will explore that in FFT section.

We will now try the above method with shifted phase of x2.

plt.subplot(211)
plt.plot(t, x2)
# we shift the x2 sample by 3 samples to the left
x2shifted = np.roll(x2, -3)
plt.subplot(212)
plt.plot(t, x2shifted)
plt.show()

apply_mfs(x2shifted, k)

4 -36.44843137107058

The above figure contains superimposition of all N // 2 frequencies to the shited x2

We will no shift x2 by 90 degrees, and apply the superimposition.

plt.subplot(211)
plt.plot(t, x2)
# shifting to left by 6 samples gives something similar to a cosine wave
x2shifted = np.roll(x2, -6)
plt.subplot(212)
plt.plot(t, x2shifted)
plt.show()

apply_mfs(x2shifted, k, cosine=True)

4 49.90133642141358

Since the current shifted x corresponds to a cosine wave, we pass in True for the cosine value, which takes the projection of position_vector_of_k in the horizontal axis. Justification for this is as such that since sample ZERO starts from 90 degrees, we can simply shift the perspective of the complex plane of k by 90 degrees.

Apply the superimposition to the compound wave

apply_mfs(x, k)

1 -150.0
4 -50.000000000000014
7 -25.000000000000068

FFT

Thanks to: https://pythonnumericalmethods.berkeley.edu