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
$\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$
  • $\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$ – Pedro Martinez Lopez Oct 24 '17 at 14:18
  • 1
    $\begingroup$ Thank you so much @daguiam, again you are right, I was not understanding the chirp signal properly. $\endgroup$ – Pedro Martinez Lopez Oct 24 '17 at 22:55

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.