In [1]:
import matplotlib.pyplot as plt
import numpy as np

Review wave basics in python¶

In [2]:
%matplotlib inline
In [3]:
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¶

In [4]:
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¶

In [5]:
# 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()
In [6]:
# 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.

In [7]:
# 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()
In [8]:
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
In [9]:
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¶

In [10]:
# number of samples
N = len(x1)
N
Out[10]:
100
In [11]:
# sample frequencies
k = np.arange(N)
# time component
t = k / N
# take the sample frequencies as column
k = k.reshape((N, 1))
In [12]:
# 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()
In [13]:
[(i, item) for i, item in enumerate(X1) if abs(item) > 1]
Out[13]:
[(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

In [14]:
# 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.

In [15]:
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