function resol(Nx,Nh,Nb,c);
% function resol(Nx,Nh,Nb,c);
% This function illustrates the difference between localization 
% and resolution on several linear FM chirp crossings.
%
% This function computes the reassigned spectrograms of three 
% signals made of two linear FM components with decreasing chirp 
% rates, and compares them with the ideal instantaneous
% frequency laws.
%
% input parameters:
%     NX  number of points
%     NB  frequency axis, freq. bin number
%     NH  (Gaussian) window length
%     C   color selection flag (C=0: bw, C=1, color)
%
% example:
%  Nx = 128;        % signal
%  Nh =  31;        % window
%  Nb = 128; c = 0; % TFD and display
%  resol(Nx,Nh,Nb,c);

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

% (time-frequency) window
h = window(Nh,'Gauss'); dh = (Nh-1)/2;

% corresponding time-frequency (ellipsoid) window base and height
a = dh; b = .5/dh; 
M = 500; phi = linspace(0,2*pi,M);

t = 1:Nx;
fact = 100;

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

% reassigned spectrogram, 1st config.
fmin = .15; fmax = .35;
[x1,if1] = fmlin(Nx,fmin,fmax);
[x2,if2] = fmlin(Nx,fmax,fmin);
x = x1 + x2;

[sp,rsp] = tfrrsp(x,t,2*Nb,h,1); clear sp;
rsp = rsp/max(max(rsp));
rsp = max(rsp,max(max(rsp))/fact);

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

subplot(2,3,1)
plot(1:Nx,if1,...
     1:Nx,if2,...
     Nx/2 + a*cos(phi),.25 + b*sin(phi),'r');
axis('square');
axis([1 Nx 0 0.5]);
set(gca,'XTick',[],'YTick',[])
xlabel(['window/signal = ',int2str(Nh),'/',int2str(Nx),' points']);
ylabel('model');

% plot reassigned spectrogram
subplot(2,3,4)
imagesc(flipud(log(rsp(1:Nb,:))))
axis('square');
set(gca,'XTick',[],'YTick',[]);
ylabel('reassigned spectrogram');

% reassigned spectrogram, 2nd config.
fmin = .2; fmax = .3;
[x1,if1] = fmlin(Nx,fmin,fmax);
[x2,if2] = fmlin(Nx,fmax,fmin);
x = x1 + x2;
[sp,rsp] = tfrrsp(x,t,2*Nb,h,1); clear sp;
rsp = rsp/max(max(rsp));
rsp = max(rsp,max(max(rsp))/fact);

subplot(2,3,2)
plot(1:Nx,if1,...
     1:Nx,if2,...
     Nx/2 + a*cos(phi),.25 + b*sin(phi),'r');

axis('square')
axis([1 Nx 0 0.5])
set(gca,'XTick',[],'YTick',[])

% plot reassigned spectrogram
subplot(2,3,5)
imagesc(flipud(log(rsp(1:Nb,:))))
axis('square')
set(gca,'XTick',[],'YTick',[]);

% reassigned spectrogram, 3rd config.
fmin = .22; fmax = .28;
[x1,if1] = fmlin(Nx,fmin,fmax);
[x2,if2] = fmlin(Nx,fmax,fmin);
x = x1 + x2;
[sp,rsp] = tfrrsp(x,t,2*Nb,h,1); clear sp;
rsp = rsp/max(max(rsp));
rsp = max(rsp,max(max(rsp))/fact);

subplot(2,3,3)
plot(1:Nx,if1,...
     1:Nx,if2,...
     Nx/2 + a*cos(phi),.25 + b*sin(phi),'r');
axis('square')
axis([1 Nx 0 0.5]);
set(gca,'XTick',[],'YTick',[]);

% plot reassigned spectrogram
subplot(2,3,6);
imagesc(flipud(log(rsp(1:Nb,:))));
axis('square');
set(gca,'XTick',[],'YTick',[]);
