# How to change the samplimg rate with interpolation and decimation?

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
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: • Can you attach a sample of the test signal that causes this? Also that plot would be a lot more useful in dB if you wouldn't mind. Sep 26, 2021 at 6:29
• Done, I added the file and the rest of my code to the question @Keegs Sep 26, 2021 at 21:29
• I appreciate it, but I need the low_cut_filter code too, and that file you posted is an mp3... Sep 27, 2021 at 0:39
• Sorry, that site converted it. I reuploaded it to another site and added the low_cut_filter code @Keegs Sep 27, 2021 at 10:22
• Ok, so are you trying to do this as some sort of experiment or assignment? Because if all you want is a shortcut to something that does this for you, you can use scipy.signal.upfirdn docs.scipy.org/doc/scipy/reference/generated/… Sep 28, 2021 at 7:07

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.

• Thank you for the answer, it is for educational purposes only and I'm trying to keep it as simple as I can that's why I don't use slicing. You are right about changing the Fs from new_fs to old_fs * L however it still didn't solve my problem, but made it smaller, it seems like there's only one extra spike now (see updated post) Sep 28, 2021 at 11:19