%% testGabMulExp_new recreates Figure 1 from 'A survey on the 
% unconditional convergence and the invertibility of frame 
% multipliers with implementation
%
% This file reads in a predetermined symbol from a PNG file and 
% creates and applies the associated multiplier.
%
% References:
% [1]  Diana T. Stoeva and Peter Balazs, A survey on the 
%       unconditional convergence and the invertibility of 
%      frame multipliers with implementation, submitted
%
% [2]  P. Soendergaard, B. Torresani, P. Balazs, "The Linear Time
%        Frequency Analysis Toolbox", IJWMIP, Vol. 10 (4), (2012) 


ltfatstart
% This scirpt uses the LTFAT toolbox, see [2] 

% The symbol was created using the MulacLab GUI for multiplier. 
% 
% If you want to do the same for another signal, in MulacLab save 
% the original frame, and laod it in this file:
% Orig = load('Orig.mat');
%
% wintype = Orig.frame.def.winType;
% winlen = Orig.frame.def.winLen; 
% a = Orig.frame.def.hop; 
% M = Orig.frame.def.nbFreq;
% sig = Orig.savedSig.ori;
% SR = Orig.savedSig.sampFreq;
%
% Here we use
wintype = 'hann'; % window type, here: Hanning window
winlen = 1024;    % window length in samples
a = 256;          % hop size, i.e time shift parameter
M = 2048;         % Number of frequency bins, i.e. length of FFT
[sig,SR] = audioread('originalsignal.wav'); % load sound file

%  If you have used MulacLab 
% sig = Orig.savedSig.ori
% Export the symbol to a PNG file, open it in a graphics program, 
% seperate the layers, and save one of them (this should be black 
% and white then.
% 
% Also choose the following according to your  saved ymbol
inv_mult = 1; % invert symbol?
set_nonzero = 1; % set the minimum value above zero?

% length of signal and number of channels:
[Ls,Channels] = size(sig);
% restrict to one channel:
sig = sig(:,1);
sound(sig,SR);

% Choose parmeters for a nice plot
maxfreq = 6000; % maximal displayed frequency
dynrange = 50; % dynamical range

%% choose a nice colormap
xcm = [         0         0    0.2000
         0         0    0.2076
         0         0    0.2152
         0         0    0.2229
         0         0    0.2305
         0         0    0.2381
         0         0    0.2457
         0         0    0.2533
         0         0    0.3778
         0         0    0.5022
         0         0    0.6267
         0         0    0.7511
         0         0    0.8756
         0         0    1.0000
         0    0.0857    0.9143
         0    0.1714    0.8286
         0    0.2571    0.7429
         0    0.3429    0.6571
         0    0.4286    0.5714
         0    0.5143    0.4857
         0    0.6000    0.4000
    0.0909    0.6364    0.3636
    0.1818    0.6727    0.3273
    0.2727    0.7091    0.2909
    0.3636    0.7455    0.2545
    0.4545    0.7818    0.2182
    0.5455    0.8182    0.1818
    0.6364    0.8545    0.1455
    0.7273    0.8909    0.1091
    0.8182    0.9273    0.0727
    0.9091    0.9636    0.0364
    1.0000    1.0000         0
    1.0000    0.8750         0
    1.0000    0.7500         0
    1.0000    0.6250         0
    1.0000    0.5000         0
    1.0000    0.3750         0
    1.0000    0.2500         0
    1.0000    0.1250         0
    1.0000         0         0
    1.0000    0.0196    0.0196
    1.0000    0.0392    0.0392
    1.0000    0.0588    0.0588
    1.0000    0.0784    0.0784
    1.0000    0.1936    0.1936
    1.0000    0.3088    0.3088
    1.0000    0.4240    0.4240
    1.0000    0.5392    0.5392
    1.0000    0.6544    0.6544
    1.0000    0.7696    0.7696
    1.0000    0.8848    0.8848
    1.0000    1.0000    1.0000
    1.0000    1.0000    1.0000
    1.0000    1.0000    1.0000
    1.0000    1.0000    1.0000
    1.0000    1.0000    1.0000
    1.0000    1.0000    1.0000
    1.0000    1.0000    1.0000
    1.0000    1.0000    1.0000
    1.0000    1.0000    1.0000
    1.0000    1.0000    1.0000
    1.0000    1.0000    1.0000
    1.0000    1.0000    1.0000
    1.0000    1.0000    1.0000];


% Do a frame analysis, here using the DGT real
[FA,FS] = framepair('dgtreal',{wintype,winlen},'dual',a,M);
Coeff_Orig = frana(FA,sig); % frame analysis,
% Depending on the frame, sometimes the length has to be adapted:
L=framelength(FA,Ls);

%% plot the original multipliers
figure(1);
% Original    
subplot(2,2,1);
plotframe(FA,Coeff_Orig,SR,'dynrange',dynrange);
set(gca,'YLim',[0 maxfreq]); % restrict the frequency range
CCC = get(gca,'CLim'); % get range for values
title(sprintf('TF representation of original signal'));
colormap(xcm);

 %% multiplier creation
% choose a minimal number to have a semi-normalized symbol:
 mmss = 0.01;
    
% read symbol from file
mask_orig=double(imread('Symbol6_BW.png'));

mask = mask_orig(:,:,1); % restircting to one color channel
mask = flipud(mask); % Flip upside down 
%                                                (orientation was changed in the PNG)
MMMM = max(mask(:)); % Maximum value of symbol
mask = mask/MMMM; % normalize to one

% Invert symbol if wanted:
 if inv_mult == 1
    mask = (1-mask);
 end

% Set small values to minimal values (if wanted)
if set_nonzero == 1
       iii = find(mask<mmss);
       mask(iii) = mmss;
end

% transfer from frame-specific shape to column
mask = framenative2coef(FA,mask); 

%% plot the symbol
subplot(2,2,2);
plotframe(FA,mask,SR,'dynrange',dynrange);
set(gca,'YLim',[0 maxfreq]);
title(sprintf('Symbol'));
colormap(xcm);

%% plot the multiplication in the TF domain
 % do the multiplication in the TF domain: 
MulAppl = Coeff_Orig.*(mask);
% plot
subplot(2,2,3);
plotframe(FA,MulAppl,SR,'dynrange',dynrange);
set(gca,'YLim',[0 maxfreq])
set(gca,'CLim',CCC);
title('Multiplication in TF domain');
 
%% plot  the analysis of the multiplier
% this is equivalent to a projection of the multplication above to the 
% range of the analysis operator
subplot(2,2,4);
% apply the frame multiplier:
filtered_signal = framemul(sig,FA,FS,mask);
% this could be done equivalenty by applying the synthesis on the 
% above multiplication
%
% Write to file:
audiowrite('modifiedsignal.wav',filtered_signal,SR);
% Play modified signal:
sound(filtered_signal,SR);

% Anaöyse the result: 
Coeff_Filt = frana(FA,filtered_signal);
plotframe(FA,Coeff_Filt,SR,'dynrange',dynrange);
set(gca,'YLim',[0 maxfreq])
set(gca,'CLim',CCC);
title('TF representation of modified signal'); 
% Save Figure as PNG-file
  saveas(gcf,'testGabMulExp_new.png','png');    
