Nyquist Frquency isn't what it is supposed to be 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()
The closer the sampling frequency gets to half the (highest spectral) frequency of the signal, the longer one has to sample to get a sufficient bounded fit (for some bound) of a reconstruction using those finite number of samples.  This required sampling time extends to infinity as one approaches Fs/2.  

So, you simply didn't sample for long enough in your experiment.  Try 5 thousand and 5 million samples for sampling a 12.9 Hz sinusoid at 26 sps, and you will see a closer fitting looking plot over most of an X axis of 5 thousand or better million reconstructed sample points.  Then try a sinusoid even closer to Fs/2 and you will find you need to sample even longer to get the same visual quality of fit.