Validating W&F analysis using MI as reference

[aoliver]

Hello,
First of all, thank you for this outstanding software. I purchased the full package a few days ago and I’m absolutely thrilled with the range of features it offers.

I’m developing a Python script to analyze unweighted wow and flutter, using the methodology described in the “Wow and Flutter Measurement using Multi-Instrument” guide and the MI user manual as a reference. I'm trying to validate the script’s output by comparing it against measurements obtained from Multi-Instrument, but I’m getting inconsistent results.

The script is fairly straightforward and self-explanatory. These are the main steps I'm following:

1. Apply a full-length Kaiser 8 window to the audio signal

2. Perform a bandpass prefilter (312.9 Hz – 5945 Hz)

3. Extract the instantaneous frequency via Hilbert demodulation

4. Normalize and remove the DC component

5. Apply a 1565 Hz anti-alias low-pass filter

6. Compute the absolute value of the deviation (abs)

7. Calculate two-sigma (2 × standard deviation) and 95th percentile

As you'll see, at first everything seems to work well. When analyzing a synthetic signal modulated at 4 Hz with a deviation of 0.997%, the 95th percentile value matches exactly, while—just as expected—the 2-sigma value is invalid due to the signal's non-Gaussian distribution.




However, when using a real signal recorded from my turntable, MI reports 0.45%, while the script gives 0.54% using the 2-sigma method.




Why do I get such different readings? Am I missing an additional processing step?
Any guidance pointing me in the right direction would be greatly appreciated.

Best regards,

Alvaro Oliver

Code: Select all

import numpy as np
from scipy.io import wavfile
from scipy.signal import butter, filtfilt, hilbert, get_window
import matplotlib.pyplot as plt

# -------------------- CONFIG --------------------
#FILE_NAME = "Record1.wav"
FILE_NAME = "04 3150HzModulatedBy4HzAt0.997Percent31.5Hz.wav"
BPF_HI = 5945.0     # bandpass high-pass frequency
BPF_LO = 312.9      # bandpass low-cut frequency
LPF_CUTOFF = 1565.0 # FM demodulation antialias filter (1565.0)
TRIM_MS = 20        # trim 20 ms each side
BETA = 8.0          # Kaiser window parameter

# -------------------- FUNCTIONS --------------------
def bandpass_filter(signal, fs, lowcut=BPF_LO, highcut=BPF_HI, order=4):
    nyq = 0.5 * fs
    b, a = butter(order, [lowcut / nyq, highcut / nyq], btype='band')
    return filtfilt(b, a, signal)

def lowpass_filter(signal, fs, cutoff, order=4):
    nyq = 0.5 * fs
    b, a = butter(order, cutoff / nyq, btype='low')
    return filtfilt(b, a, signal)

# -------------------- MAIN --------------------
# Load WAV audio
fs, data = wavfile.read(FILE_NAME)
TRIM = int(fs * TRIM_MS / 1000)
audio = data.astype(np.float32)
if audio.ndim > 1:
    audio = audio[:, 0]  # just left chanel

# Fit window to the full length of the audio
audio_len = len(audio)
window = get_window(('kaiser', BETA), audio_len)
block = audio * window

# Instantaneous Frequency
prefiltered = bandpass_filter(block, fs)
analytic = hilbert(prefiltered)
phase = np.unwrap(np.angle(analytic))
inst_freq = np.diff(phase) * fs / (2 * np.pi)
inst_freq = inst_freq[TRIM:-TRIM]

# Normalization
mean_freq = np.mean(inst_freq)
normalized = (inst_freq - mean_freq) / mean_freq # normalization and DC removal
dealiased = lowpass_filter(normalized, fs, LPF_CUTOFF) # antialias postfilter

# Unweighted Wow and Flutter
magnitude = np.abs(dealiased) # deviation magnitude
sigma_2 = 2 * np.std(magnitude) # unweighted peak WnF (2σ)
percentil_95 = np.percentile(magnitude, 95) # percentile 95 peak

t = np.arange(len(dealiased)) / fs # time axis

# -------------------- PLOT RESULTS --------------------
fig = plt.figure(figsize=(14, 10))
fig.suptitle(FILE_NAME, fontsize=12)

