using PyPlot

freq = 440.0
N = 256
duration = 8 / freq
t = linspace(0, duration, N)
x = cos(2 * pi * freq * t)

plot(t, x)

X = rfft(x)
sampRate = N / duration
f = linspace(0, sampRate / 2, int(N / 2) + 1)
magX = abs(X)
plot(f, magX)

chord = cos(2 * pi * freq * t) + 0.5 * cos(4 * pi * freq * t)

plot(t, chord)

Chord = abs(rfft(chord))
plot(f, Chord)

noise = randn(length(x)) / 20
noisy_x = x + noise
plot(t, noisy_x)

mavg_n = 16
mavg_h = ones(Float64, mavg_n) / mavg_n
# conv is the convolution function
filtered_x = conv(noisy_x, mavg_h)
t = [0:length(filtered_x) - 1] / sampRate
plot(t, filtered_x)

sinc_w = 75
t = [-sinc_w:sinc_w] / sampRate
freqCrit = 660.0
sinc_h = sin(2 * pi * freqCrit * t) ./ (2 * pi * freqCrit * t)
# need to fix the singularity at the origin
sinc_h[sinc_w + 1] = 1
# normalize so that the sum is 1.0
sinc_h ./= sum(sinc_h)
plot(t, sinc_h)

sinc_H = abs(rfft(sinc_h))
f = linspace(0, sampRate / 2, sinc_w + 1)
plot(f, sinc_H)

M = 2 * sinc_w
i = [0:M]
hamming_w = 0.54 - 0.46 * cos(2 * pi * i / M)
plot(t, hamming_w)

win_sinc_h = hamming_w .* sinc_h
win_sinc_h ./= sum(win_sinc_h)
plot(t, win_sinc_h)

Win_Sinc_H = abs(rfft(win_sinc_h))
plot(f, Win_Sinc_H)

filtered_x = conv(chord, win_sinc_h)
t = [0:length(filtered_x)-1] / sampRate
plot(t, filtered_x)

filtered_X = abs(rfft(filtered_x))
f = linspace(0, sampRate / 2, length(filtered_X))
plot(f, filtered_X)

delta = zeros(Float32, 2 * sinc_w + 1)
delta[sinc_w + 1] = 1.0
t = [-sinc_w:sinc_w] / sampRate
plot(t, delta)

Delta = abs(rfft(delta))
# set the edges to different numbers so you can clearly see the scale
f = linspace(0, sampRate / 2, sinc_w + 1)
ylim([0.0, 1.1])
plot(f, Delta)

highp_h = delta - win_sinc_h
plot(t, highp_h)

Highp_H = abs(rfft(highp_h))
plot(f, Highp_H)

filtered_x = conv(chord, highp_h)
t = [0:length(filtered_x) - 1] / sampRate
plot(t, filtered_x)

Filtered_X = abs(rfft(filtered_x))
f = linspace(0, sampRate / 2, length(Filtered_X))
plot(f, Filtered_X)