% reprInvMult.m: Gabor Multipliers applied on sound files 
function MultMat = reprInvMult(filename)
%
% This test file demonstrates frame multipliers for sound processing. 
% Furthermore it shows in how far multipliers can be inverted. We reproduce
% the figures from [2].
% 
% For the Gabor transform the painless approach is used. This file uses the
% LTFAT frame approach to analysis and synthesis, to allow easy adaptation
% to other frames (than Gabor) in the future.
%
% For the creation of the symbol/mask, we have chosen to use an external
% graphical program. In the first Step (set Step to 1), the Spectrogram is
% saved to a JPG file. From that, using for example any graphical program
% another JPG file has to be created, which represents the multiplier's
% symbol, i.e. the mask.
%
% In Step 2 this file is read and applied.
% 
% There exists a GUI for multipliers in LTFAT, named MULACLAB. 
% The approach with a 3rd party graphical program was chosen to give a
% more complete set of tools for the generation of the symbol/mask.
%
% other code:
% - using the linear Time Frequency Toolbox (LTFAT) v.1.2.0
% 
% References:
% [1] P. Soendergaard, B. Torresani, P. Balazs, "The Linear Time Frequency 
%   Analysis Toolbox",  International Journal of Wavelets, Multiresolution 
%   and Information Processing, Vol. 10 (4), 1250032 (2012)
% [2] D. T. Stoeva, P. Balazs, "Representation of the inverse of a 
%   multiplier as a multiplier", submitted (2013)
%
% Authors:
% Peter Balazs, XXL: contact: peter.balazs@oeaw.ac.at
% last modification 20.03.2013
%
% This work is licensed under the Creative Commons Attribution-NonCommercial-ShareAlike 3.0 Unported License.
% To view a copy of this license, visit http://creativecommons.org/licenses/by-nc-sa/3.0/

%% version

ver_RINV = '10th July 2013 Final version';
disp(sprintf('File version of %s.',ver_RINV));


%% internal options

inoutfiles = 1; % save new JPG and wav files ? 
playfiles = 0; % play the sound files?

inv_mult = 1; % invert symbol?

calc_condnum = 0; % calculate the condition number? (very time-consuming) 

Step = 2; % Step:
% Step 1: A picture of the TF-representation is stored, to be used in a 
%         graphical program (e.g. GIMP) to generate a mask  
%
% Step 2: The picture of the mask is loaded and processed! 

%% intialization of outside code
ltfatstart; % intialize LTFAT

%% loading the signal
if nargin < 1
    filename = 'jumpshortorig';
end
[sig,fs,nBits,OPTS] = wavread(sprintf('%s.wav',filename));

% structure for signal
signal.fs = fs;              %sampling rate
signal.nBits = nBits;        % bit rate
signal.nCh = OPTS.fmt.nChannels; % number of channels
signal.name = filename;      %filename

[Ls mmm] = size(sig);
if mmm ~= signal.nCh
    disp('Iternal error: number of channels');
    clear signal;
    return
end
signal.sig = sig;             % samples
signal.len = Ls;             % length of signal
signal.unit.Temp = 'Samples';
signal.unit.Amp  = 'linear'; % linear from -1 to 1

if playfiles == 1
    wavplay(signal.sig,signal.fs)
end


%% parameters for Gabor transform and plot
Nwin_ms = 80;                    % length of window in ms
Nwin = round(Nwin_ms*fs/1000);   % length of window
winname = 'hanning';             % window type
% g = firwin(winname,Nwin);       % analysis window
% gs = gabdual(g,a,M); % synthesis window, here: dual window

overlap = 8;                %overlap, redundancy in time
a = round(Nwin/overlap);         % time shift, hop size
zeropad = 1;           % use zeropadding? redundancy in frequency.
M = Nwin*(zeropad);            % number of channels, length of FFT (do it painless)

