1function W=ref_wignervilledist(f,g)
2%-*- texinfo -*-
3%@deftypefn {Function} ref_wignervilledist
4%@verbatim
5%REF_WIGNERVILLEDIST Reference wigner-Ville distribution
6%   Usage:  W=ref_wignervilledist(f)
7%
8%   REF_WIGNERVILLEDIST(f,g) computes the Wigner-Ville distribution of f and g.
9%@end verbatim
10%@strong{Url}: @url{http://ltfat.github.io/doc/reference/ref_wignervilledist.html}
11%@end deftypefn
12
13% Copyright (C) 2005-2016 Peter L. Soendergaard <peter@sonderport.dk>.
14% This file is part of LTFAT version 2.3.1
15%
16% This program is free software: you can redistribute it and/or modify
17% it under the terms of the GNU General Public License as published by
18% the Free Software Foundation, either version 3 of the License, or
19% (at your option) any later version.
20%
21% This program is distributed in the hope that it will be useful,
22% but WITHOUT ANY WARRANTY; without even the implied warranty of
23% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
24% GNU General Public License for more details.
25%
26% You should have received a copy of the GNU General Public License
27% along with this program.  If not, see <http://www.gnu.org/licenses/>.
28
29% AUTHOR: Jordy van Velthoven
30
31if ~all(length(f)==length(g))
32  error('%s: f and g must have the same length.', upper(mfilename));
33end;
34
35L = length(f);
36H = floor(L/2);
37R = zeros(L,L);
38W = zeros(L,L);
39
40% Compute the analytic representation of f
41if (nargin == 1)
42  if isreal(f)
43   z = fft(f);
44   z(2:L-H) = 2*z(2:L-H);
45   z(H+2:L) = 0;
46   z1 = ifft(z);
47   z2 = z1;
48  else
49   z1 = f;
50   z2 = z1;
51  end
52elseif (nargin == 2)
53  if isreal(f) || isreal(g)
54   z1 = fft(f);
55   z1(2:L-H) = 2*z1(2:L-H);
56   z1(H+2:L) = 0;
57   z1 = ifft(z1);
58   z2 = fft(g);
59   z2(2:L-H) = 2*z2(2:L-H);
60   z2(H+2:L) = 0;
61   z2 = ifft(z2);
62  else
63    z1 = f;
64    z2 = g;
65  end
66end
67
68% Compute instantaneous autocorrelation matrix R
69for l = 0 : L-1;
70  for m = -min([L-l, l, round(L/2)-1]) : min([L-l, l, round(L/2)-1]);
71    R(mod(L+m,L)+1, l+1) =  z1(mod(l+m, L)+1).*conj(z2(mod(l-m, L)+1));
72  end
73end
74
75% Compute the Wigner-Ville distribution
76for hh=0:L-1
77  for ii=0:L-1
78    for jj = 0:L-1
79      W(hh+1, ii+1) = W(hh+1, ii+1) + R(jj+1, ii+1) .* exp(-2*pi*i*hh*jj/L);
80    end
81  end
82end
83
84