0
$\begingroup$

I'm performing a STFT of a sine wave in which the frequency increases (up-chirp signal) using python numpy.fft as shown in the following code:

import numpy as np
import matplotlib.pyplot as plt

def stft(x, L, overlap, window):
    opoints = int(L*overlap)
    X = np.array([np.abs(np.fft.fft(window*x[i:i+L])[:L//2])
                 for i in range(0,len(x)-L,opoints)])
    return X

t = np.linspace(0, 5, 11025*5)

xf = np.linspace(0, 2000, 11025*5)
w = 2*np.pi*xf
x = np.sin(w*t)

f = np.linspace(0, 11025, 1024)

window = np.hanning(1024)

stft = stft(x,1024,0.5,window)

plt.imshow(stft[:,:1024//2].T, cmap = plt.get_cmap('gist_yarg'), origin='lower', aspect='auto',
                 interpolation='nearest', extent=(0,5, 0,11025/2))

plt.xlabel("Time [s]")
plt.ylabel("Frequency [Hz]")
plt.show()

It works for a constant frequency sinewave. However it shows not an offset but a gain error if using the up-chirp signal as input.

STFT of increasing frequency sine wave

I'm using:

  • Hanning window
  • Overlap = 0.5
  • Window size 1024
  • Sampling frequency 11025Hz
  • Total sample time: 5 seconds

What is causing the gain error that leads the 2 kHz signal appearing at double frequency at the end of the time window?

Thank you in advance, Pedro

$\endgroup$

1 Answer 1

1
$\begingroup$

A sinusoidal signal is

$$y(t) = sin(\phi(t))$$

The phase of a linear chirp, is

$$\phi(t) = \phi_0 + 2\pi\left(f_0t+\frac{k}{2}t^2\right)$$ where $$ k = \frac{f_1-f_0}{t_1-t_0}$$

So your chirp signal is missing a 0.5 correction for the frequency variation.

Use the scipy.signal.chirp:

import numpy as np
import matplotlib.pyplot as plt
from scipy.signal import chirp

def stft(x, L, overlap, window):
    opoints = int(L*overlap)
    X = np.array([np.abs(np.fft.fft(window*x[i:i+L])[:L//2])
                 for i in range(0,len(x)-L,opoints)])
    return X

tmax = 5
fs = 11025
N = fs*tmax

t = np.linspace(0, tmax, N)
x = chirp(t,0,tmax,2000)
windowsize = 1024

f = np.linspace(0, fs, windowsize)
window = np.hanning(windowsize)
stft = stft(x,windowsize,0.5,window)

plt.imshow(stft.T, cmap = plt.get_cmap('gist_yarg'), origin='lower', aspect='auto',
           interpolation='nearest', extent=(0,t[-1], 0,fs/2))
plt.xlabel("Time [s]")
plt.ylabel("Frequency [Hz]")
plt.show()

Chirp evolution

Sources:

1: https://en.wikipedia.org/wiki/Chirp

2: Basic question about chirp signal

$\endgroup$
2
  • $\begingroup$ You are right, the bug was causing the gain error not to show when using overlap < 0.05 since all the signal data was not being displayed. I edited the first message. Thank you for the contribution, still trying to figure out what is causing the gain error. $\endgroup$ Oct 24, 2017 at 14:18
  • 1
    $\begingroup$ Thank you so much @daguiam, again you are right, I was not understanding the chirp signal properly. $\endgroup$ Oct 24, 2017 at 22:55

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

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