function quasilin(Nx,Nh,Nb,Nhp,c);
% function quasilin(Nx,Nh,Nb,Nhp,c);
%
% computes and displays the reassigned spectrogram of a 
% sinusoid with sinusoidal frequency modulation.
%
% input parameters:
%  NX  number of points
%  NH  (Gaussian) window length
%  NB  frequency axis, freq. bin number
%  NHP half-period of the sinusoidal FM  
%  C   color selection flag (C=0: bw, C=1, color)
%
% example:
%  Nx = 256; Nhp = 21; % signal
%  Nh =  19;           % window
%  Nb = 128;    c = 0; % TFD and display
%  quasilin(Nx,Nh,Nb,Nhp,c);

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

% signal generation (sinus FM)
fmins = 0.1; fmaxs = 0.4;
SinPeriod = 2*Nx/Nhp;
[xs,fxs] = fmsin(Nx,fmins,fmaxs,SinPeriod);

% TF model plot
tt = 1:Nx;

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

subplot(2,1,1);
plot(tt,fxs(tt));
axis([1 Nx 0 0.5]);
set(gca,'XTick',[],'YTick',[]);
xlabel('time'); ylabel('frequency');
title('signal model');

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

% TF computation
dt = 1; td = 1; tf = Nx; 
t = td:dt:tf; 
h = window(Nh,'Kaiser'); % define the window

% define the corresponding time-frequency (ellipsoid) 
% window base and height
dh = (Nh-1)/2; a = dh; b = 0.5/dh; 
M = 500; phi = linspace(0, 2*pi, M);

% compute the reassigned spectogram of signal x
% at time intants t with 2*Nb frequency points 
% [i.e., 2*Nb samples for the FFT time-base]
% using window h. The last parameter is a progression 
% indicator flag [1=shown, 0=not shown]

[s,rs] = tfrrsp(xs,t,2*Nb,h,1);

% change the value of rs at one point so that
% the colormap does not vary with time.
rs(1,1) = max(max(rs));

% shows a progressive plot of the reassigned spectrogram
for k = td:dt:tf;   %ceil(15*Nx/Nhp)
 subplot(211)
 plot(tt,fxs(tt),':',...
 td:k,fxs(tt(td:k)),...
 tt(k) + a*cos(phi),fxs(tt(k)) + b*sin(phi),'r',...
 tt(k),fxs(tt(k)),'r+');

 axis([1 Nx 0 0.5])
 set(gca,'XTick',[],'YTick',[])
 xlabel('time'); ylabel('frequency')
 title('model')

 subplot(2,1,2)
 imagesc(flipud(rs(1:Nb,1:1+(k-td)/dt)));
 axis([1 max(size(t)) 1 Nb])
 set(gca,'XTick',[],'YTick',[])
 xlabel('time'); ylabel('frequency')
 title('reassigned spectrogram')

 pause(.3)

end
