How to change the samplimg rate with interpolation and decimation?

Question

I have a wav file at a 44.1KHz. I'm trying to change the sampling frequency to 1.26MHz. For that, I need to use interpolation and decimation, and so I did, but I'm getting odd results back. It seems like there are a few spikes when I do the FFT although there shouldn't be.

Here is my function:

import scipy
import numpy as np
import matplotlib.pyplot as plt
import scipy.io.wavfile
from scipy.fftpack import fft
import scipy.signal as sig

def chunks(lst, n):
    """Yield successive n-sized chunks from lst."""
    for i in range(0, len(lst), n):
        yield lst[i:i + n]

def SamplerateConversion(samples, new_fs, old_fs):
    L = int(new_fs / 100)
    M = int(old_fs / 100)

    hcf = compute_hcf(L,M)
    L = int(L / hcf)
    M = int(M / hcf)

    if L > 1 and M == 1:
        inter = interpolation(samples,L,old_fs)
        w = low_cut_filter(inter,new_fs,22050)
        return w
    elif L == 1 and M > 1:
        dec = decimation(samples,M,old_fs)
        return dec
    elif L > 1 and M > 1:
        inter = interpolation(samples,L,old_fs)
        filtered_inter = low_cut_filter(inter,old_fs * L,22050)
        dec = decimation(filtered_inter,M,old_fs)
        return dec

def interpolation(samples,num, old_fs):
    inter_samples = []

    chunks_samp = list(chunks(samples,old_fs))
    for chunk in chunks_samp:
        for j in chunk:
            inter_samples.append(j)
            for h in range(num-1):
                inter_samples.append(0)
    return inter_samples

def decimation(samples,num, old_fs):
    inter_samples = []
    i = num - 1
    chunks_samp = list(chunks(samples,old_fs))
    for chunk in chunks_samp:
        for j in chunk:
            if i == num-1:
                inter_samples.append(j)
                i = 0
            else:
                i += 1
    return inter_samples

def compute_hcf(x, y):
# choose the smaller number
    if x > y:
        smaller = y
    else:
        smaller = x
    for i in range(1, smaller+1):
        if((x % i == 0) and (y % i == 0)):
            hcf = i 
    return hcf

def low_cut_filter(x, fs, cutoff=70):
    """FUNCTION TO APPLY LOW CUT FILTER

    Args:
        x (ndarray): Waveform sequence
        fs (int): Sampling frequency
        cutoff (float): Cutoff frequency of low cut filter

    Return:
        (ndarray): Low cut filtered waveform sequence
    """

    nyquist = fs // 2
    norm_cutoff = cutoff / nyquist

    # low cut filter
    fil = sig.firwin(255, norm_cutoff, pass_zero=True)
    lcf_x = sig.lfilter(fil, 1, x)

    return lcf_x 

def main():
    SAMPLE_FOR = 1 # in seconds
    samplerate, data = scipy.io.wavfile.read(r'Recording.wav')
    time = np.arange(0,SAMPLE_FOR,1/samplerate) #time vector
    data = data[0:int(samplerate*SAMPLE_FOR)]

    BW = 6300
    w = low_cut_filter(data,samplerate,BW)

    samples = SamplerateConversion(w,1260000,44100)
    #samples = generateSignalAM(samples)
    fft_out = fft(samples)
    freq_vector = np.arange(0, 1260000, 1)
    logsn = 20*np.log10(np.abs(fft_out))
    plt.plot(freq_vector, logsn)
    plt.show()

main()

I tested both interpolation and decimation separately and they seem to work fine but not together, I don't know what's causing it.

If it's relevant, my bandwidth is 12.6KHz

UPDATED: changed new_fs to old_fs * L:

Answer 1

I'm assuming this is for an assignment or other educational tasks; as Keegs already pointed out, there are already plenty of ways to resample in Python, in SciPy or other packages (which usually have specific use-cases).

Assuming efficiency is not the driving consideration here, resampling is simple and you have most of the pieces there already: upsampling by interpolation, filtering with an antialias filter, then downsampling.

Doing it efficiently removes the need to store the intermediate upsampled data and filtering that by using a polyphase filter (which does all three steps in one), but that is messy and best left to well-tested packages.

Looking at your code (witout detailed analysis) i notice a few problems:

1. (most importantly) You filter the upsampled signal before downsampling. That is correct. BUT, remember that the sampling frequency of that signal is not new_fs but olf_fs * L!

2. Your upsampling and downsampling functions are really awkward and inefficient. Use slice operators! Upsampling, create a zero array of the needed length, then do out[::L] = in. Downsampling is even easier: out = in[::M].

Edit: I just noticed another thing: Your first plot is in linear scale, your second in log scale (not dB though, if you're using 'abs()', it should be '20*log10()'). So your noise is a good 60 dB below signal, which is probably the quantization noise in your wav file.