Search code examples
pythonnumpysignalssignal-processingfrequency

Interpolating between two frequencies in numpy


I want to create a sine wave that starts from the frequency f1 and ends at the frequency f2. Here is the code I used:

import matplotlib.pyplot as plt
import numpy as np

def freq_interp(dur,f1,f2,fs=44100):
    num_samples = fs*dur
    t = np.linspace(0,dur,num_samples)
    a = np.linspace(0,1,num_samples)
    f = (1-a)*f1+a*f2 # interpolate
    samples = np.cos(2*np.pi*f*t)
    return samples,f

When I try to generate a WAV file or just plot the STFT of the signal, I get an unexpected result. For example I used the code below:

def plot_stft(sig,fs=44100):
    f, t, Zxx = signal.stft(sig,fs=fs,nperseg=2000)
    plt.pcolormesh(t, f, np.abs(Zxx), vmin=0, vmax=0.1)
    plt.ylim(0,2000)
    plt.title('STFT Magnitude')
    plt.ylabel('Frequency [Hz]')
    plt.xlabel('Time [sec]')
    plt.show()

s,f = freq_interp(dur=2,f1=1,f2=1000)
plt.plot(f)
plt.show()
plot_stft(s)

s,f = freq_interp(dur=2,f1=1000,f2=1)
plt.plot(f)
plt.show()
plot_stft(s)

I get these plots: enter image description here

The problem is more evident in the second row. Where the frequency has bounced back at t=1s. Also in the first row you can see that the frequency has gone up to 2000Hz which is wrong. Any idea why this happens and how I can fix it?


Solution

  • A sin wave is sin(p(t)) where p(t) is the phase function. And frequency function is f(t) = d p(t) / dt, to calculate p(t), you first calculate f(t) and then integrate it. The simplest method of integration is by using cumsum().

    def freq_interp(dur,f1,f2,fs=44100):
        num_samples = int(fs*dur)
        t = np.linspace(0,dur,num_samples)
        f = np.linspace(f1, f2, num_samples)
        phase = 2 * np.pi * np.cumsum(f) / fs
        samples = np.cos(phase)
        return t, samples