% plot parameters
maxfreq = 6000; % maximal displayed frequency
dynrange = 50; % dynamical range
xcm = nicecolormap; % choose the nice colormap

%% frame creation
% using the frame concept in LTFAT allows an easy extension later.
% e.g. to wavelets, NSGT or general frames
[Fa,Fs]=framepair('dgtreal',{winname,Nwin},'dual',a,M); % create the frames  

start_dual = 1;          % we started out with a dual pair!                               
% for other pair of windows use:
 % Fa=frame('dgtreal',g,a,M);
 % Fs=frame('dgtreal',gs,a,M);

% calculate the missing TF parameter 
L=framelength(Fa,Ls);
% b = L/M;
% N = L/a;



Coeff_Orig = frana(Fa,signal.sig); % frame analysis,
% for the discrete Gabor transform equivalent to dgtreal(signal.sig,g,a,M); 

%% plot the original signal                                    

figure(1);
if Step == 2     
    subplot(2,2,1)
end

plotframe(Fa,Coeff_Orig,signal.fs,'dynrange',dynrange);
set(gca,'YLim',[0 maxfreq]);
CCC = get(gca,'CLim');
title(sprintf('TF representation of original signal'));
colormap(xcm);

%% mask/symbol for the multiplier 
if Step == 1 
    %% Write to file and then stop, to be manipulated outside to create the mask.   
    % Standard program is called to open JPG file.
    saveas(gcf,sprintf('%s_testGabMult.jpg',signal.name),'jpg');   
    winopen(sprintf('%s_testGabMult.jpg',signal.name));
