function f=dwv(s,M);

%                       f=dwv(s,M);
%
% Calcule la transformée de Wigner-Ville discrète d'un signal digital s
% donné en colonne. Le signal est prolongé par des zéros de part et d'autre
% de son intervalle de définition, après avoir été préalablement multiplié
% par une fenêtre de Hamming. L'axe des ordonnées est l'axe des fréquences,
% la plage fréquentielle utile étant [0, 2*pi], avec la convention qu'un
% réel omega entre pi et 2*pi correspond à la fréquence omega - 2*pi. Le
% calcul de transformée de Fourier est un calcul de M-dft avec M supérieur
% ou égal à la longueur N de s.  

N=length(s);
y=[zeros(N,1); s.*hamming(N); zeros(N,1)];
A=[];
T1=ones(N,1);
T2=ones(N,1);
for k=1:N,
  T1=y((N+2-k):(2*N+1-k));
  T2=conj(y((N+k):(2*N+k-1)));
  T=T1.*T2;
  A=[A T];
end
f=fft(A');