There are some Python classes I’m finishing that know how to use the calibration data inside the QA40x so that you can send/receive absolute values to/from the QA40x hardware. This will come in subsequent post. This post is a continuation of the study in the post linked below.
The focus for this first post is on DSP first principles in Python so that you can round-trip a waveform in Python and have it make sense. It’s pretty easy to generate a sine and see the spectrum. But there’s a lot of confusion around how to make noise measurements and amplitude measurements that make sense AND have those measurements make sense when you change FFT sizes. So, for starters here’s a simple framework helping to understand this in Python. I’m not a Python guy, but with ChatGPT I can muddle through. Feel free to point out improvements.
import numpy as np
import matplotlib.pyplot as plt
from scipy.signal.windows import flattop, boxcar, hamming
from scipy.fft import rfft, rfftfreq
# Define the parameters
N = 4096
fs = 48000 # Sampling frequency in Hz
bin_resolution = fs / N
f_target = 1000
f_adjust_strategy = 'snap' # Set to 'none', 'snap', or 'straddle'
# Adjust f_target based on the provided strategy
if f_adjust_strategy == 'snap':
f = round(f_target / bin_resolution) * bin_resolution
elif f_adjust_strategy == 'straddle':
bin_index = round(f_target / bin_resolution)
f = (bin_index + 0.5) * bin_resolution
else:
f = f_target
t = np.arange(N) / fs # Time vector
# Create a sinusoidal signal with a peak amplitude of 1.41 = 0 dBV
amplitude = 1.41
signal = amplitude * np.sin(2 * np.pi * f * t)
# Choose the window function
window_function = 'hamming' # Change to 'flattop', 'boxcar', or 'hamming'
if window_function == 'flattop':
window = flattop(N)
elif window_function == 'boxcar':
window = boxcar(N)
elif window_function == 'hamming':
window = hamming(N)
else:
raise ValueError("Unknown window function")
# Compute the Amplitude Correction Factor (ACF). This factor is used to correct the
# FFT results so that the amplitudes are displayed as RMS (Root Mean Square) values.
# The ACF ensures that the peak amplitude of a 0 dBV sine wave (for example) will
# be accurately represented in the frequency domain plot.
mean_w = np.mean(window)
ACF = 1 / mean_w
# Compute the Energy Correction Factor (ECF). This factor is used to ensure accurate
# energy calculations when integrating across the frequency domain. In the time domain,
# the RMS (Root Mean Square) value of a signal can be calculated directly using the formula:
# RMS = sqrt(mean(signal**2)). However, in the frequency domain, to accurately calculate
# the RMS energy over a range of frequencies, we need to apply the ECF
rms_w = np.sqrt(np.mean(window**2))
ECF = 1 / rms_w
# Print the window function, correction factors, frequencies, and bin size
print(f"Window function: {window_function}")
print(f"Amplitude Correction Factor (ACF): {ACF:.4f} ({20 * np.log10(ACF):.2f} dB)")
print(f"Energy Correction Factor (ECF): {ECF:.4f} ({20 * np.log10(ECF):.2f} dB)")
print(f"f_target: {f_target} Hz")
print(f"f (adjusted): {f} Hz")
print(f"Bin size: {bin_resolution:.2f} Hz")
print(f"Frequency adjustment strategy: {f_adjust_strategy}")
# Apply the selected window
windowed_signal = signal * window
# Perform the FFT and normalize. Normalizing means that, regardless of the FFT size,
# we'll still arrive at the same amplitude for a given time-domain signal Note the rfft
# assumes real input data, and computes only the positive freq components
fft_result = rfft(windowed_signal)
fft_normalized = (np.abs(fft_result) / (N / 2)) / np.sqrt(2)
# Apply the ACF to the entire FFT magnitude. With this step, a 0 dBV sine will have
# a peak of 0 dBV
fft_magnitude_acf = fft_normalized * ACF
fft_magnitude_acf_dBV = 20 * np.log10(fft_magnitude_acf) # Convert to dBV
# Compute the frequencies corresponding to the positive half of the FFT
frequencies = rfftfreq(N, d=1/fs)
# Find the peak in the corrected FFT magnitude
peak_index = np.argmax(fft_magnitude_acf) # Consider only the positive frequencies
plot_peak_frequency = peak_index * bin_resolution
peak_amplitude = fft_magnitude_acf[peak_index]
peak_amplitude_dBV = 20 * np.log10(peak_amplitude)
# Compute the RMS energy of the entire spectrum (excluding DC and Nyquist bins)
rms_energy_entire_spectrum = np.sqrt(np.sum(fft_normalized[1:]**2)) * ECF
rms_energy_entire_spectrum_dBV = 20 * np.log10(rms_energy_entire_spectrum)
# Compute the RMS energy from start frequency to stop frequency
start_freq = 900
stop_freq = 1100
# Find the indices corresponding to start and stop frequencies
start_index = np.searchsorted(frequencies, start_freq, side='left')
stop_index = np.searchsorted(frequencies, stop_freq, side='right')
# Compute the RMS energy within the specified frequency range
rms_energy_range = np.sqrt(np.sum(fft_normalized[start_index:stop_index]**2)) * ECF
rms_energy_range_dBV = 20 * np.log10(rms_energy_range)
# Print frequency domain peak values
print(f"Plot Peak Frequency: {plot_peak_frequency:.2f} Hz")
print(f"Plot Peak Amplitude: {peak_amplitude:.4f} ({peak_amplitude_dBV:.2f} dBV)")
print(f"RMS Energy of the Entire Spectrum (excluding DC and Nyquist): {rms_energy_entire_spectrum:.4f} ({rms_energy_entire_spectrum_dBV:.2f} dBV)")
print(f"RMS Energy from {start_freq} Hz to {stop_freq} Hz: {rms_energy_range:.4f} ({rms_energy_range_dBV:.2f} dBV)")
# Plot the window function, signal, and its FFT
plt.figure(figsize=(12, 6))
# Plot the window function
plt.subplot(2, 2, 1)
plt.plot(t, window, label='Window Function')
plt.xlabel('Time [s]')
plt.ylabel('Amplitude')
plt.title('Window Function')
plt.legend()
plt.grid()
# Plot the time domain signal
plt.subplot(2, 2, 2)
plt.plot(t, signal, label='Original Signal')
plt.plot(t, windowed_signal, label='Windowed Signal')
plt.xlabel('Time [s]')
plt.ylabel('Amplitude')
plt.title('Time Domain Signal')
plt.legend()
plt.grid()
# Plot the FFT magnitude
plt.subplot(2, 2, 3)
ax3 = plt.gca()
ax3.plot(frequencies, fft_magnitude_acf, label='Magnitude', color='b')
ax3.set_xlabel('Frequency [Hz]')
ax3.set_ylabel('Magnitude')
ax3.set_title('Frequency Domain (FFT)')
ax3.grid()
# Add a secondary y-axis to show dBV
ax3_right = ax3.twinx()
ax3_right.plot(frequencies, fft_magnitude_acf_dBV, color='r', alpha=0.5, linestyle='--', label='Magnitude (dBV)')
ax3_right.set_ylabel('Magnitude (dBV)')
ax3_right.grid(False)
# Combine legends
lines, labels = ax3.get_legend_handles_labels()
lines2, labels2 = ax3_right.get_legend_handles_labels()
ax3_right.legend(lines + lines2, labels + labels2, loc='upper right')
# Plot a zoomed-in view of +/- 10 bins from the center frequency
plt.subplot(2, 2, 4)
ax4 = plt.gca()
zoom_start = max(0, peak_index - 10)
zoom_end = min(len(frequencies), peak_index + 11)
ax4.plot(frequencies[zoom_start:zoom_end], fft_magnitude_acf[zoom_start:zoom_end], label='Magnitude', color='b')
ax4.set_xlabel('Frequency [Hz]')
ax4.set_ylabel('Magnitude')
ax4.set_title(f'Zoomed-in View (+/- 10 bins) around {plot_peak_frequency:.2f} Hz')
ax4.grid()
# Add a secondary y-axis to show dBV
ax4_right = ax4.twinx()
ax4_right.plot(frequencies[zoom_start:zoom_end], fft_magnitude_acf_dBV[zoom_start:zoom_end], color='r', alpha=0.5, linestyle='--', label='Magnitude (dBV)')
ax4_right.set_ylabel('Magnitude (dBV)')
ax4_right.grid(False)
# Combine legends
lines, labels = ax4.get_legend_handles_labels()
lines2, labels2 = ax4_right.get_legend_handles_labels()
ax4_right.legend(lines + lines2, labels + labels2, loc='upper right')
plt.tight_layout()
plt.show()
The code above can be run in JupyterLab, and the following output will show:
Note in the above we’re using a hamming window with a frequency adjustment strategy of “snap”. There are 3 frequency adjustment strategies in the code:
none When this is selected, the f_target you select will be the frequency that is used for generating a sine.
snap When this is selected, the f_target you select will be nudged slighted such that the frequency falls in the center of an fft bin. This is ideal for signal processing, because it usually means you don’t need to use a window.
straddle When this is selected, the f_target will be nudged so that the frequency falls between two bins. This is the often the worst-case for dealing with amplitude issues related to windowing. We’ll see this in more detail below.
In the first example we ran, a snap strategy was used, which means the frequency was adjusted from our target of 1000 Hz to 996.09375 Hz. This adjustment put the tone right in the middle of a bin. Note in the Zoomed-in View how the peak is precisely 1.0 Vrms (our source sine was a 1.41V peak, which is 1.0Vrms. So, our frequency domain graphs are accurately reflecting the rms of the sine waves.
Now we can change the Window to boxcar (rectangle) and see how the output changes:
By and large, it looks very similar. And this is expected, because when your signal falls in the middle of a bin, your amplitude errors are minimized.
Now, let’s change the adjustment strategy to ```straddle`` and run it again:
Hmmm. Now we’ve got a problem. The peak shown in the frequency domain looks very low. And instead of a single peak in the zoomed-in view, we now have two peaks. Since we asked for a tone to be placed between two bins, this makes sense. Each bin will have half the energy. Finally, notice the skirts. At 1050 Hz, we can the amplitude is quite high–around -23 dBV.
Now let’s switch the window to flattop:
This is better! Notice the amplitudes in the freq plots are back to 0 dBV and the skirts are also improved.
RMS Energy Measurements
The trials above include RMS energy measurements. These are measurements made on the spectral data. It’s well known how to compute RMS from the time domain series. But it’s very useful to compute RMS from a spectrum too. Notice that the RMS measurements are much less sensitive to window selection. When boxcar/straddle is used, there’s a bit of error in the RMS energy calc due to the skirts being high outside the of the 900 to 1100 Hz. But that could be corrected by opening the window a bit. When the entire spectrum is considered, the RMS is correct regardless of tone frequency and/or window type.
ACF and ECF
In the code, note that we have two adjustments that need to be made to the raw spectrum data that is returned by the rfft function. The first is Amplitude Correction Factor, or ACF. The second is Energy Correction Factor or ECF.
A way to think about this is as follows: The ACF must be applied to data you want to plot, and the ECF must be applied to data you want to integrate (to learn RMS, for example).
Thankfully, both can be learned from the window shape with some simple math.
# Compute the Amplitude Correction Factor (ACF). This factor is used to correct the
# FFT results so that the amplitudes are displayed as RMS (Root Mean Square) values.
# The ACF ensures that the peak amplitude of a 0 dBV sine wave (for example) will
# be accurately represented in the frequency domain plot.
mean_w = np.mean(window)
ACF = 1 / mean_w
# Compute the Energy Correction Factor (ECF). This factor is used to ensure accurate
# energy calculations when integrating across the frequency domain. In the time domain,
# the RMS (Root Mean Square) value of a signal can be calculated directly using the formula:
# RMS = sqrt(mean(signal**2)). However, in the frequency domain, to accurately calculate
# the RMS energy over a range of frequencies, we need to apply the ECF
rms_w = np.sqrt(np.mean(window**2))
ECF = 1 / rms_w
Adding Another Tone
We can tweak the code above to add another tone that is 3X the adjusted frequency at an amplitude that is 1/100th the indicated amplitude (-40 dB below the fundamental).
# Create a sinusoidal signal with a peak amplitude of 1.41 = 0 dBV
amplitude = 1.41
signal = amplitude * np.sin(2 * np.pi * f * t)
# Add a secondary frequency component at 3 times the target frequency with 1/100 the amplitude
secondary_frequency = 3 * f
secondary_amplitude = amplitude / 100
signal += secondary_amplitude * np.sin(2 * np.pi * secondary_frequency * t)
We can see the added tone appear in the frequency domain at the -40 dBV level:
Summary
It can be very challenging to get FFT data that matches expectations as you change FFT size, window types and signal frequency. The code above makes it easy to adjust key parameters and see how the results match your expectations.
Next, we’ll pull in the some classes that can handle the various calibrations required to get absolute data in and out of the QA40x hardware.
I wanted to reach out and ask about the calibration, I took a look last time we spoke, but didn’t get too far!
