% sparsetfr_test.m 
%
% Written by: Patrick Flandrin and Pierre Borgnat, CNRS-ENS de Lyon
% Email: {flandrin,pborgnat}@ens-lyon.fr
% Created: April 2015

%% 1. creating the signal

% up-to-6-tone synthetic signal 
n = 256 ; % signal length
[s1,fs1] = fmlin(n,0.07,0.19) ; a1 = 1*amgauss(n,n/2,n/2) ; %  ; %fmconst(n,0.07) ,0.19
[s2,fs2] = fmconst(n,0.14) ; a2 = 0*amgauss(n,n/2,n/2)/1 ; % fmlin(n,0.05,0.4) 
[s3,fs3] = fmconst(n,0.21) ; a3 = 1*amgauss(n,n/2,n/2)/1 ; % fmlin(n,0.3,0.1)
[s4,fs4] = fmconst(n,0.28) ; a4 = 0*amgauss(n,n/2,n/2)/1 ;
[s5,fs5] = fmsin(n,0.3,0.4,n/2) ; a5 = 1*amgauss(n,n/2,n/2)/1 ; %  ; %fmconst(n,0.35)
[s6,fs6] = fmconst(n,0.42) ; a6 = 0*amgauss(n,n/2,n/2)/1 ;
s = a1.*s1 + a2.*s2 + a3.*s3 + a4.*s4 + a5.*s5 + a6.*s6 ; 

% noise (if any)
s = s + randn(size(s))/3 ; %3
% 
% % normalization to unit energy
% s = s/sqrt(sum(abs(s).^2)) ;

% 2. Computing TF references and plotting for check
N = n ;
    
% time-frequency model
fs1d = min(N,1+round(2*(N)*fs1)) ; 
fs2d = min(N,1+round(2*(N)*fs2)) ; 
fs3d = min(N,1+round(2*(N)*fs3)) ; 
fs4d = min(N,1+round(2*(N)*fs4)) ; 
fs5d = min(N,1+round(2*(N)*fs5)) ; 
fs6d = min(N,1+round(2*(N)*fs6)) ; 
TF1 = zeros(N) ;
TF2 = zeros(N) ;
TF3 = zeros(N) ;
TF4 = zeros(N) ;
TF5 = zeros(N) ;
TF6 = zeros(N) ;
for kk = 1:N
    TF1(fs1d(kk),kk) = a1(kk).^2 ;
    TF2(fs2d(kk),kk) = a2(kk).^2 ;
    TF3(fs3d(kk),kk) = a3(kk).^2 ;
    TF4(fs4d(kk),kk) = a4(kk).^2 ;
    TF5(fs5d(kk),kk) = a5(kk).^2 ;
    TF6(fs6d(kk),kk) = a6(kk).^2 ;
end
TF = TF1 + TF2 + TF3 + TF4 + TF5 + TF6 ;

% (ressigned) spectrogram
[Xsp,Xrsp]=tfrrsp(s,1:N,2*N,tftb_window(49,'Gauss'),0);

% Wigner-Ville distribution
X = tfrwv(s) ;  
x = X(:) ;

% ambiguity function
AF = fliplr(fftshift(fft2(X)')) ; 
AF = AF/max(max(abs(AF))) ;

% check plot
figure(1)
dyn = 20 ; % dynamic range (in dB)
subplot(2,2,1); imageTF(TF,dyn,1,0); title('model')
subplot(2,2,2); imageTF(Xsp(1:N,:),dyn,1,0); title('spectrogram')
subplot(2,2,3); imageTF(X,dyn,1,0); title('Wigner-Ville distribution')
subplot(2,2,4); imageTF(Xrsp(1:N,:),dyn,1,0); title('reassigned spectrogram')

%% 3. computing TF sparse solutions

larg = 12 ;
epsilon = 5e-1 ;

tic

[ TF_l1, TF_l2, epsilon ] = sparsetfr( s, larg, epsilon )
    
toc

% 4. plotting results

figure(2)

% square mask in the central region of the ambiguity function
M = zeros(N) ;
M(1+n/2-larg/2:1+n/2+larg/2,1+n/2-larg/2:1+n/2+larg/2) = ones(larg+1,larg+1) ;

% display
subplot(3,2,1); imageTF(TF,dyn,1,0); title('model')
subplot(3,2,5); imageTF(Xrsp(1:N,:),dyn,1,0); title('reassigned spectrogram')
subplot(3,2,3); imageTF(X,dyn,1,0); title('Wigner-Ville distribution')
subplot(3,2,2); imageAF(abs(AF),dyn,1,0) ; 
hold on
plot([n/2-larg/2 n/2+larg/2],[n/2-larg/2 n/2-larg/2],'k','LineWidth',3)
plot([n/2-larg/2 n/2+larg/2],[n/2+larg/2 n/2+larg/2],'k','LineWidth',3)
plot([n/2-larg/2 n/2-larg/2],[n/2-larg/2 n/2+larg/2],'k','LineWidth',3)
plot([n/2+larg/2 n/2+larg/2],[n/2-larg/2 n/2+larg/2],'k','LineWidth',3)
plot([N/2 N/2],[1 N],'k')
plot([1 N],[N/2 N/2],'k')
hold off
title(['ambiguity function and mask (d = ',int2str(larg+1),'x',int2str(larg+1),')'])
subplot(3,2,4); imageTF(TF_l2,dyn,1,0); title('l_2 backprojection')
subplot(3,2,6); imageTF(TF_l1,dyn,1,0); title(['l_1 reconstruction (eps = ',num2str(epsilon),')']);

