Algorithm - return aliasing frequency

Question

In Python, I'm trying to write an algorithm alias_freq(f_signal,f_sample,n), which behaves as follows:

def alias_freq(f_signal,f_sample,n):
    f_Nyquist=f_sample/2.0
    if f_signal<=f_Nyquist:
        return n'th frequency higher than f_signal that will alias to f_signal
    else:
        return frequency (lower than f_Nyquist) that f_signal will alias to

The following is code that I have been using to test the above function (f_signal, f_sample, and n below are chosen arbitrarily just to fill out the code)

import numpy as np
import matplotlib.pyplot as plt

t=np.linspace(0,2*np.pi,500)

f_signal=10.0
y1=np.sin(f_signal*t)
plt.plot(t,y1)

f_sample=13.0
t_sample=np.linspace(0,int(f_sample)*(2*np.pi/f_sample),f_sample)
y_sample=np.sin(f_signal*t_sample)
plt.scatter(t_sample,y_sample)

n=2
f_alias=alias_freq(f_signal,f_sample,n)
y_alias=np.sin(f_alias*t)
plt.plot(t,y_alias)

plt.xlim(xmin=-.1,xmax=2*np.pi+.1)
plt.show()

My thinking is that if the function works properly, the plots of both y1 and y_alias will hit every scattered point from y_sample. So far I have been completely unsuccessful in getting either the if statement or the else statement in the function to do what I think it should, which makes me believe that either I don't understand aliasing nearly as well as I want to, or my test code is no good.

My questions are: Prelimarily, is the test code I'm using sound for what I'm trying to do? And primarily, what is the alias_freq function that I am looking for?

Also please note: If some Python package has a function just like this already built in, I'd love to hear about it - however, part of the reason I'm doing this is to give myself a device to understand phenomena like aliasing better, so I'd still like to see what my function should look like.

Answer 1

Here is a Python aliased frequency calculator based on numpy

def get_aliased_freq(f, fs):
    """
    return aliased frequency of f sampled at fs
    """
    import numpy as np

    fn = fs / 2
    if np.int(f / fn) % 2 == 0:
        return f % fn
    else:
        return fn - (f % fn)

Answer 2

As far as I understood the question correctly, the frequency of the aliased signal is abs(sampling_rate * n - f_signal), where n is the closest integer multiple to f_signal.

Thus:

n = round(f_signal / float(f_sample))
f_alias = abs(f_sample * n - f_signal)

This should work for frequencies under and over Nyquist.

Answer 3

I figured out the answer to my and just realized that I forgot to post it here, sorry. Turns out it was something silly - Antii's answer is basically right, but the way I wrote the code I need a f_sample-1 in the alias_freq function, where I just had an f_sample. There's still a phase shift thing that happens sometimes, but just plugging in either 0 or pi for the phase shift has worked for me every time, I think it's just due to even or odd folding. The working function and test code is below.

import numpy as np
import matplotlib.pyplot as plt

#Given a sample frequency and a signal frequency, return frequency that signal frequency will be aliased to.
def alias_freq(f_signal,f_sample,n):
    f_alias = np.abs((f_sample-1)*n - f_signal)
    return f_alias

t=np.linspace(0,2*np.pi,500)

f_signal=13
y1=np.sin(f_signal*t)
plt.plot(t,y1)

f_sample=7
t_sample=np.linspace(0,int(f_sample)*(2*np.pi/f_sample),f_sample)
y_sample=np.sin((f_signal)*t_sample)
plt.scatter(t_sample,y_sample)

f_alias=alias_freq(f_signal,f_sample,3)
y_alias=np.sin(f_alias*t+np.pi)#Sometimes with phase shift, usually np.pi for integer f_signal and f_sample, sometimes without.
plt.plot(t,y_alias)

plt.xlim(xmin=-.1,xmax=2*np.pi+.1)
plt.show()