% MTFR_comp.m 
%
% J. Xiao and P. Flandrin, 2006
%
% Computes and displays time-frequency estimates based on 
% Hermite multitaper reassigned spectrograms in the case 
% of synthetic signals embedded in noise.
%
% Allows to reproduce Figures 4, 6 and 7 of the manuscript
% "Multitaper Time-Frequency Reassignment for Nonstationary
% Spectrum Estimation and Chirp Enhancement" by J. Xiao and
% P. Flandrin, submitted for publication (2006), cf.
% http://perso.ens-lyon.fr/patrick.flandrin/multitfr.html
%
% calls  - tfrrsp_hm.m
%        - imageTF.m
%        - tfrrsp_h.m
%        - hermf.m
%        - mean_hm.m
%        - tvnoise2.m
%
% needs  - Time-frequency ToolBox http://tftb.nongnu.org/
%        - WaveLab http://www-stat.stanford.edu/~wavelab/

% 3 possible types of signal

    % Stype = 1 : 'pure tone'
    % Stype = 2 : 'sinus FM'
    % Stype = 3 : 'sinus FM + linear chirp'
    
% 2 possible types of noise

    % Ntype = 1 : 'white Gaussian noise'
    % Ntype = 2 : 'band-pass time-varying white Gaussian noise'
    % Ntype = 3 : 'band-pass time-limited white Gaussian noise'
    
% parameters

    % Nx   : signal length
    % t    : time vector
    % Nfft : number of frequency bins
    % Mh   : number of tapers
    % Nh   : tapers length
    % tm   : tapers support
    % SNR  : signal-to-noise ratio (dB)
    % Td   : threshold on "difference" estimates
    % fact : factor of proportion for threshold on multitaper estimate
    % dyn  : dynamic range of time-frequency representations display (dB)
  
% nominal choices in the manuscript

    % Figures 4, 6 and 7 : 
    %   Nx = 512 ; t = 1:1:Nx ; 
    %   Nfft = Nx ; 
    %   Nh = 155 ; tm = 6 ; 
    %   Td = 1.5 ; fact = 7 ; 
    %   dyn = 30 ; 

    % Figure 4 :
    %   Stype = 1 ; Ntype = 1 ;
    %   SNR = 10 ; Mh = 4 ; 
    
    % Figure 6 : 
    %   Stype = 2 ; Ntype = 2 ;
    %   SNR = 5 ; Mh = 8 ; 
    
    % Figure 7 : 
    %   Stype = 3 ; Ntype = 3 ; 
    %   SNR = 10 ; Mh = 8 ; 

% actual user's choice -----------------%

    Nx = 512 ; t = 1:1:Nx ;             %
    Nfft = Nx ;                         %
    Nh = 155 ; tm = 6 ;                 % 
    Td = 1.5 ; fact = 7 ;               %
    dyn = 30 ;                          %

    Stype = 3 ; Ntype = 3 ;             %  
    SNR = 10 ; Mh = 8 ;                 %

%---------------------------------------%
    
% signal, noise and signal+noise generation

if Stype == 1

    signal = real(fmconst(Nx,1/4)) ;

elseif Stype == 2

    [sig,iflaw] = fmsin(Nx,0.15,0.45,2*Nx,Nx/4,0.3,-1) ;
    signal = real(sig).*amgauss(Nx,Nx/2,10*Nx) ;

elseif Stype == 3

    [sig,iflaw] = fmsin(Nx,0.15,0.45,2*Nx,Nx/4,0.3,-1) ;
    d = 0.1 ;
    signal = real(sig).*amgauss(Nx,Nx/2,10*Nx) + real(fmlin(Nx,0.25+d,0.25-d)) ;

end

if Ntype == 1

    tdep  = 1 ; 
    tarr = Nx ; 
    fcdep = 1/4 ;
    fcarr = 1/4 ;
    B = 1/2 ;
    noise = randn(Nx,1) ; 

elseif Ntype == 2

    tdep  = 1 ; 
    tarr = Nx ; 
    fcdep = 1/8 ;
    fcarr = 3/8 ;
    B = 1/8 ;
    noise = tvnoise2(Nx,tdep,tarr,fcdep,fcarr,B) ;

elseif Ntype == 3

    tdep  = 3*Nx/4 ; 
    tarr = 7*Nx/8 ; 
    fcdep = 3/8 ;
    fcarr = 3/8 ;
    B = 1/8 ;
    noise = tvnoise2(Nx,tdep,tarr,fcdep,fcarr,B) ;

end

x = sigmerge(signal,noise,SNR) ;

% initialization

S = [] ; RS = [] ; Smean = [] ; RSmean = [] ;
RSN = [] ; RSD = [] ; RSDmean = [] ; 
Mask = [] ; Pre = [] ; RSM = [] ; RSP = [] ; RSF = [] ; 

% computation of spectrograms (S), reassigned spectrograms (RS) and their
% means (Smean and RSmean)

[S,RS,hat,tt] = tfrrsp_hm(x,t,Nfft,Nh,Mh,tm) ;
[Smean, RSmean] = mean_hm(S,RS,hat,1) ;

MMM = max(max(max(RS))) ;

for k = 1:Mh   % normalization (RSN) of RS
    
    RSN(:,:,k) = max(RS(:,:,k),10^(-dyn/10)*MMM) ;
    
end 

% computation of "differences", masks and enhanced estimates

for k = 1:Mh-1   % "difference" (RSD), mean (RSDmean) and normalization (RSDmean) of RSN
        
    RSD(:,:,k) = exp(log(RSN(:,:,k+1)) - log(RSN(:,:,k))) ; 
    
