% function frqtrack(x, gm, tao, prs, t, nf)
if 1
    global ubnd;
%    fs = 100;
    fs = 20;
    tt = (-2:(1/fs):2);
%   tt = (-1:(1/fs):1);
    ubnd = fs/2;
%     fs = 50;
   %  tt = 0:(1/fs):0.5;
    % tt = -2:(1/fs):2;
    % perc = 0.85;
    perc = 1;
    tr = randperm(length(tt));
    tc = sort(tr(1:floor(length(tt)*perc))); % truncated index of tt
     xns = 0; %1e-3; % 1e-3;
    t = tt(tc) + xns*rand(size(tc));
    w = 4;
%    aplmd = cos(2*pi*t*1);% + 2.5;
    % aplmd = exp(-t*0.8);
     xfg = 1;
     if xfg == 1
         amp = linspace(1,3,length(t));
         nslevel = 1e-1;
        % x = chirp(t,2,2,18,'q')'.*amp' ... %+  chirp(t,10,2,14)'...
         x = chirp(t,2,2,8,'q')'.*amp' ... %+  chirp(t,10,2,14)'...
        + nslevel*randn(size(t))';
        % x = chirp(t,15,2,35,'q')' ...
        %    + 1e-3*randn(size(t))';
    elseif xfg == 2
         x = chirp(t,15,1,35,'q')'  +  chirp(t,25,1,45)'...
          + 1e-3*randn(size(t))';
    elseif xfg == 3
        t = t+ 2;
         x = chirp(t,5,4,15)' + chirp(t,15,4,5)'...
             + 1e-3*randn(size(t))';        
         a=psd(x(1:10),[],fs);
    elseif xfg == 4
        x = (sin(2*pi*t*w).*aplmd)' ... % + (sin(2*pi*t*30).*aplmd)'...
        + 1e-4*randn(size(t))';
    elseif xfg == 5
         t = t+ 2;
         x = chirp(t,8,4,15)' + chirp(t,5,4,12)'...
             + 1e-3*randn(size(t))';        
    elseif xfg == 6
         t = t+ 2;
         x = cos(1./(3.01-t) + 8*pi*t)' + cos(1./(4.01-t)+8*pi*t)' ... % chirp(t,5,4,12)'...
             + 0*randn(size(t))';        
    elseif xfg == 7
        x = (sin(2*pi*t*3).*aplmd)' + (sin(2*pi*t*8).*aplmd)'...
        + 1e-4*randn(size(t))';
    end
    figure,specgram(x,30,fs,16); % Display the spectrogram
else
    n = 200
    s(:,1)= ones(2,1); x =[];
    s(1)=5; s(2)=10;
    for i = 2:n
        s(:,i) = s(:,i-1) + 1*randn(2,1); 
     %    x(i) = s(1,i)* log(s(2,i)^2) + randn(1);
         x(i) = s(1,i)* obs_h(s(2,i)) + randn(1);
     end
     x = colize(x);
     figure,subplot(313), plot(x);
     subplot(311), plot(s(1,:));
     subplot(312), plot((s(2,:)));
end

global bsz nf na np ns;
nf = 1; 
na = nf; % na = nf * 2 + 1;
% ns = nf * 3; % ns = nf * 3 + 1;
np = na + nf;
ns = np + nf;
bsz = 1;
[n,dx]=size(x);

nf3= 3*nf;
inda = [];  indp = []; indf = [];
inda = vec(1:na); 
indp = vec((na+1):(na+nf)); 
indf = vec(na+nf+1:na+2*nf);

tao = (4*(t(2)-t(1)))^2*eye(nf3) ; % state noise variance
% tao(inda,inda) = 0;
tao(inda,inda) = 1e-3*eye(nf);
tao(indp,indp) = 0;
% grnd = 4*t.^2 +2;
grnd = 1.5*t.^2 +2;
% tao(indf,indf) = mean(abs(grnd(2:end)-grnd(1:end-1)))^2*eye(nf);

% gm = 1e-1* eye(bsz); % observation noise variance
gm = nslevel^2 * eye(1); % observation noise variance

ms0 = zeros(nf3,1); %
ms0(indf) = fs/4;
ms0(indp) = 0;
ms0(inda) = 5;
vs0 = 1e2 * eye(nf3);
vs0(1,1)= 1e1;

prs= struct('m0',ms0,'v0',vs0);
mthd = 0;  nodes = 9;
%mthd = 1; nodes = 7*1e2;
x = vec(x); t = rowize(t);
wshift = 1; % actually overlap length = bsz-wshift-1
nfig = 1;
if 0
    [mst,tg] = frqtrackmc(x, gm, tao, prs, t, nf, nodes,1,mthd, wshift);
else
    [mst,tg] = frqtrackmc3(x, gm, tao, prs, t, nf, nodes, nfig, mthd, wshift);
end

nhd = 1;
plotmst(mst,tg, 5, nhd);
hold on,
subplot(311),plot(t,amp,'--')
hold on,
% plot(t,5*t.^2 +15,'--')
subplot(313),
plot(t,grnd,'--')


if 0

% QR for multinomial distri
mn = [2 2]';
nodes =9;
[ndps, w] = quadrl(nodes); % ndps is row vec.
ndps2=[ndps; zeros(size(ndps))];
ndps2=[ndps2 flipud(ndps2)];
ndps2 = ndps2 + repmat(mn,1,size(ndps2,2));
w2 = [w; w];
if 0
    lik = ( exp(-(ndps2(1,:))).^x(1) .* ...
     exp(-(ndps2(2,:))).^x(2) ) ./ ...
     (exp(-(ndps2(1,:))) + ...
    exp(-(ndps2(2,:))) ).^(sum(x)) * w2;
else
    allp =logsumexp( -ndps2, 1);
    mulp = exp(-ndps2 - repmat(allp,size(ndps2,1),1)) ;
    lik = (mulp(1,:).^x(1) .* mulp(2,:).^x(2) ) * w2
end

nsp = 1*1e5;
sp = normrnd_mt(nsp, mn, eye(length(mn)));
% sp = normrnd(2,1,2,nsp);
    allp =logsumexp( -sp, 1);
    mulp = exp(-sp - repmat(allp,size(sp,1),1)) ;
    lik = (mulp(1,:).^x(1) * mulp(2,:)'.^x(2) ) / nsp
    
end