I am new here and also am not very knowledgeable about DSP so this might be dumb and easy. I am aware of the fact that in order to reconstruct a signal, I need to sample it at a frequency that is more than twice its highest frquency.
In the example below, I use N=50
samples and expect the reconstruction to fail if I use freq = 25
or so (perhaps due to numerical round-off). But I get that the reconstruction breaks around freq = 12.5
. I feel I am missing something fundamental.
Here are plots and the code to generate them:
from numpy.fft import rfft as rfft
from numpy.fft import irfft as irfft
import numpy as np
import matplotlib.pyplot as plt
import math
def interpolant( f_hat, pts ):
'''
This is the trigonometric polynomial interpolating
the signal that is encoded in the values of f_hat.
it is evaluated at the points pts
The commented part is the true calculation to be done,
the uncommented just does it more accurately.
'''
N = len( f_hat )
f = f_hat[0]
for k in range( 1, N/2 ):
f = f + 2 * f_hat[k].real * np.cos( k * pts )
f = f - 2 * f_hat[k].imag * np.sin( k * pts )
#f = f + f_hat[k] * np.exp( 1j * k * pts )
#f = f + f_hat[k].conjugate() * np.exp(-1j * k * pts )
f = f + f_hat[N/2] * np.cos( N/2 * pts )
return f / N / 2
def f( pts, freq ):
'''
this is "the conitinuous" signal, evaluated (sampled)
at pts
'''
return np.sin( freq * pts )
# Number of sampled points
N = 50
# Where we sample the signal
pts = np.linspace(0, 2 * np.pi, num = N, endpoint = False )
# A much finer grid, used solely to display results
oversampled = np.linspace( 0, 2 * np.pi, 20 * N, endpoint = False )
# The grid we use for the plotting
grid = oversampled
# The frequncy of the sine wave below
freq = 7
# Do the FFT on the sampled signal
f_hat = rfft( f( pts, freq ) )
plt.plot( grid, interpolant( f_hat, grid ) , color = "g" )
plt.plot( grid, f( grid, freq ) , color = "r" )
title1 = str(N) + " samples. Signal frequency is " + str(freq)+" \n"
title2 = "Red is true, green is interpolant. Reconstruction succeeds"
plt.title( title1 + title2 )
plt.show()
# The frequncy of the sine wave below
freq = 12.9
# Do the FFT on the sampled signal
f_hat = rfft( f( pts, freq ) )
plt.plot( grid, interpolant( f_hat, grid ) , color = "g" )
plt.plot( grid, f( grid, freq ) , color = "r" )
title1 = str(N) + " samples. Signal frequency is " + str(freq)+" \n"
title2 = "Red is true, green interpolant. Reconstruction fails"
plt.title( title1 + title2 )
plt.show()