end

RSDmean = mean(RSD,3) ; 
    
for m = 1:Nfft    % threshold Td on RSDmean
    for n = 1:Nx
        if RSDmean(m,n) >= Td
            Mask(m,n) = 0 ;
        else 
            Mask(m,n) = 1 ; 
        end
    end
end

ts = mean(mean(RSmean));

for m = 1:Nfft    % threshold Ts on RSmean
    for n = 1:Nx
%         if RSmean(m,n) >= ts*fact
        if abs(RSmean(m,n)-1)+1 >= ts*fact
            Pre(m,n) = 1 ;
        else 
            Pre(m,n) = 0 ; 
        end
    end
end

RSM = RSmean.*Mask ;   % "difference" average
RSP = RSmean.*Pre ;    % "pre-denoised" average
RSF = sqrt(RSP.*RSM) ; % combined estimate

% visualization

c = 0 ; % black and white colormap ('jet' : c = 1)
b = 0 ; % no colorbar (yes : b = 1)

figure(1)

subplot(231) % spectrogram
imageTF(S(1:Nfft/2,:,1),dyn,c,b) ;
title(['S1'])
hold on
plot([tdep tarr],[Nfft*(1/2-fcdep+B/2) Nfft*(1/2-fcarr+B/2)],'k')
plot([tdep tarr],[Nfft*(1/2-fcdep-B/2) Nfft*(1/2-fcarr-B/2)],'k')
plot([tdep tdep],[Nfft*(1/2-fcdep+B/2) Nfft*(1/2-fcdep-B/2)],'k')
plot([tarr tarr],[Nfft*(1/2-fcarr+B/2) Nfft*(1/2-fcarr-B/2)],'k')
hold off
xlabel('') ; ylabel('') ;

subplot(232) % reassigned spectrogram
imageTF(RS(1:Nfft/2,:,1),dyn,c,b) ;
title(['RS1'])
hold on
plot([tdep tarr],[Nfft*(1/2-fcdep+B/2) Nfft*(1/2-fcarr+B/2)],'k')
plot([tdep tarr],[Nfft*(1/2-fcdep-B/2) Nfft*(1/2-fcarr-B/2)],'k')
plot([tdep tdep],[Nfft*(1/2-fcdep+B/2) Nfft*(1/2-fcdep-B/2)],'k')
plot([tarr tarr],[Nfft*(1/2-fcarr+B/2) Nfft*(1/2-fcarr-B/2)],'k')
hold off
xlabel('') ; ylabel('') ;

subplot(233) % multitaper reassigned spectrogram
imageTF(RSmean(1:Nfft/2,:),dyn,c,b) ;
title(['RS'])
hold on
plot([tdep tarr],[Nfft*(1/2-fcdep+B/2) Nfft*(1/2-fcarr+B/2)],'k')
plot([tdep tarr],[Nfft*(1/2-fcdep-B/2) Nfft*(1/2-fcarr-B/2)],'k')
plot([tdep tdep],[Nfft*(1/2-fcdep+B/2) Nfft*(1/2-fcdep-B/2)],'k')
plot([tarr tarr],[Nfft*(1/2-fcarr+B/2) Nfft*(1/2-fcarr-B/2)],'k')
hold off
xlabel('') ; ylabel('') ;

subplot(234) % "pre-denoised" average 
imageTF(RSP(1:Nfft/2,:),dyn,c,b) ;
title(['RSP'])
hold on
plot([tdep tarr],[Nfft*(1/2-fcdep+B/2) Nfft*(1/2-fcarr+B/2)],'k')
plot([tdep tarr],[Nfft*(1/2-fcdep-B/2) Nfft*(1/2-fcarr-B/2)],'k')
plot([tdep tdep],[Nfft*(1/2-fcdep+B/2) Nfft*(1/2-fcdep-B/2)],'k')
plot([tarr tarr],[Nfft*(1/2-fcarr+B/2) Nfft*(1/2-fcarr-B/2)],'k')
hold off
xlabel('') ; ylabel('') ;

subplot(235) % "difference" average
imageTF(RSM(1:Nfft/2,:),dyn,c,b) ;
title(['RSM'])
hold on
plot([tdep tarr],[Nfft*(1/2-fcdep+B/2) Nfft*(1/2-fcarr+B/2)],'k')
plot([tdep tarr],[Nfft*(1/2-fcdep-B/2) Nfft*(1/2-fcarr-B/2)],'k')
plot([tdep tdep],[Nfft*(1/2-fcdep+B/2) Nfft*(1/2-fcdep-B/2)],'k')
plot([tarr tarr],[Nfft*(1/2-fcarr+B/2) Nfft*(1/2-fcarr-B/2)],'k')
hold off
xlabel('') ; ylabel('') ;

subplot(236) % combined estimate
imageTF(RSF(1:Nfft/2,:),dyn,c,b) ;
title(['RSF'])
hold on
plot([tdep tarr],[Nfft*(1/2-fcdep+B/2) Nfft*(1/2-fcarr+B/2)],'k')
plot([tdep tarr],[Nfft*(1/2-fcdep-B/2) Nfft*(1/2-fcarr-B/2)],'k')
plot([tdep tdep],[Nfft*(1/2-fcdep+B/2) Nfft*(1/2-fcdep-B/2)],'k')
plot([tarr tarr],[Nfft*(1/2-fcarr+B/2) Nfft*(1/2-fcarr-B/2)],'k')
hold off
xlabel('') ; ylabel('') ;

