The FFT of a sine wave with varying frequency and constant amplitude does not show a plateau, why?

Question

I want to use python to see what the audio response of a speaker is. This would be done by comparing the output of the speaker measured with a microphone to the given input. That means that for both of these signals that will be in the time domain I need to perform a FFT. Eventually I want to play a target curve close to a Harman curve. I tried to generate the signal to put it in a wav file and when checking if the signal does what I ask with the FFT I got some weird behaviour. So then I tried to perform an FFT on a sine wave that has a constant amplitude and goes from 21-19700 Hz in 5 seconds with a sample rate of 44100 Hz. I would like to use the welch method where I took the window to be the full length of the sample as in this case the signal I generated does not have any noise. This got me the following result:

I used the following python code

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

samplerate = 44100
duration = 5
time = np.linspace(0,duration,samplerate*duration, endpoint= False)
frequency = np.linspace(21,19700,len(time),endpoint= False)
Const_wave = np.sin(np.pi*2*frequency*time)
f, check = welch(Const_wave,samplerate,nperseg=len(Wave))
plt.plot(f,check)
plt.xscale('log')
plt.xlabel('Frequency [Hz]')
plt.yscale('log')
plt.ylabel('Power')
plt.grid()
plt.show()

Afterwards, I tried the scipy.fft module which yielded a different but also not expected result:

for this I used the following code:

f = rfftfreq(samplerate*duration,1/samplerate)
check = rfft(Const_wave)
plt.plot(f,check) 
plt.xscale('log')
plt.xlabel('Frequency [Hz]')
plt.yscale('log')
plt.ylabel('Power')
plt.grid()
plt.show()

I would expect to see a flat plateau of equal power on all frequencies, but that is not happening. I also tried to do the same with a linear sweep wav file generated by REW and got similar results, so is my expectation of a flat plateau incorrect or is there something else at play here?

Answer 1

Your Python code has two weak spots, rather typos, but these are not the source of your failure to "show a plateau" in PSD graph. First we review these "typos".

The line

f, check = welch(Const_wave,samplerate,nperseg=len(Wave))

must read

f, check = welch(Const_wave,samplerate,nperseg=len(Const_wave))

or else, the variable Const_wave must be Wave. Without correcting this "typo", which might appear when you copied the code to the post, your code does not compile.

Then,

time = np.linspace(0,duration,samplerate*duration, endpoint= False)
frequency = np.linspace(21,19700,len(time),endpoint= False)
Const_wave = np.sin(np.pi*2*frequency*time)

results in a frequency sweep up to 19700*2, not to 19700, but maybe you did it intentionally. If not, notice that the instantaneous frequency is the derivative of the phase. To reach a frequency of $f_{max}$ in the frequency sweep, the formulas are

$$ phase(t) = f_{min}·t + (f_{max} - f_{min})/duration·t^2/2 \\ sweep\_wave(t) = \sin(2π·phase(t)) \\ instantaneous\_freq(t) = f_{min} + ((f_{max} - f_{min})/duration)·t $$

and you should write the code

from scipy import signal 
...
fmin = 21
fmax = 19700
time = np.linspace(0,duration,samplerate*duration, endpoint= False)
Const_wave = np.sin(np.pi*2*(fmin*time + (fmax-fmin)/duration*time*time/2.0))

---

Back to your question why "a sine wave with varying frequency and constant amplitude does not show a plateau": the answer can be found in your choice of the Welch method's parameters. You apply the Welch method to the signal in the entire range processed with a Hann window (see the scipy.signal.welch reference), and your Hann window "leaks" into the PSD. You have two options to reach the desired "plateau".

1st: Specify the named parameter window with the array-like value, like this

f, check = signal.welch(Const_wave,samplerate,window=np.ones(len(Const_wave)))

2nd: Retaining the default (Hann) window (omitting the named parameter window), decrease the value of the named parameter nperseg to a lower value commensurate with a maximum wave period (minimum frequency) in the sweep, like this

f, check = signal.welch(Const_wave,samplerate,nperseg=samplerate/frequency[0])

you can also try to adjust the noverlap parameter value (defaulted to half the nperseg value).

In both cases, the graphs would not yet become ideal "plateaus", and to understand discrepancies--quite substantial in the second case--you should delve deeper into the Welch's method working in particular and into the study of signal processing algorithms in general.

Any way, you have to, since the audio engineering requires a quite solid background in signal processing. In the future, your involvement with more advanced AE projects may (and often will) require a subscription to the Journal of the Audio Engineering Society.

In the meantime, let me recommend you references related to the issues you raised in your question (including Harman target curves):

Transfer-Function Measurement with Sweeps by Swen Müller and Paulo Massarani

Simultaneous measurement of impulse response and distortion with a swept-sine technique by Angelo Farina

Github's project AutoEQ

Answer 2

This is a common trap in generating frequency ramps. For the ramp given by $\cos(2\pi f(t) t)$ with the ramp function as $f(t)$, $f(t)$ is NOT the instantaneous frequency desired. Please see this post which fully details this:

Simulation of a Frequency ramp

For optimized frequency chirp generation for channel estimation that accounts for this as well as minimizes FFT related aliasing effects, please see this post:

How can I plot the frequency response on a bode diagram with Fast Fourier Transform?