else    
    %% reads the symbol from outside 
    mmss_dB = -40; % minimal number in dB
    mmss = 10^(mmss_dB/20);
    % read symbol from file
    % save it to the right filename !!
    mask_orig=double(imread(sprintf('%s_mask_bw.png',filename)));
    
    mask = mask_orig;
    if inv_mult == 1
        mask = (1-mask);
    end
   
    % raise to the mimimum
    iii = find(mask<mmss);
    mask(iii) = mmss;
    
    mask = mask(:); % for the frame framework a vector instead of a matrix is needed.
    
    % plot the mask
    subplot(2,2,2);
    plotframe(Fa,mask,signal.fs);
    set(gca,'YLim',[0 maxfreq]);
    title('Gabor mask: symbol of multiplier');
                                                                             
    %% condition number of multiplier                               
    
    if calc_condnum == 1
        B=framemuleigs(Fa,Fs,mask,1,'iter','print','tol',1e-3)

        opts.disp=1;
        opts.isreal = 1;
        opts.maxit  = 100;
        opts.tol    = 1e-3;

        Ainv=eigs(@(x) iframemul(x,Fa,Fs,mask,'maxit',100,'tol',1e-3),L,1,'LM',opts)

        disp(sprintf('Condition number: %e',B*Ainv));
    end
        
    %% TF plot of the result of the multiplier 
    subplot('Position', [.35, 0.1, .35, .35]);
    filtered_signal = framemul(signal.sig,Fa,Fs,mask);
   
    Coeff_Filt = frana(Fa,filtered_signal);
    plotframe(Fa,Coeff_Filt,signal.fs,'dynrange',dynrange);
    set(gca,'YLim',[0 maxfreq])
    set(gca,'CLim',CCC);
    title('TF representation of modified signal');
    
    % Play result!
    if playfiles == 1
         xyz = 'y';
         while xyz == 'y'
           wavplay(filtered_signal,signal.fs);
            xyz = input('Replay filtered signal? y/n [y]: ', 's');
            if isempty(xyz)
                 xyz = 'y';
             end
         end
    end
    if inoutfiles == 1
        wavwrite(real(filtered_signal),signal.fs,signal.nBits,sprintf('filtersignal_%s',filename));
        saveas(gcf,sprintf('multipliers_%s.eps',filename),'eps');
        saveas(gcf,sprintf('multipliers_%s.png',filename),'png');
        saveas(gcf,sprintf('multipliers_%s.jpg',filename),'jpg');
        saveas(gcf,sprintf('multipliers_%s.fig',filename),'fig');
    end
 
    %% multiplier inversion 
    % for the naive approach use the canonical duals of Fa and Fs.    
    
    if start_dual == 1
        Fad = Fs;
        Fsd = Fa;
    else
        Fad = framedual(Fa);
        Fsd = framedual(Fs);
    end
    
    % apply naive inverse multiplier 
    recon_naive = framemul(filtered_signal,Fsd,Fad,1./mask);
    recon_naive = recon_naive(1:Ls); % cut to original length
    Coeff_naive = frana(Fa,recon_naive);
    
    figure(2);
    
    subplot(2,2,1);
    plotframe(Fa,Coeff_naive,signal.fs,'dynrange',dynrange);
    set(gca,'YLim',[0 maxfreq])
    set(gca,'CLim',CCC);
    title('TF representation of naive inversion');
    colormap(xcm); 
    
    Err_naive = recon_naive - signal.sig;
    
    rel_err_naive = norm(Err_naive)/norm(signal.sig);
    
        % Play result!
    if playfiles == 1
         xyz = 'y';
         while xyz == 'y'
           wavplay(recon_naive,signal.fs);
            xyz = input('Replay naive inversion signal? y/n [y]: ', 's');
            if isempty(xyz)
                 xyz = 'y';
             end
         end
    end
    if inoutfiles == 1
        wavwrite(real(recon_naive),signal.fs,signal.nBits,sprintf('naiveinvsignal_%s',filename));
    end

    % apply iterative inversion of the multiplier 
    recon_inv=iframemul(filtered_signal,Fa,Fs,mask);
    recon_inv = recon_inv(1:Ls); % cut to original length
    Coeff_inv = frana(Fa,recon_inv);
    
     subplot(2,2,3);
     plotframe(Fa,Coeff_inv,signal.fs,'dynrange',dynrange);    
     set(gca,'YLim',[0 maxfreq])
     set(gca,'CLim',CCC);
     title('TF representation of iterative inversion');
    
    Err_inv = recon_inv - signal.sig;
    
    rel_err_inv = norm(Err_inv)/norm(signal.sig);
    
    subplot(2,2,2);
    plotframe(Fa,frana(Fa,Err_naive),signal.fs,'dynrange',dynrange);
    set(gca,'YLim',[0 maxfreq]);
    set(gca,'CLim',CCC);
    title(sprintf('TF representation of error of naive inversion \n rel. err.: %e',rel_err_naive));

    %     CCC2 = get(gca,'CLim');

    subplot(2,2,4);
    plotframe(Fa,frana(Fa,Err_inv),signal.fs,'dynrange',dynrange);
    set(gca,'YLim',[0 maxfreq]);
    set(gca,'CLim',CCC);
    
    title(sprintf('TF representation of error of iterative inversion \n rel. err.: %e',rel_err_inv));
      
    % Play result!
    if playfiles == 1
         xyz = 'y';
         while xyz == 'y'
           wavplay(recon_inv,signal.fs);
            xyz = input('Replay inversion signal? y/n [y]: ', 's');
            if isempty(xyz)
                 xyz = 'y';
             end
         end
    end
    if inoutfiles == 1
        wavwrite(real(recon_inv),signal.fs,signal.nBits,sprintf('invsignal_%s',filename));
        saveas(gcf,sprintf('inv_mult_%s.eps',filename),'eps');
        saveas(gcf,sprintf('inv_mult_%s.png',filename),'png');
        saveas(gcf,sprintf('inv_mult_%s.jpg',filename),'jpg');
        saveas(gcf,sprintf('inv_mult_%s.fig',filename),'fig');
    end
end

return


function xcm = nicecolormap()
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];


