function varwindow(Nx,Nh1,Nh2,Nh3,Nb,c);
% function varwindow(Nx,Nh1,Nh2,Nh3,Nb,c)
% compares the spectrogram and its reassigned version
% for three window lengths.
%
% input parameters:
%     NX  number of points
%     NH1 (Gaussian) window 1 length
%     NH2 (Gaussian) window 2 length
%     NH3 (Gaussian) window 3 length
%     NB  frequency axis, freq. bin number
%     C   color selection flag (C=0: bw, C=1, color)
%
% example:
%  Nx  = 128 ;                      % signal
%  Nh1 =  21 ; Nh2 = 63; Nh3 = 127; % windows
%  Nb  = 128 ;   c =  0;            % TFD and display
%  varwindow(Nx,Nh1,Nh2,Nh3,Nb,c);

% P. Flandrin, Feb. 12, 2001
% [Time-frequency reassignment: from principles to algorithms]
% Figure 4

% signal generation
fminl = .05; fmaxl = .45;
[xl,fxl] = fmlin(Nx,fminl,fmaxl);
x = xl.*amgauss(Nx,Nx/2,Nx/3);

% model plot
subplot(334); 
tt = fix(Nx/3):fix(2*Nx/3);
plot(tt,fxl(tt))
axis([1 Nx 0 0.5])
axis('square')
set(gca,'XTick',[],'YTick',[])
xlabel([int2str(Nx),' points']); ylabel('signal model')


% set the colormap
if c == 1
 colormap(jet)
else
 colormap(flipud(gray))
end

% define time vector 
% [time instants at which the spectrogram will be computed]
t = 1:Nx;

% compute the Wigner-Ville distribution of signal x
% at time intants t, with Nb (positive) frequency points.
% The last parameter is a progression indicator flag 
% [1=shown, 0=not shown]
wv = tfrwv(x,t,Nb,1);

threshold = 10^(-2);
wv = wv/max(max(wv)); 
wv=max(wv,threshold);

figure(gcf); % put the figure in the forefront

subplot(336);
imagesc(flipud(wv));
axis('square')
set(gca,'XTick',[],'YTick',[])
ylabel('Wigner-Ville')

% spectro (+ reass. spectro) with window Nh1
% with Nb (positive) frequency points and 2*Nb
% frequency bins (negative and positive)
h = window(Nh1,'Gauss');
[sp,rsp] = tfrrsp(x,t,2*Nb,h,1);

subplot(331); 
imagesc(flipud(sp(1:Nb,:)))
axis('square')
set(gca,'XTick',[],'YTick',[])
xlabel(['window = ',int2str(Nh1)]); ylabel('spectro')

subplot(337); 
imagesc(flipud(rsp(1:Nb,:)))
axis('square')
set(gca,'XTick',[],'YTick',[])
xlabel(['window = ',int2str(Nh1)]); ylabel('reass. spectro')

% spectro (+ reass.) with window Nh2
h = window(Nh2,'Gauss');
[sp,rsp] = tfrrsp(x,t,2*Nb,h,1);

subplot(332); 
imagesc(flipud(sp(1:Nb,:)))
axis('square')
set(gca,'XTick',[],'YTick',[])
xlabel(int2str(Nh2))

subplot(338); 
imagesc(flipud(rsp(1:Nb,:)))
axis('square')
set(gca,'XTick',[],'YTick',[])
xlabel(int2str(Nh2))

% spectro (+ reass.) with window Nh3
h = window(Nh3,'Gauss');
[sp,rsp] = tfrrsp(x,t,2*Nb,h,1);

subplot(333); 
imagesc(flipud(sp(1:Nb,:)))
axis('square')
set(gca,'XTick',[],'YTick',[])
xlabel([int2str(Nh3),' points'])

subplot(339); 
imagesc(flipud(rsp(1:Nb,:)))
axis('square')
set(gca,'XTick',[],'YTick',[])
xlabel([int2str(Nh3),' points'])
