Drastically speed up spectrums by zero-padding to a power-of-two.

common.py

@@ -100,3 +100,4 @@ def parabolic_polyfit(f, x, n):
     xv = -0.5 * b/a
     yv = a * xv**2 + b * xv + c
     return (xv, yv)
+

frequency_estimator.py

@@ -2,15 +2,16 @@
 # -*- coding: utf-8 -*-
 
 from __future__ import division
-from common import analyze_channels
-from common import parabolic as parabolic
+from common import analyze_channels, parabolic
 from numpy.fft import rfft
 from numpy import argmax, mean, diff, log, copy, arange
 from matplotlib.mlab import find
 from scipy.signal import fftconvolve, kaiser, decimate
+from scipy.signal import _next_regular # TODO: will be named _next_opt_len() soon. See https://github.com/endolith/waveform-analyzer/pull/6#1
 from time import time
 
 
+
 def freq_from_crossings(signal, fs):
     """Estimate frequency by counting zero crossings
 
@@ -53,7 +54,7 @@ def freq_from_fft(signal, fs):
 
     # Compute Fourier transform of windowed signal
     windowed = signal * kaiser(N, 100)
-    f = rfft(windowed)
+    f = rfft(windowed, _next_regular(N))
     # Find the peak and interpolate to get a more accurate peak
     i_peak = argmax(abs(f)) # Just use this value for less-accurate result
     i_interp = parabolic(log(abs(f)), i_peak)[0]
@@ -103,7 +104,7 @@ def freq_from_hps(signal, fs):
     windowed = signal * kaiser(N, 100)
 
     # Get spectrum
-    X = log(abs(rfft(windowed)))
+    X = log(abs(rfft(windowed, _next_regular(N))))
 
     # Downsample sum logs of spectra instead of multiplying
     hps = copy(X)

thd_analyzer.py

@@ -5,9 +5,11 @@
 from scipy.signal import kaiser
 from numpy.fft import rfft, irfft
 from numpy import argmax, mean, log10, log, ceil, concatenate, zeros
+from scipy.signal import _next_regular # TODO: will be named _next_opt_len() soon. See https://github.com/endolith/waveform-analyzer/pull/6#1
 from common import analyze_channels, rms_flat, parabolic
 from A_weighting import A_weight
 
+
 def THDN(signal, sample_rate):
     """Measure the THD+N for a signal and print the results
 
@@ -32,7 +34,7 @@ def THDN(signal, sample_rate):
     total_rms = rms_flat(windowed)
 
     # Find the peak of the frequency spectrum (fundamental frequency)
-    f = rfft(windowed)
+    f = rfft(windowed, _next_regular(len(windowed)))
     i = argmax(abs(f))
     true_i = parabolic(log(abs(f)), i)[0]
     print 'Frequency: %f Hz' % (sample_rate * (true_i / len(windowed)))
@@ -43,7 +45,7 @@ def THDN(signal, sample_rate):
     f[lowermin: uppermin] = 0
 
     # Transform noise back into the time domain and measure it
-    noise = irfft(f)
+    noise = irfft(f)[:len(windowed)]
     THDN = rms_flat(noise) / total_rms
 
     # TODO: RMS and A-weighting in frequency domain?
@@ -73,7 +75,7 @@ def THD(signal, sample_rate):
     windowed = signal * kaiser(len(signal), 100)
 
     # Find the peak of the frequency spectrum (fundamental frequency)
-    f = rfft(windowed)
+    f = rfft(windowed, _next_regular(len(windowed)))
     i = argmax(abs(f))
     true_i = parabolic(log(abs(f)), i)[0]
     print 'Frequency: %f Hz' % (sample_rate * (true_i / len(windowed)))