# Subplot Layout
ax1 = plt.subplot2grid((3, 2), (0, 0), colspan=2)
ax2 = plt.subplot2grid((3, 2), (1, 0), colspan=2)
ax3 = plt.subplot2grid((3, 2), (2, 0))
ax4 = plt.subplot2grid((3, 2), (2, 1))

# Plot 1: Instantaneous Frequency
ti = np.arange(len(inst_freq)) / fs
ax1.plot(ti, inst_freq, linewidth=0.3)
ax1.axhline(y=mean_freq, color='red', linestyle='--', label=f"Mean: {mean_freq:.2f} Hz", linewidth=1)
ax1.set_xlabel("Time (s)")
ax1.set_title("Instantaneous Frecuency (Hz)")
ax1.legend()
ax1.grid()

# Plot 2: Magnitude
ax2.plot(t, magnitude, linewidth=0.3)
ax2.axhline(y=sigma_2, color='red', linestyle='--', label=f"Unw. WnF (2σ): {sigma_2 * 100:.4f}% peak", linewidth=1)
ax2.axhline(y=percentil_95, color='blue', linestyle='--', label=f"Unw. WnF (P95): {percentil_95 * 100:.4f}% peak", linewidth=1)
ax2.set_title("Normalized Absolute Deviation Magnitude")
ax2.set_xlabel("Time (s)")
ax2.legend()
ax2.grid()

# Plot 3: Spectrum
spectrum = np.abs(np.fft.rfft(dealiased))
freqs = np.fft.rfftfreq(len(dealiased), d=1/fs)
ax3.semilogx(freqs, spectrum / np.max(spectrum), color='purple', alpha=0.7) # eje x logarítmico
ax3.set_title("Spectrum")
ax3.set_xlabel("Frequency (Hz)")
ax3.grid()
ax3.grid(which="minor", color="0.9")
ax3.set_xlim([freqs[1], 1000])

# Plot 4: Histogram
ax4.hist(dealiased, bins=1024, density=True, color='gray', alpha=0.6)
ax4.set_title("Histogram")
ax4.set_xlim(-0.05,0.05)
ax4.grid()

plt.tight_layout(rect=[0, 0, 1, 0.96]) # leave space for filename
plt.show()


[VirtinsTech]

Thank you for your questions.

In the wow and flutter measurement of Multi-Instrument, a window function (e.g., Kaiser 8) is used only to determine the average frequency of the test signal in the frequency domain. It is not used in FM demodulation. Alternatively, the average frequency of the test signal can be obtained using frequency counting in the time domain.

Only the 95th percentile is used, as specified by the AES standard. This method is also known as the 2-sigma method.

[aoliver]

Thank you very much for your helpful response.

Following your guidance, I’ve removed the Kaiser 8 window from my Python implementation, as it is not part of the FM demodulation process. As expected, this change had no noticeable effect on the results, which remain consistent with the previous version of the script.

May I ask a follow-up question:
Is this signal processing chain (bandpass filtering → Hilbert demodulation → DC removal → anti-alias LPF → absolute deviation → P95) equivalent to the one used by the Unweighted Peak WnF DDP in Multi-Instrument? Or does the DDP apply additional steps I may be missing?

Assuming the processing is correct and complete, what could explain the difference in results between the script and Multi-Instrument when analyzing real-world recordings? On synthetic signals, the results match perfectly, but with real audio, MI typically reports a lower W&F (P95) value than my script.

Thank you again for your support. I'm very interested in understanding the underlying methodology more deeply and ensuring that my implementation aligns as closely as possible with yours.

[VirtinsTech]

Yes, MI basically uses the same signal processing chain. However, there may be some differences in the detailed implementation, such as those in the following aspects:
(1) the cutoff frequencies of the bandpass filter and low pass filter. In MI, they are not fixed but based on the measured carrier frequency.
(2) the digital filtering methods
(3) the suppression / removal methods for boundary effects resulting from filtering and FM demodulation.
(4) the differencing method used in FM demodulation. This is particular important when the ratio of sampling rate to carrier frequency is low.

[aoliver]

Thank you very much for the responses. I completely understand that you cannot share more details about the implementation of MI. Regards.