Here it's about calculating the Wigner-Ville-Distribution (WVD) with Python. I found no Python-libary for this, so I ported it from the Time-Frequency Toolbox for GNU Octave and Matlab. I added calculating the analytic signal (for avoiding interferences between the negative and positive frequency components) and a filtering of the WVD by multiplying it with the short-time Fourier transform (STFT). It is only about the auto-WVD, beacause – honestly – I have no idea how to apply the cross-WVD to audio in a meaningful way. If you have some ideas or any other comments or improvements, please write me!
The WVD is a time-frequency-transformation with the finest resolution, but it has the disadvantage of interferences. Furthermore one has to pay the high-resolution with calculation time. Half a second in CD-quality induces a matrix of 22050 × 22050 ≈ 486 million values. Below you find the classic example – the superposition of two opposed chirp signals – by comparison: the STFT, WVD and filtered WVD.
Short-time Fourier transform:
The code for the STFT and for plotting is to be found here.
#coding: utf-8 """ Filtered and non-filtered Wigner Ville Distribution This work is licensed under a Creative Commons Attribution 3.0 Unported License. Frank Zalkow, 2013 ported from tfrwv from the matlab/octave Time-Frequency Toolbox (http://tftb.nongnu.org/) differences to tfrwv: - only takes one signal -> only the auto wigner distribution, no cross wigner distribution - argument make_analytic for calculating the wigner distribut- ion of the analytic signal (less crossterms because negative frequency content is zero) """ from time import clock from numpy import abs, arange, shape, array, ceil, zeros, conj, ix_, transpose, append, fft, real, float64, linspace, sqrt from scipy.signal import hilbert from scipy import interpolate from math import log, ceil, floor """ get next power of 2 """ def nextpow2(n): return 2 ** ceil(log(n, 2)) """ auxilliary function for wvd if trace is true """ def disprog(i, N, steps): global begin_time_disprog if i == 0: begin_time_disprog = clock() if i == (N-1): print "100 %% complete in %f seconds." % (clock() - begin_time_disprog) del begin_time_disprog elif (floor(i * steps / float(N)) != floor((i-1) * steps / float(N))) : print "%d" % (floor(i * steps / float(N)) * ceil(100.0 / float(steps))), """ calculate the wigner ville distribution of an audio file """ def wvd(audioFile, t=None, N=None, trace=0, make_analytic=True): if make_analytic: x = hilbert(audioFile) else: x = array(audioFile) if x.ndim == 1: [xrow, xcol] = shape(array([x])) else: raise ValueError("Signal x must be one-dimensional.") if t is None: t = arange(len(x)) if N is None: N = len(x) if (N <= 0 ): raise ValueError("Number of Frequency bins N must be greater than zero.") if t.ndim == 1: [trow, tcol] = shape(array([t])) else: raise ValueError("Time indices t must be one-dimensional.") if xrow != 1: raise ValueError("Signal x must have one row.") elif trow != 1: raise ValueError("Time indicies t must have one row.") elif nextpow2(N) != N: print "For a faster computation, number of Frequency bins N should be a power of two." tfr = zeros([N, tcol], dtype='complex') if trace: print "Wigner-Ville distribution", for icol in xrange(0, tcol): ti = t[icol] taumax = min([ti, xcol-ti-1, int(round(N/2.0))-1]) tau = arange(-taumax, taumax+1) indices = ((N+tau)%N) tfr[ix_(indices, [icol])] = transpose(array(x[ti+tau] * conj(x[ti-tau]), ndmin=2)) tau=int(round(N/2))+1 if ((ti+1) <= (xcol-tau)) and ((ti+1) >= (tau+1)): if(tau >= tfr.shape): tfr = append(tfr, zeros([1, tcol]), axis=0) tfr[ix_([tau], [icol])] = array(0.5 * (x[ti+tau] * conj(x[ti-tau]) + x[ti-tau] * conj(x[ti+tau]))) if trace: disprog(icol, tcol, 10) tfr = real(fft.fft(tfr, axis=0)) f = 0.5*arange(N)/float(N) return (transpose(tfr), t, f ) """ get the filtered wvd by multiplying the wvd and the stft """ def filtered_wvd(wvd, stft): qstft = abs(stft) qstft = float64(qstft * qstft) bigstft = zeros(shape(wvd), float64) x = arange(0, shape(qstft)) y = arange(0, shape(qstft)) xx = linspace(x.min(), x.max(), shape(wvd)) yy = linspace(y.min(), y.max(), shape(wvd)) interpolator = interpolate.RectBivariateSpline(x,y,qstft, kx=1,ky=1) bigstft = interpolator(xx,yy) return (sqrt(abs(bigstft * wvd)), wvd, wvd)