In [ ]:
import numpy as np
from scipy import fftpack as fft
import matplotlib.pyplot as plt
plt.style.use('ggplot')
%matplotlib inline

DTMF: Linear combination of two sinusoids


In [ ]:
Fs = 32768
duration = 0.25
t = np.linspace(0, duration, duration * Fs)
f1, f2 = 697, 1336
y1 = np.sin(2 * np.pi * f1 * t);
y2 = np.sin(2 * np.pi * f2 * t);
y = (y1 + y2) / 2
plt.plot(t, y)

In [ ]:
from IPython.display import Audio
Audio(y, rate=44100)

Repeating the signal ten times


In [ ]:
t = np.linspace(0, duration * 10, duration * 10 * Fs)
f1, f2 = 697, 1336
y1 = np.sin(2 * np.pi * f1 * t);
y2 = np.sin(2 * np.pi * f2 * t);
y = (y1 + y2) / 2
Audio(y, rate=44100)

In [ ]:
# Recreate the original signal for simplicity

t = np.linspace(0, duration , duration * Fs)
f1, f2 = 697, 1336
y1 = np.sin(2 * np.pi * f1 * t);
y2 = np.sin(2 * np.pi * f2 * t);
y = (y1 + y2) / 2

n = y.shape[0]
p = np.abs(fft.fft(y));
f = fft.fftfreq(n, d=1/Fs)
plt.plot(f,p);

Q: Why did this happen?

Exercise: Re-plot the spectrum for only positive frequencies, and limit the X-axis to only 2 kHz


In [ ]:
# enter code here

Check that peaks are at the correct frequencies


In [ ]:
max_power_index = np.argsort(p)[::-1]
max_frequency_index = f[max_power_index]
print(max_frequency_index[:5])

Processing a noisy signal


In [ ]:
time_step = 0.02
period = 5.
time_vec = np.arange(0, 20, time_step)
sig = np.sin(2 * np.pi / period * time_vec) + \
      0.5 * np.random.randn(time_vec.size)
plt.plot(time_vec, sig)

In [ ]:
sample_freq = fft.fftfreq(sig.size, d=time_step)
sig_fft = fft.fft(sig)
pidxs = np.where(sample_freq > 0)
freqs = sample_freq[pidxs]
power = np.abs(sig_fft)[pidxs]
plt.plot(freqs, power)
plt.xlim(0, 5)

In [ ]:
# denoising
freq = freqs[power.argmax()]
sig_fft[np.abs(sample_freq) > freq] = 0
# Reconstruction
recons = fft.ifft(sig_fft)
plt.plot(time_vec, sig, time_vec, recons)

exercise: 2D FFT - Denoise the follwing image

Hint: Look for 2D FFT functions in scipy.fftpack module

1. Visualize the frequency spectrum of the image

2. Threshold on a suitable value


In [ ]:
x = plt.imread("moonlanding.png")
plt.imshow(x, cmap=plt.cm.gray)
plt.yticks([])
plt.xticks([])
plt.grid()

In [ ]:
# enter code here