0
$\begingroup$

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.

enter image description here

enter image description here

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:

enter image description here

$\endgroup$
6
  • $\begingroup$ 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. $\endgroup$
    – Keegs
    Sep 26 '21 at 6:29
  • $\begingroup$ Done, I added the file and the rest of my code to the question @Keegs $\endgroup$ Sep 26 '21 at 21:29
  • $\begingroup$ I appreciate it, but I need the low_cut_filter code too, and that file you posted is an mp3... $\endgroup$
    – Keegs
    Sep 27 '21 at 0:39
  • $\begingroup$ Sorry, that site converted it. I reuploaded it to another site and added the low_cut_filter code @Keegs $\endgroup$ Sep 27 '21 at 10:22
  • 1
    $\begingroup$ 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/… $\endgroup$
    – Keegs
    Sep 28 '21 at 7:07
1
$\begingroup$

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.

$\endgroup$
1
  • $\begingroup$ 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) $\endgroup$ Sep 28 '21 at 11:19

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.